Differential Equations Lab II

This lab explores more 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, shown in bolder type.  To execute a command, click the mouse in the cell containing the command then do a Shift-Enter (or Shift-Return), that is, hold down the shift key and hit the Enter key (or Return key).

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  Explicit Solutions for Second Order Equations. 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 !

Explicit Solutions for Second Order Equations

To find the exact formula for the solution of simple scalar 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.  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}];

And how about some Legendre functions   ?

\(n = 2; 

    rules = DSolve[\ 
        \((1 - t^2)\)*\(\(y'\)'\)[t]\  - 2*t*\(y'\)[t]\  + \ 
            n*\((n + 1)\) y[t] == 0, \ y[t], \ t]\)

This time both solutions are real. The first solution is a polynomial and the second is a special function.

ygen  =   y[t] /. rules[[1]]ypart1 = ygen /. {C[1] -> 1, C[2] -> 0}ypart2 = ygen /. {C[1] -> 0, C[2] -> 1}Plot[{ypart1, ypart2}, {t, -1, 1}];

Try going back and changing the value of n for the Legendre equation to some other integer.  If you try non-integers, there won't be a polynomial solution.

Eigenvalues, Eigenvectors -- Coupled Linear Oscillators

Mathematica has lots of built-in linear algebra functions which make it easy to work with linear systems.  A matrix is represented as a list of lists.  Each list represents one row of the matrix.  For example

A  =  {{1, 1}, {1, 0}}MatrixForm[A]

As you can see, the function  MatrixFormwill print the matrix in the familiar way.  Once you have entered the matrix you can find its determinant, characteristic polynomial, eigenvalues and eigenvectors.

Det[A]CharacteristicPolynomial[A, z]Eigenvalues[A]Eigenvectors[A]

Here is a more complicated example in four-dimensions.  Two coupled oscillators  are modelled by coupled second-order equations:  x1'' + omega1^2 x1 = eps(x2-x1) and x2'' + omega2^2 x2= eps(x1-x2).  Here x1 and x2 represent the positions of the oscillators.  If we introduce the velocities v1 and v2 as new variables then x1' = v1 and x2' = v2.  The second derivatives of the positions become the first derivatives of the velocities.  In the end we have a four-dimensional  first-order linear system with the matrix below.

\(A\  = 
      \ {{0, 1, 0, 0}, 
		{\(-omega1^2\) - eps, 0, eps, 0}, 

        		{0, 0, 0, 1}, 
		{eps, 0, \(-omega2^2\) - eps, 0}}; 

    MatrixForm[A]\)

You can compute the eigenvalues  even with the variable entries in the matrix.

CharacteristicPolynomial[A, z]Eigenvalues[A] // Simplify

This is fairly complicated so we will pick some particular values for the frequencies of the two oscillators and for the coupling constant.  First we will take the two frequencies equal.

\(omega1 = 1; 
omega2 = 1; 
eps\  = \ 0.1; 
MatrixForm[A]\)lambda  =  Eigenvalues[A]V  =  Eigenvectors[A]

All of the eigenvalues are imaginary numbers.  Each pair of imaginary eigenvalues determines a simple periodic motion of the coupled oscillator system.  These are the so-called normal modes of oscillation. You can find these periodic solutions as follows. First generate the complex solution corresponding to one of the complex eigenvectors, then extract its real part to get a real solution representing  a normal mode of oscillation.

Z[t]  =  Exp[lambda[[1]]*t]*V[[1]] // ComplexExpandX[t]  =  Re[Z[t]] // ComplexExpand

Now the four variables in our solution vector are {x1[t],v1[t],x2[t],v2[t]}.  We will plot the first and third components which represent the positions of the oscillators so we can visualize the normal mode. The first oscillator is shown in red and the second in green.

x1  =  X[t][[1]]x2 =   X[t][[3]]Plot[{x1, x2}, {t, 0, 10}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}];

As you can see, for this normal mode of oscillation, the two oscillators are exactly out of phase with one another.  When one moves to the left, the other moves to the right.  We can find the other normal mode in the same way using the other eigenvalue and corresponding eigenvector.

Z[t]  =  Exp[lambda[[3]]*t]*V[[3]] // ComplexExpandX[t]  =  Re[Z[t]] // ComplexExpandx1  =  X[t][[1]]x2 =   X[t][[3]]Plot[{x1, x2}, {t, 0, 10}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}];

For this normal mode the two oscillators are exactly in synch with one another.  These normal modes are the ones for the case when the two frequencies omega1 and omega2 are equal.  Now we are in a position to look at other case, too complicated to be worked out by hand.

In the cell below you can choose different values for the frequencies of the separate oscillators and then look at the two normal modes. For example omega1 = 1 and omega2 = 2,  the normal modes are entirely different than the ones above.  In one mode, the first oscillator moves with large amplitude while the second is almost motionless.  In the second mode, the situation is reveresed.  Try several other choices of omega1 and omega2.  For example, you could keep  omega1=1 and take a sequence of values for omega2 decreasing from 2 back down to 1.

\(omega1 = 1; 
omega2 = \ 2; 
eps\  = \ 0.1; 
MatrixForm[A]\)lambda  =  Eigenvalues[A]\(V\  = \ Eigenvectors[A]; 

    X1[t]\  = Re[\ Exp[lambda[\([1]\)]*t]*V[\([1]\)]] // ComplexExpand; 

    X2[t]\  = Re[\ Exp[lambda[\([3]\)]*t]*V[\([3]\)]] // ComplexExpand; 

    x1\  = \ \(X1[t]\)[\([1]\)]; x2 = \ \ \(X1[t]\)[\([3]\)]; 

    Plot[{x1, x2}, {t, 0, 10}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; 

    x1\  = \ \(X2[t]\)[\([1]\)]; x2 = \ \ \(X2[t]\)[\([3]\)]; 

    Plot[{x1, x2}, {t, 0, 10}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; \)

The normal modes are just the simplest kinds of solutions for the coupled oscillators. A typical solution looks like a linear combination of the two normal modes and is quasiperiodic rather than periodic.  To plot one of these solutions you just have to take a combination of the two complex solutions.

\(omega1 = 1; 
omega2 = 1; 
eps\  = \ 0.1; 
MatrixForm[A]\)lambda  =  Eigenvalues[A]\(V\  = \ Eigenvectors[A]; 

    Z[t]\  = \ 
      Exp[lambda[\([1]\)]*t]*V[\([1]\)]\  + 
        Exp[lambda[\([3]\)]*t]*V[\([3]\)]; 

    X[t]\  = Re[Z[t]] // ComplexExpand; 
x1\  = \ \(X[t]\)[\([1]\)]; 
    x2 = \ \ \(X[t]\)[\([3]\)]; 

    Plot[{x1, x2}, {t, 0, 60}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; \)

For the case of equal frequencies shown above, the resulting motion is very interesting.  First one oscillator has large amplitude while the other is relatively quiet, then later the situation is reversed.

Linear Flows -- Mapping the Cat

For a linear system X' = AX we have seen that the solution is a linear function of its ititial condition, namely,
X(t) = exp(At) X(0) where exp(At) is the matrix exponential.  Thinking of all initial vectors X(0) moving simultaneously in this way we have the concept of the flow of the system. In this section we will compute exp(tA) for some simple linear systems in the plane and use it to animate the flow.

There is a built-in  Mathematica function called  MatrixExp which can actually compute the  flow matrix in simple cases.
For example,

\(A\  = \ {{1, 1}, {1, 0}}; 
MatrixForm[A]\)\(etA\  = \ MatrixExp[t*A] // Simplify; 
MatrixForm[etA]\)

Next we will use the flow matrix to transform a simple shape.  The simplest shapes are curves represented by   Mathematica Line objects.  A Line is just a list of points to be connected with line segments.  For example, a square could be represented

square  =   Line[{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}}]Show[Graphics[square], Axes -> True, AspectRatio -> 1];

To transform the square we just need to multiply each of the points in the Line by the flow matrix.  We will define a function to do this for us.

Transform[l_Line, A_]  :=  Line[Map[A . #&, l[[1]]]]

To compute the image of the square under the time 1 map of the flow, set

\(t = 1; 
sq1\  = \ Transform[square, etA]; 

    Show[Graphics[{square, sq1}], Axes -> True]; \)

Of course an animation is even better.

Do[
      Show[Graphics[Transform[square, etA]], Axes -> True, 
        PlotRange -> {{0, 6}, {0, 5}}], {t, 0, 1, 0.1}];

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 plots showing how the square is transformed.  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 in the Edit menu. Select GraphicsOptions and then  Animation.  You will find an item called AnimationDisplayTime.  Click on the old value, enter a new one and hit Enter (or Return).  Increasing it will slow the movie down.

It would be much more interesting to animate the images of some more complicated shape.  The shape below consists of several Line objects which together form an approximation to a cat face.

\(c1\  = \ Line[Table[{2*Cos[t], 1.5*Sin[t]}, {t, 0, 2*Pi, Pi/20}]]; 

    c2\  = \ Line[
        Table[{1 + 0.3*Cos[t], 0.5 + 0.2*Sin[t]}, {t, 0, 2*Pi, Pi/20}]]; 

    c3\  = \ Line[
        Table[{\(-1\) + 0.3*Cos[t], 0.5 + 0.2*Sin[t]}, {t, 0, 2*Pi, Pi/20}]]; 
    
c4\  = \ 
      Line[{{2*Cos[Pi/6], 1.5*Sin[Pi/6]}, {1.5, 2}, {2*Cos[Pi/3], 
            1.5*Sin[Pi/3]}}]; 

    c5\  = \ Line[{{2*Cos[5*Pi/6], 1.5*Sin[5*Pi/6]}, {\(-1.5\), 2}, {
            2*Cos[2*Pi/3], 1.5*Sin[2*Pi/3]}}]; 

    c6\  = \ Line[{{2.2, \(-0.8\)}, {1, \(-0.6\)}, {2.5, \(-0.5\)}}]\)\(c7\  = \ 
      Line[{{\(-2.2\), \(-0.8\)}, {\(-1\), \(-0.6\)}, {\(-2.5\), \(-0.5\)}}]

    \)\(cat\  = \ {c1, c2, c3, c4, c5, c6, c7}; 


    Show[Graphics[cat], PlotRange -> {{\(-3\), 3}, {\(-3\), 3}}, 
      Axes -> True, AspectRatio -> 1]; \)

Since the cat is made of a list of lines, we need to define a function for tranforming lists of lines by transforming each line in the list.

Transform[l_List, A_]  :=  Map[Transform[#, A]&, l]\(t = 1; 
cat1 = \ Transform[cat, etA]; 

    Show[Graphics[cat1], Axes -> True]; \)

And finally,  CAT -- The Movie.

\(t =. ; 
A\  = \ {{1, 1}, {1, 0}}; 
MatrixForm[A]\)\(etA\  = \ MatrixExp[t*A] // Simplify; 


    Do[Show[Graphics[Transform[cat, etA]], Axes -> True, AspectRatio -> 1, 
        PlotRange -> {{\(-10\), 10}, {\(-10\), 10}}], {t, 0, 1, 0.1}]; \)

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 plots showing how the cat is transformed.  To stop the cat just click the mouse.  

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

You might want to go back and try mapping the cat with some other matrices A which have different kinds of flows.   For example A = {{1,-5},{5,1}} has a unstable spiral at the origin while A = {{2,0},{0,1}} has an unstable node.

2D Phase Portraits

For nonlinear systems, it is not generally possible to find explicit formulas for the solutions. However,   NDSolve can be used to generate an approximate solution corresponding t o a given initial condition.  You have to specify the ODE and the initial conditions as indicated below. This example is a simple linear system in the plane.

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

To extract the solution from these substitution rules use the following

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

The answer is in the form of a vector X of two functions, representing the x and y 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_, {x0_, y0_}, time_]\  := 

      	{x[t], y[t]} /. 
        \(NDSolve[
            ode~Join~{x[0] == x0, y[0] == y0}, {x[t], y[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\  = 
      \ {\(x'\)[t] == \(-x[t]\)/2 - y[t], \(y'\)[t] == \ x[t] - y[t]/2}; 

    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\  = \ {\(x'\)[t] == y[t], \(y'\)[t] == \(-4\)*x[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}}]; \)

Now for some nonlinear systems.  Here is a competing species system.  By judiciously choosing the intial conditions, you can get a pretty comprehensive phase portrait.   The choice of a and b below leads to coexistence of the two species.  Try some other values such as:   a=0.8,b=0.6 and   a=1.5,b=1.2

a=1.2;
b=1.5;

ode = {x'[t]==x[t]*(1-x[t]-y[t]),y'[t]==y[t]*(a-x[t]-b*y[t])}

iclist = {{1.5,0.1},{1.5,0.5},{1.5,1},{1.5,1.5},
		  {1,1.5},{0.5,1.5},{0.1,1.5},
		  {0.01,0.1},{0.05,0.1},{0.1,0.1},
		  {0.1,0.05},{0.1,0.01}};
orbs = OrbitList[ode,iclist,20];

ParametricPlot[Evaluate[orbs],{t,0,20},
	Axes->True,AspectRatio->1,
	PlotRange->{{0,1.5},{0,1.5}}];

Now for the famous Van der Pol equations. There is a single attracting periodic orbit to which most other orbits converge.  Try a few other values of  mu, such as  mu = 0.1, 1, 10

mu = 0.6;

ode := {x'[t]==y[t],y'[t]==-x[t]+mu*(1-x[t]^2)*y[t]}

iclist = {{4,4},{4,-4},{-4,4},{-4,-4},{0.1,0.1}};
orbs = OrbitList[ode,iclist,20];

ParametricPlot[Evaluate[orbs],{t,0,20},
	Axes->True,AspectRatio->1,PlotRange->All];

The nonlinear oscillation of the Van der Pol equation is far from the usual sinusoidal oscillations produced by linear systems.  To illustrate this it is better to look at the graphs of the solutions rather than their orbits.  We can use our   Orbit function to compute a solution and then plot the x and y components of this solution as functions of t.

\(mu\  = \ 10; 
{x, y}\  = \ Orbit[ode, {2, 0}, 20]; 

    Plot[{x, y}, {t, 0, 20}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; \)

Try several different  mu values in the cell above.  For   mu=10 the oscillation looks more like a step function than a sine function.

3D Phase Portraits -- Lorenz and Rossler Attractors

To handle 3D nonlinear systems,  only a small changeto the functions defined in the last section is required.  We need a version of the Orbit function which takes a 3D initial vector.

\(Orbit[ode_, {x0_, y0_, z0_}, time_, opts___]\  := 

      	{x[t], y[t], z[t]} /. 
        \(NDSolve[
            ode~Join~{x[0] == x0, y[0] == y0, z[0] == z0}, {x[t], y[t], 
              z[t]}, {t, 0, time}, opts]\)[\([1]\)]\)

Also, when plotting the results you have to use  ParametricPlot3D. As a simple first example, we will try a linear equation with two spiralling eigenvalue and one real one.

ode = {	x'[t]==-0.2*x[t]-y[t],
		y'[t]==x[t]-0.2*y[t],
		z'[t]==-0.1*z[t]};

iclist = { {1,1,0},{1,1,1},{-1,-1,-1} };
orbs = OrbitList[ode,iclist,30];

pp = ParametricPlot3D[Evaluate[orbs],{t,0,30},
	BoxRatios->{1,1,1},PlotRange->All];

This shows three orbits spiralling in toward the origin, one in the (x,y) plane, one above it and one below it.  You can look at this from a different point of view in 3D by using the viewpoint option.  Since  the pictureis named pp, you can just use the Show command to redisplay it.

Show[pp,ViewPoint->{1,2,0}];
Show[pp,ViewPoint->{0,0,1}];

In 3D it is possible for ODE's to have attractors which are more complicated than restpoints or limit cycles.  Probably the most famous ODE with a "strange" attractor is the Lorenz equation.  One  long orbit of this is sufficient since all orbits end up on the attractor anyway.

Unprotect[Out];Clear[Out];  (*this is to free up memory*)

ode = {	x'[t]== 10*(y[t]-x[t]),
		y'[t]== -x[t]*z[t]+28*x[t]-y[t],
		z'[t]== x[t]*y[t]-(8/3)*z[t]};
orb = Orbit[ode,{1,0,10},50,MaxSteps->5000];

ParametricPlot3D[orb,{t,0,50},PlotPoints->5000];

Next we will make a movie by generating a sequence of pictures of longer and longer segments of the orbit. This requires lots of memory so to prevent crashes there is some code in the loop to remove old output variables.

range = {{-30,30},{-30,30},{0,50}};
Do[	ParametricPlot3D[orb,{t,0,5*i},
AspectRatio->1,PlotRange->range,PlotPoints->3000,Compiled->False];Clear[Out];,
{i,1,10}];

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 plots showing how the orbit gradually fills in the attractor.  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 in the Edit menu. Select GraphicsOptions and then  Animation.  You will find an item called AnimationDisplayTime.  Click on the old value, enter a new one and hit Enter (or Return).  Increasing it will slow the movie down.

Another three-dimensional system with a strange attractor is the Rossler system.  This attractor looks like a fractal Mobius band.

ode = {	x'[t]== -y[t]-z[t],
		y'[t]== x[t]+ 0.398*y[t],
		z'[t]== 2 + z[t]*(x[t]-4)};
orb = Orbit[ode,{1,0,0},200,MaxSteps->5000];

ParametricPlot3D[orb,{t,0,200},PlotPoints->2000,PlotRange->All];