Ordinary Differential Equations

George E. Hrabovsky

MAST

Introduction

I trust that  everyone taking this course knows what a ordinary differential equation (ODE) is. This is not the place for an introduction. I am going to proceed on the assumption that you have all encountered ordinary differential equations and have solved them in the past. What this lesson will cover is a systematic discussion of how to solve these equations using Mathematica. The chief advantages of using Mathematica are that you can solve a huge number of ODEs symbolically, you can merge numerical and symbolic solution methods, and the development time for coding is fractional as compared to other systems. If you need a quick refresher on the topic of ODEs, there is usually a pretty good survey chapter in any good mathematical methods book.

DSolve

Mathematica has a command that is designed to provide a symbolic solution to a differential equation. It can be used for ODEs and partial differential equations (PDEs), integral equations, and differential-algebraic equations (DAEs). Say we have an ordinary differential equation,

l3_1.png

(3.1)

Solving this by hand requires us to recognize that the equation satisfies a differential operator of the form

l3_2.png

(3.2)

We next need to recognize that this operator has equal powers of x and orders of D in each term. We need to think  of an ansatz to continue, we choose l3_3.png, where we rewrite the equation,

l3_4.png

(3.3)

This implies

l3_5.png

(3.4)

We solve for a and find that a=(1,2). Thus the general solution is l3_6.png. Let’s see how DSolve does.

l3_7.png

l3_8.png

l3_9.png

l3_10.png

How do we know if this is a solution? We can plug it into the original equation and make sure it is satisfied. To do this construct each of the terms using the solution.

l3_11.png

l3_12.png

If we solved this using DSolveValue, we would not need the part, [[]], specification.

l3_13.png

l3_14.png

l3_15.png

l3_16.png

We then construct the equation.

l3_17.png

l3_18.png

So this solution satisfies (3.1).

Linear ODEs

One class of common ODEs is the linear ODE. DSolve handles these quite well. For example, here is an exponential ODE.

l3_19.png

l3_20.png

If we specify an initial-value of y this becomes an initial-value problem (IVP).

l3_21.png

l3_22.png

Inhomogeneous equations utilize the Method of Variation of Parameters or the Method of Undetermined Coefficients by default (Mathematica chooses which one to use).

l3_23.png

l3_24.png

Note the two parts of the solution. The second part is the general solution of the homogeneous equation.

l3_25.png

l3_26.png

The first part of the result in sol3 is a particular solution of the inhomogeneous equation. We can get a list of these by varying the parameters.

l3_27.png

l3_28.png
l3_29.png
l3_30.png
l3_31.png
l3_32.png
l3_33.png
l3_34.png
l3_35.png
l3_36.png
l3_37.png
l3_38.png
l3_39.png
l3_40.png
l3_41.png
l3_42.png
l3_43.png
l3_44.png
l3_45.png

We can plot these.

l3_46.png

l3_47.gif

As an exercise read the documentation for DSolve noting specifically Basic Examples and both Basic Uses and Linear Differential Equations under Scope. Work your own examples based on your reading.

Read the tutorial Symbolic Differential Equation Solving with specific reference to the sections Introduction to Differential Equation Solving with DSolve, Classification of Differential Equations, and under ODEs the subsections Overview of ODEs, First-Order ODEs, Linear Second-Order ODEs, and Higher-Order ODEs. Here you will encounter most of the kinds of linear ODEs encountered in a standard course.

Nonlinear ODEs

Mathematica is also capable of solving nonlinear ODEs. We begin with a classic Riccati equation.

l3_48.png

Its solution is then.

l3_49.png

l3_50.png

Here is a second-order nonlinear ODE.

l3_51.png

l3_52.png

Sometimes the solution requires inverse functions.

l3_53.png

l3_54.png

l3_55.png

We begin by adding a set of conditions that constitute a solution for DSolve, in other words we convert a solution from DSolve into a set of conditions. This will give us a result similar to Reduce.

