Math 260 Lab #13 Lorenz equations 
(adapted from lab in Interactive Differential Equations by West, Strogatz, McDill & Cantwell)

Save your work regularly--you never know when Mathematica may crash.  
Email your completed file to tleise@amherst.edu by Monday 4pm.

In 1963, Edward Lorenz published a paper on a simplified model for convection that has become a cornerstone for the field of nonlinear dynamics.  We will investigate changes in the behavior of solutions of the Lorenz equation as the parameter r is varied.  This parameter r>0 is called the Rayleigh number and is proportional to the difference in temperature from the warm base of a convection cell (for example, in the atmosphere) to the cooler top.  The variable x(t) is the speed of rotation of the cell.  The other two variables y(t) and z(t) represent temperature distribution in the cell.

dx/dt = 10(y-x)
dy/dt = rx-y-xz
dz/dt = xy-(8/3)z

Lorenz calculated a trajectory for this system using 6 decimal places of accuracy (computers were much more primitive in the 1960s), then reran it using three decimal places as a check.  But the two runs predicted very different trajectories.  Lorenz realized that the system exhibited sensitive dependence on initial conditions: very small changes in initial conditions can make very large differences in long-term  behavior.

Exercise 1.   The origin is a fixed point of the Lorenz equations. Find the eigenvalues of the linearization at the origin to show that the origin is an attractor if 0<r<1 but unstable if r>1.

Exercise 2.   Show that, in fact, for 0<r<1, the origin is a global attractor by showing that E(x,y,z)=x^2/10+y^2+z^2 is a Lyapunov function for the Lorenz system.  That is, show that dE/dt<0 if 0<r<1 for points other than the origin. (Hint: complete the square to get rid of the xy term.)

Exercise 3.   Show that the Lorenz equations are symmetric with respect to the z-axis: if (x(t),y(t),z(t)) is a solution, then so is (-x(t),-y(t),z(t)).

Exercise 4.   Find formulas for the other two fixed points in terms of the parameter r (you may use the Solve command).  For what values of r do these fixed points exist?  (These fixed points are called (SuperPlus[C]) and (SuperMinus[C]) in the text.)

Exercise 5.    The eigenvalues of the linearization at each of these other two fixed points are complicated functions of r, but we can still get some important information from them.  Fill in the linearization at the fixed point with positive coordinates in the eigenvalue command below.  The plot command then plots the real parts of the three eigenvalues as functions of r, so we can see whether the fixed point is stable or unstable for each value of r:

ev=Eigenvalues[A]
Plot[{Re[ev[[1]]],Re[ev[[2]]],Re[ev[[3]]]},{r,1,50},PlotStyle->{Black,Red,Blue}]

Exercise 6.  You should see that there is an eigenvalue whose real part is always negative and two eigenvalues whose real parts can be positive or negative, depending on the value of r.  At what value of r does the fixed point switch from stable to unstable?  Use FindRoot[Re[ev[[2]]]==0,{r,fill in guess based on plot}] to get a good estimate (or plot over a sufficiently small interval that you can estimate the zero-crossing point to two decimal places).

Exercise 7.  Plot a trajectory of the Lorenz system (like we did for the scroll circuit in Lab #11) for r=20, 24.06, 25, 99, 100, 101 over the time interval {t,0,100}, with initial point NOT at the origin (since it's a fixed point).  Note that not all large values of r yield chaotic solutions.  There are "periodic windows" with stable limit cycles, although nearby values of r yield chaotic solutions with "strange attractors".  Which of the above values of r corresponds to a periodic window?