Differential Equations Lab

This lab explores some of the ODE solving capabilities of Mathematica. The notebook is organized into several sections.  Each section is just a group of Mathematica cells containing explanatory text like this and actual commands like the one below. To execute a command, click the mouse in the cell containing the command then do a Shift-Return, that is, hold down the shift key and hit the Return key.  Try it on the cell below, which evaluates a complicated integral.

Integrate[(1 - 2*x^2)^(-3/2), x]

To begin the actual lab work you should open the first group of cells below by double-clicking on the bracket to the far right of the title Finding Exact Solutions with DSolve. The bracket will open up to reveal lots of interesting text and commands.  Just follow through the section, reading, executing commands and (hopefully) thinking about the results.  When you get done you can close up the section again by double clicking on the expanded bracket at the right edge of the cells and move on to the next section.  Have fun !

Finding Exact Solutions with DSolve

To find the exact formula for the solution of simple one-variable ODE's,  use the DSolve command.  Of course, not all equations can be solved.  When using this command, you have to be careful to use the double equal sign and to write all the y's as y[t] in the equation.

rules = DSolve[ y'[t] == (-1/2)*y[t] + Sin[t], y[t], t]

The answer is always a list of one or more subtitution rules. A list in Mathematica is always enclosed in set brackets like this { first thing, second thing ,... }.  In this case the list contains just one thing, a substitution rule.  A substitution rule is indicated by an arrow,   -> and it is also enclosed in set brackets.  Hence the double set brackets in the output.  The first (and in this case the only) element of the list of solution rules is called rules[[1]].

rules[[1]]

To extract the solution as a function for plotting or further analysis, you have to apply one of the substitution rules to the expression y[t].  You apply a rule in Mathematica using the cryptic  /.  notation.   Give the result another name, ygen,  so you can use it later.

ygen = y[t]/.rules[[1]]

Here ygen is the general solution with the arbitrary constant, C[1].  Note that this is a linear ODE and so the general solution consists of a particular solution plus and exponential term.To pick out a particular solution, just apply another substitution rule, substituting a particular constant for C[1].

ypart = ygen/.C[1]->0
Plot[ypart,{t,-10.0,10}];

To see lots of integral curves , you can plot a list of several solutions with different values of C[1].  One quick way to do this is to use the Table command to generate the list of particular solutions.  For some reason you need the Evaluate command inside the Plot command.

sollist = Table[ ygen/.C[1]->i/5,{i,-5,5}]

Plot[Evaluate[sollist],{t,-10,10}, PlotRange->{-2,2} ];

You can use the mouse to make the plot bigger.  Note that all solutions converge to the particular solution we plotted earlier.  This is due to the exponential decay of the term involving C[1].  In a case like this you would refer to the particular solution as the steady state solution.

Go back to the top of this section and try some other solvable ODE's.  For example, 
		y'[t] == y[t]^2	 or 	y'[t] == -t*y[t]
Just replace the differential equation inside DSolve with the new equation.  You may also want to change the plotting ranges.

Some other ODE's are too difficult to solve explicitly. Try the equation   y'[t] == (t^2-y[t]^2)/(1+y[t]^2).  The program just spits the problem back at you.

You can also use DSolve for second-order ODEs.  The first example is a typical damped, unforced oscillator.

rules = DSolve[ y''[t] + y'[t] + y[t]==0, y[t], t]

As for first-order equations, the answer is a list of one or more subtitution rules.  The first (and in this case the only) element of the list of solution rules is called rules[[1]].

rules[[1]]

To extract the solution as a function for plotting or further analysis, you have to apply one of the substitution rules to the expression y[t].  You apply a rule in Mathematica using the cryptic  /.  notation.   Give the result another name, ygen,  so you can use it later.

ygen = y[t]/.rules[[1]]

Here ygen is the general solution with two arbitrary constants, C[1] and C[2].  Note that the solution is given in terms of complex exponentials.  The exponents are cube roots of -1.  To get a real solution out of this you have to select the constants to be complex conjugates of one another.   Since this is rather inconvenient, we will adopt a different strategy, namely, specify some real-valued initial conditions directly in the DSolve command.

rules  =  
      DSolve[ {y''[t]  +  y'[t]  +  y[t] == 0, 
          y[0] == 1, y'[0] == 0},  y[t],  t]
y1  = y[t] /. rules[[1]]

This still doesn't look real but it actually is.  You can force Mathematica to put it into a nicer form by using the commands ComplexExpand and Simplify.  The former tries to expand out the expression into its real and imaginary parts under the assumption that all the variables in the expression are real.

y1  =  Simplify[ ComplexExpand[y1]]

Here is a plot.

Plot[y1, {t, 0, 10}];

Another independent real solution can be obtained by changing the initial conditions.

rules  =  
      DSolve[ {y''[t]  +  y'[t]  +  y[t] == 0, 
          y[0] == 0, y'[0] == 1},  y[t],  t]
y2  = Simplify[ ComplexExpand[y[t] /. rules[[1]] ]]
Plot[{y1, y2}, {t, 0, 10}];

Perhaps the best use of DSolve is to find explicit solutions of certain special second-order linear equations with non-constant coefficients.  In the 18th and 19th centuries, many such equations were solve via power series and their solutions were names after their discovers.  These are the so-called "special functions".  Mathematica has been trained to recognize many of these equations and can give you explicit solutions to them in terms of these special functions.  Here is Airy's equation   y''[t]+t*y[t]==0 and its solution in terms of "Airy functions".

rules  =  
      DSolve[ y''[t]  + t*y[t] == 0,  y[t],  t]

Here is what some Airy functions look like.  Again, the solution is complex so we first extract real and imaginary parts before plotting.

ygen  =   y[t] /. rules[[1]]
ypart = ygen /. {C[1] -> 1, C[2] -> 0}
Plot[{Re[ypart], Im[ypart]}, {t, -5, 5}];

Finding Approximate Solutions with NDSolve

For problems which cannot be solved explicitly you can still approximate the solutions numerically on any desired interval of time.  For this you can use the Mathematica function NDSolve.  The following command brings up information about this command.

?NDSolve

We will apply it to the ODE which DSolve failed to solve in the last section.

rules = NDSolve[
			{y'[t] == (t^2-y[t]^2)/(1+y[t]^2),y[0]==1},
			y[t],{t,-5,5}
			]

As before, you get a list of substitution rules, this time not given as a formula but as a so-called Interpolating Function. This is some function (whose formulas are hidden from you) which approximates the solution over the given time interval.  You plot it just as you did the particular solutions above.

ypart = y[t]/.rules[[1]];

Plot[ypart,{t,-5,5},PlotRange->All];

To draw lots of curves, make a list of solutions with different initial conditions.   This is a bit more complicated than just choosing different values for the arbitrary constant as we did above. We will start with an empty list and append solutions to it using the Append command  The solutions are obtained by repeating the process we just went through.  We will make the machine loop through these steps automatically using a Do loop which runs through the initial conditions y[0]=-3,-2.5,-2,...,0,...,2.5,3, that is, y[0]=i/2 for i=-6 to 6.

sollist = {};
Do[rules  = NDSolve[{y'[t]==(t^2-y[t]^2)/(1+y[t]^2),y[0]==i/2},y,{t,-5,5}];
AppendTo[sollist, y[t]/.rules[[1]]],
{i,-6,6}]

sollist

Finally we plot the solutions on the list.

Plot[ Evaluate[sollist], {t,-5,5}];

Now see if you can make a picture of the integral curves of the equation y'[t]==Sin[Pi*y[t]].  Just go back two cells and make some changes.  Can you explain the resulting picture by thinking about the phase line ?

2D Phase Portraits

NDSolve can be used to produce phase portraits for second-order ODEs or first-order systems.  Here is the second-order damped oscillator y '' + y ' + y = 0 written as a first-order system

\(rules\  = \ 
      NDSolve[{\(y'\)[t] == v[t], \(v'\)[t] == \ \(-v[t]\) - y[t], 
          y[0] == 1, v[0] == 0}, 
		{y[t], v[t]}, {t, 0, 10}]\)

To extract the solution from these substitution rules use the following

X  =  {y[t], v[t]} /. rules[[1]]

The answer is in the form of a vector X of two functions, representing the y and v coordinates of the solution. Now you can plot it using the function  ParametricPlot

ParametricPlot[X, {t, 0, 10}, AspectRatio -> 1];

To get a phase portrait you just have to compute several orbits and plot them together.  Before doing this it is convenient to automate the process of computing an orbit  a bit.   The following function   Orbit takes an ODE, an initial vector and a time as input and returns the correspoding solution ready for plotting.

\(Orbit[ode_, {y0_, v0_}, 
        time_]\  := 
	{y[t], 
          v[t]} /. \(NDSolve[
            ode~Join~{y[0] == y0, v[0] == v0}, {y[t], v[t]}, {t, 
              0, time}]\)[\([1]\)]\)

With this function defined, we can compute lots of orbits and plot them together to form the phase portrait.

ode  =  {y'[t] == 
            v[t], v'[t] ==  -v[t] - y[t]};
X1  =  Orbit[ode, {1, 1}, 10];
X2  =  Orbit[ode, {1, -1}, 10];
X3  =  Orbit[ode, {-1, -1}, 10];
\(\(X4\  = \ Orbit[ode, {\(-1\), 1}, 10];\)
\)
ParametricPlot[{X1, X2, X3, X4}, {t, 0, 10}, AspectRatio -> 1, 
        PlotRange -> All];

In the never ending quest for automation, we could go further and define a function which takes a list of initial conditions and returns a list of orbits.  This is pretty easy to do.  First,  the expression  Orbit[ode,#,time]&denotes a function of one argument (the initial condition) with the ode and the time fixed.  Then one can apply this new function to each entry in a list of initial conditions using the Map command as shown below.

OrbitList[ode_, iclist_, time_]  :=  
      Map[Orbit[ode, #, time] &, iclist]

We will try this out on a different linear system.  The iclist is a list of vectors representing the initial conditions of the orbits we want to compute.  For some reason you need the Evaluate command inside the plotting function.

ode  =  {y'[t] == 
            v[t], v'[t] == -4*y[t]};
iclist  =  {{1, 0}, {2, 0}, {3, 0}, {0.5, 0}};
\(\(orbs\  = \ OrbitList[ode, iclist, 10];\)
\)
\(\(ParametricPlot[Evaluate[orbs], {t, 0, 10}, 
	AspectRatio -> 1, 
        PlotRange -> {{\(-6\), 6}, {\(-6\), 6}}];\)\)

The Euler and Runge-Kutta Methods

Although the NDSolve command implements a complicated numerical method for you, it is easy and instructive to implement some of the simple methods discussed in class.  Below we define a function which carries out one step of Euler's method for any ODE.  As input it takes the function f[t,y] defining the ODE, the initial point {t,y} and the step size h.  As output it returns the point {t+h,y+ h*f[t,y]}.

EulerStep[f_,{t_,y_},h_] := {t+h,y+h*f[t,y]}

We can apply this method over and over again to step forward in time as far as we want.  At each step we want to start at the point obtained from the previous step.  The Mathematica function NestList automates this in a simple way.  It is designed to take a function and feed its outputs back into its inputs n times producing a list of all values obtained.  Throughout the process we want to keep the ODE f and the stepsize h fixed.  Only the point {t , y} is changing.  The expression EulerStep[f,#,h]& inside the NestList command represents a function with one argument (denoted by the # symbol) into which we repeatedly feed back values.

EulerList[f_,{t_,y_},h_,n_] := 
	NestList[ EulerStep[f,#,h]&,{t,y},n]

We will try it out on the simple ODE y ' = y with initial condition y[0]=1.  Of course we know the exact solution is the exponential function.  We will try to approximate it over the time interval from 0 to 1 by taking 4 steps of size h = 0.25.  The y value at the end of the list should be e = 2.71828 ...

f[t_,y_] := y
h = 0.25;
n = 4; 
elist = EulerList[f,{0,1},h,n]

To improve the approximation you need to make h smaller and n larger.  Go back to the last command and try taking  h=0.1 and n=10 and also h=0.01 and n=100.  

The lists are rather long and to continue, it is better not to print them all out.  You can simply print out the last element of the list using the Last command.  The command below shows the result of using 1000 steps.  Experiment to see how many steps it takes to get an approximation for  e  which is correct in the first four decimal places.

h = 0.001;
n = 1000;
Last[ EulerList[f,{0,1},h,n ] ]

Next we will compute an approximate solution to a more difficult problem and plot it.  Plotting just involves applying the Line command to the list produced by EulerList .  We will try to compute one of the solutions to the ODE y'[t]==(t^2-y[t]^2)/(1+y[t])^2 that we got in the previous section of this notebook using NDSolve . We will work on the time interval from -4 to 4.

f[t_,y_] := (t^2-y^2)/(1+y^2)
n = 4;
h = 8.0/n;
t0 = -4;
y0 = -2; 
elist = EulerList[f,{t0,y0},h,n];
Show[ Graphics[Line[elist]], Axes->True ];

You should go back to the last cell and try a sequence of larger n values so you can see the method converging to the actual solution.

For fun, you can make a movie showing the convergence of Euler's method by making a table of pictures and then using Animate Selected Graphics command in the Cell menu.  First we produce a table of 20 pictures with n values ranging from 4 to 80 and appropriately chosen h values.

\(Table[

      Show[Graphics[Line[EulerList[f, {t0, y0}, 8.0/\((4*n)\), 4*n]]], 
	
        Axes -> True, \ PlotRange -> {{\(-4\), 4}, {\(-2\), 4}}], 
{n, 1, 
        20}]\)

To animate them you have to select them by clicking with the mouse on the blue bracket on the right of the cell.  Once they are selected, go into the Cell menu and choose the item called Animate Selected Graphics.  You should see a movie of your 20 plots showing how Euler's method converges.  To stop the movie just click the mouse.  

If your movie goes too fast or too slow you can fix it as follows.  Go to the Preferences item (depending on your version of Mathematica, it may be in the Edit menu or somewhere else). Select GraphicsOptions and then Animation.  You will find an item called AnimationDisplayTime.  Click on the old value, enter a new one and hit Return.  Increasing it will slow the movie down.

We will now write our own version of the Runge-Kutta method. The only difference is that the method for taking one step is more complicated.  To carry it out we will use a Module command.  This allows us to use temporary variable to hold the four necessary slope values.

RKStep[f_,{t_,y_},h_] := Module[{m1,m2,m3,m4},
			m1 = f[t,y];
			m2 = f[t+h/2,y+h/2*m1];
			m3 = f[t+h/2,y+h/2*m2];
			m4 = f[t+h,y+h*m3];
			{t+h, y+(m1 + 2*m2 + 2*m3 + m4)*h/6.0 }
			]
			
RKList[f_,{t_,y_},h_,n_] := 
	NestList[ RKStep[f,#,h]&,{t,y},n]

This can be used in the same way as EulerList to produce an approximate solution,  The difference is that it does a better job for a given step size, h.  Here is a comparison of the two methods for the ODE we have been looking at with Euler in red and Runge-Kutta in green.  After you see the plot, try increasing n and watching them converge (the formula for h is adjusted automatically when you change n).

n = 4;
h = 8.0/n;

elist = EulerList[f,{t0,y0},h,n];
rklist = RKList[f,{t0,y0},h,n];

Show[ Graphics[{RGBColor[1,0,0],Line[elist],
		RGBColor[0,1,0],Line[rklist]}],Axes->True ];

Try increasing the number, n, of steps in the last cell from 4 to higher values.  For large n the two methods both give a good approximations to the actual solution.

Finally, we will compare the errors made by the two methods in approximating e by solving the differential equation  y ' = y with y[0] = 1 for the time interval from 0 to 1.  We will plot the differences between the values of y[1] produced by the two methods and the actual value of e  for various values of h.  To read off the value of y[1] we apply the Last method to get the pair {1 , y [1] } and then apply Last again.

e = N[Exp[1]];

f[t_,y_] := y

errortable = Table[
		{1.0/2^n,	
		e - Last[Last[ EulerList[f,{0,1},1.0/2^n,2^n]]],
		e - Last[Last[ RKList[f,{0,1},1.0/2^n,2^n]]]},
		{n,0,10}];

PrependTo[errortable,{"Stepsize","Euler","Runge-Kutta"}];

TableForm[errortable]

Recall that the three methods have errors of orders h and h^4 respectively.  This is reflected in the table where dividing h by 2  causes the errors to be divided by factors near  2 and 16 respectively.  To explore this we will make a table of the ratios of the errors as the step size is repeatedly cut in half.  We need to extract the second and third columns of the table above and make new list of ratios.

errortable  =  
        Drop[errortable, 1];

    eulererrors  =  Transpose[
            errortable][[2]];

    rkerrors  =  Transpose[errortable][[3]];

    eulererrorratios  =  
      Table[eulererrors[[i]]/eulererrors[[i + 1]], {i, 1, 
          Length[eulererrors] - 1}]

    rkerrorratios  =  
      Table[rkerrors[[i]]/rkerrors[[i + 1]], {i, 1, 
          Length[rkerrors] - 1}]

As you can see, the error for the Euler method changes by a factor of about 2 when the stepsize does whereas the Runge-Kutta errors change by a factor of 16.

We can also use the Euler and Runge-Kutta methods for first-order systems of differential equations.  The only change is that the function f[t , y] needs to be replaced by a vector valued function of a vector argument  f[t_, {y_,v_}] as shown below.  Here is the harmonic oscillator y'' + 2 y = 0 as a first order system.  Note that the initial condition is also a vector.

f[t_, {y_, v_}]  :=  {v,  -2*y}
elist = EulerList[f, {0, {1, 0}}, 0.1,  100]

To plot the phase portrait we need to extract the list of points without the time values.  The tricky way to do this is to take the second row of the transposed array.

ListPlot[Transpose[elist][[2]], PlotJoined -> True, 
      Axes -> True, AspectRatio -> 1, 
      PlotRange -> {{-3, 3}, {-3, 3}}]

The real solution moves on an ellipse,  Here we get a curve which appears to spiral out.  Try increasing the number of steps in the cell below to get a more accurate result.  You will find that it is rather hard to get something which looks right.

n = 100;
h  =  10.0/n;
elist = 
        EulerList[f, {0, {1, 0}}, h,  n];
ListPlot[Transpose[elist][[2]], PlotJoined -> True, 
      Axes -> True, AspectRatio -> 1, 
      PlotRange -> {{-3, 3}, {-3, 3}}]

Here is a similar test of the Runge-Kutta method for systems.  Try gradually increasing n until you get a reasonable picture.

n = 10;
h  =  10.0/n;
rklist = RKList[f, {0, {1, 0}}, h,  n];
ListPlot[Transpose[rklist][[2]], PlotJoined -> True, 
      Axes -> True, AspectRatio -> 1, 
      PlotRange -> {{-3, 3}, {-3, 3}}]

If you get to this point you could try replacing the harmonic oscillator by some other second-order ODE.  Or, for a challenge, try making a movie showing the convergence of the Euler method phase plots as n gets larger.