l3_56.png

Here we attempt to use the command Internal`InheritedBlock.

l3_57.png

l3_58.png

There is no advice on what it is or how to use it in the system. It is similar to Block, however one difference is that the function itself is changed locally, but not globally. Here, inside of our InheritedBlock we unprotect Solve and change it so that it uses Reduce. DSolve internally calls to Solve, so this may be a way to allow us to get results similar to Reduce.

We will also use the undocumented Internal`WithLocalSettings. This command allows to perform a cleanup so that things do not become jumbled. Specifically this applies preprocessing, then code, and then postprocessing where the postprocessing takes place before returning from aborts or similar jumps.

This method of solving the ODE can take a very long time, if it works. Be sure to copy your ODE into the expression sol= below.

l3_59.gif

l3_60.png

l3_61.png

l3_62.png

l3_63.png

l3_64.png

l3_65.png

l3_66.png

l3_67.png

l3_68.png

l3_69.png

l3_70.png

l3_71.png

l3_72.png

l3_73.png

l3_74.png

l3_75.png

l3_76.png

l3_77.png

l3_78.png

l3_79.png

l3_80.png

l3_81.png

l3_82.png

l3_83.png

These are not error messages or warnings, they just tell you that something is on or off.

l3_84.png

l3_85.png

l3_86.png

l3_87.png

l3_88.png

l3_89.png

l3_90.png

l3_91.png

l3_92.png

l3_93.png

l3_94.png

l3_95.png

l3_96.png

l3_97.png

l3_98.png

l3_99.png

l3_100.png

l3_101.png

l3_102.png

l3_103.png

l3_104.png

This code is adapted from the work of the user Michael E2 on the Mathematica Stack Exchange.

As an exercise read the documentation for DSolve noting specifically nonlinear Differential Equations under Scope. Work your own examples based on your reading.

Read the tutorial Symbolic Differential Equation Solving with specific reference to the subsections under ODEs, Nonlinear Second-Order ODEs.

Systems of ODEs

We begin with a simple diagonal system

l3_105.png

l3_106.png

We can solve an inhomogeneous system having constant coefficients.

l3_107.png

l3_108.png

l3_109.png

Here we have a linear system having rational coefficients.

l3_110.png

l3_111.png

Here we have a nonlinear system.

l3_112.png

l3_113.png

l3_114.png

We can also perform DSolve on vector equations. Here we define a matrix and a vector function of a variable.

l3_115.png

l3_116.png

l3_117.png

Not every element of the system needs to be a differential equation.

l3_118.png

l3_119.png

Such a system is called a differential-algebraic equation (DAE).

As an exercise read the documentation for DSolve noting specifically Systems of Differential Equations and Differential-Algebraic Equations under Scope. Work your own examples based on your reading.

Read the tutorial Symbolic Differential Equation Solving with specific reference to the subsection under ODEs, Systems of ODEs. Then read the section on DAEs.

Using Symbolic Vectors and Matrices

Beginning with version 9 in 2012, a new way was introduced to abstractly represent of vectors in a specified number of dimensions. This is the Vectors command.

l3_120.png

l3_121.png

We can thus declare a symbol, or collection of symbols, to be a vector.

l3_122.png

l3_123.png

Note the assumption that the domain is the complex numbers. We can, assign this, too.

l3_124.png

l3_125.png

Beginning in version 13 a new type of specification in ODEs can be made. We can state that a dependent variable is an abstract vector. For example, say in sol10 above we did not specify the components of the vector V[t]. We can rename it v[t] here.

l3_126.png

l3_127.png

We can use DSolveValue to get the formulas without the transformation rules.

l3_128.png

l3_129.png

Delay Equations

A delay-type differential equation is one where the time-derivative at some time is written in terms of the function at earlier times. In general, if our function is l3_130.png where we define l3_131.png, then the equation has the form,

l3_132.png

(3.5)

We can specify this condition using the condition symbol /;.

l3_133.png

1
l3_134.png
2
l3_135.png
3
l3_136.png
4
l3_137.png
5
l3_138.png

As an exercise read the documentation for DSolve noting specifically Delay Differential Equations under Scope. Work your own examples based on your reading.

Read the tutorial Symbolic Differential Equation Solving with specific reference to the subsection under Working with DSolve. Play with some of your own examples.

NDSolve

It is inevitable that you will encounter a situation where you cannot solve an equation analytically. In such cases you will need to solve the differential equation numerically. It is used in a manner very much like DSolve, except you must specify all parameters, constants, initial-values, and boundary values. The result of any NDSolve is an interpolating function that captures the behavior of the numerical solution. It is important to realize that you should not extend the interpolating function beyond the parameters of the case that you solved. For example, if you numerically solve an equation for values of t from 0 to 25, do not expect the interpolating function to be accurate if you set t to be 100.

For example, say we have the equation and its initial condition to solve.

l3_139.png

l3_140.gif

We can plot this equation.

l3_141.png

l3_142.gif

We can also plot its phase portrait.

l3_143.png

l3_144.gif

Of course, if we solve a second-order equation we need a system of two initial values.

l3_145.png

l3_146.gif

Here we see this plot.

l3_147.png

l3_148.gif

We can even zoom in, if we suspect small oscillations.

l3_149.png

l3_150.gif

We can even apply it to systems of ODEs.

l3_151.png

l3_152.gif

We can plot the dynamics of this system.

l3_153.png

l3_154.gif

Read the documentation for the NDSolve command, specifically Basic Examples and Ordinary Differential Equations under Scope. Develop some examples and play with them.

Read the Introduction to the tutorial Advanced Numerical Differential Equation Solving in the Wolfram Language, make up some examples and play with them.

NDSolve Methods for ODEs

Unless you state otherwise, Mathematica will choose what it considers to be the appropriate solution method for your problem, but it will not tell you what it did. As with many of the commands in Mathematica, you can choose the method to use. Unlike many, there is extensive documentation of these methods in the Advanced Numerical Differential Equation Solving in the Wolfram Language tutorial. I will show off a few of these methods below.

Every course in numerical methods talks about the Euler method. It was one of the first methods and it is easy to implement. It is not the most accurate, and as we will see, it can lead to instabilities. It is an explicit method to integrate an evolution equation, thus we call it a time integration scheme. What we are doing is discretizing the time into a sequence of points and integrating at each point. Here we apply the Euler method to the harmonic oscillator problem.

l3_155.png

l3_156.gif

We can plot the InterpolatingFunction.

l3_157.png

l3_158.gif

One way to improve accuracy is to reduce the initial step size.

l3_159.png

l3_160.gif

l3_161.png

l3_162.gif

Note that this causes an increase in the amplitude of oscillation. We can see this on a phase plot.

l3_163.png

l3_164.gif

This is what I meant earlier about instability. We can improve this by using the explicit midpoint method, often called the improved Euler method.

l3_165.png

l3_166.png

l3_167.gif

We have gotten a warning. It seems that we need to increase the number of steps to cover our domain.

l3_168.png

l3_169.gif

l3_170.png

l3_171.gif

Within the domain of interest the stability problem seems to have been solved. Here is the phase plot.

l3_172.png

l3_173.gif

There are many more accurate methods. Here is an explicit Runge-Kutta method.

l3_174.png

l3_175.gif

Here is the plot.

l3_176.png

l3_177.gif

Read the documentation for NDSolve with particular reference to Details and Options, and Options>Methods. Work up some examples of your own and play with them.

Read the section of the tutorial Advanced Numerical Differential Equation Solving in the Wolfram Language ODE Integration Methods and play with some examples.

Stability and Stiffness

We have mentioned stability and even seen an example of instability. What follows is based on the online documentation.

What is stability? Let's say that we have an initial-value problem characterized by the equation

l3_178.png

(3.6)

with the conditions l3_179.png and l3_180.png. We then apply the Dahlquist test equation

l3_181.png

(3.7)

where λ∈C and Re(λ)<0. This is named after the Swedish mathematician Germund Dahlquist (1925-2005), who made significant progress in the theory of numerical analysis. We can write our numerical solution in the form

l3_182.png

(3.8)

where, for some step size h,

l3_183.png

(3.9)

and r(z) is a rational function representing stability. That being the case there is a region,

l3_184.png

(3.10)

that we can call the boundary of absolute stability.

For a more concrete example we can look at the explicit Euler method described by,

l3_185.png

(3.11)

Applying this to (3.8) we can see that

l3_186.png

(3.12)

To continue with this we need to load the Function Approximations package.

l3_187.png

We then examine an OrderStarPlot.

l3_188.png

l3_189.png

l3_190.png

l3_191.gif

The shaded region represents instability. If we look to where the boundary of the white region intersections the negative axis, we call this the linear stability boundary (LSB). For this method we see this value is -2.

For the eigenvalue λ=-1 we need a step size such that h<2. Not hard to do. On the other hand, say we want l3_192.png, then we require a step size where l3_193.png, very tough to do.

When you attempt to numerically solve a problem that has instabilities, we call that stiffness. Even if the point you are concerned with might be stable, if there are nearby points that are not then the problem is stiff.  What is the problem? Stiff problems take longer for non-stiff methods to compute. It is an issue of the efficiency of computation. If you do not care how long your system takes to solve the problem, then you don’t have to worry about it. On the other hand, if time is a concern (say you are trying to develop a numerical forecast of the weather before the tornado occurs) then you will need some sort of stiffness handling. Most explicit methods in Mathematica have the ability to perform a StiffnessTest to allow for this possibility.

To illustrate, we return to an example from above.

l3_194.png

l3_195.gif

l3_196.png

l3_197.gif

Note that this has taken care of the instability of the method.

l3_198.png

l3_199.gif

Read the section of the tutorial Advanced Numerical Differential Equation Solving in the Wolfram Language ODE Integration Methods and play with some examples.

There is an extensive discussion of this on the Mathematica Stack Exchange at the website:

https://mathematica.stackexchange.com/questions/39028/understanding-ndsolvendsz/237989#237989

Events

Sometimes we need to consider changes to the state of a system that are conditional in some way. For example, if we drop a ball and it hits the ground, it most likely will bounce. How do we capture that in a model? We use the WhenEvent[] command. For example, if we start with the equation of motion under the action of the acceleration due to gravity we write it this way, with initial conditions. We also assume that the ball bouncing absorbs 20% of the force with each bounce.

l3_200.png

We can plot this

l3_201.png

l3_202.gif

Go ahead and read the documentation for the WhenEvent command. Pay particular attention to the Details and Options, and then to Applications. Make up some examples and play with them.

A More Traditional Numerical Approach

There are those people who are, for one reason or another, wanting to do things the old fashioned way. In a later lesson we will explore more deeply the ideas of finite differences. Here we will give an overview of how to duplicate the results of older, tried and true, methods. We begin with a differential equation and an initial condition.

l3_203.png

l3_204.png

We apply the DifferenceQuotient command to make this into a difference equation.

l3_205.png

l3_206.png

We have already seen the solution of this in ivp1. We can use DSolveValue to get the solution as a formula.

l3_207.png

l3_208.png

We use RSolveValue to get the formula to solve the difference equation.

l3_209.png

l3_210.png

We can compare the plots of these for any case where we assign a value of a.

l3_211.png

l3_212.gif

We will explore this more in a later lesson. I will note that as of version 13 RSolve is now able to solve all linear systems with polynomial coefficients.

Created with the Wolfram Language