Mathematics 304, section 1 – Ordinary Differential Equations Discussion/Lab 4 – The Lorenz Equations December 1-3, 2004

Background

In the last 25 years or so, it has become widely recognized that the behavior of solutions of 3- and higher-dimensional systems of first order ODE can be much, much more interesting and complicated than systems in the plane. From the historical point of view, this probably should not have come as a surprise. Some work of the famous French mathematician and physicist Henri Poincaré in the 1880’s on ODE problems in celestial mechanics already hinted at what could happen in more dimensions. However, Poincaré’s work was “ahead of its time” in several ways and was not widely appreciated at the time or indeed for a long time after it appeared. Moreover, “ordinary” mathematicians (i.e. people without Poincaré’s prodigious intuition^{1}) did not have any good way to visualize these systems and understand just how strange they could be before the advent of modern computers and good software for 3D graphing.

In this final discussion/lab of the course, we will investigate one of the examples that alerted mathematicians to the wonderfully complex behavior possible in higher dimensional systems of ODE – the 3-dimensional Lorenz system:

(LS)  \{ \begin{matrix}x′ = f(x, y, z) = −10x + 10y \\ y′ = g(x, y, z) = 28x − y − xz \\ z′ = h(x, y, z) = \frac{−8}{3}z + xy \end{matrix}

This slightly unlikely-looking system of ODE was studied first by a meteorologist named Edward N. Lorenz (born in 1917 and still living) as part of an atmospheric model for weather forecasting. The exact interpretations of the variables and the meanings of the constant coefficients are too complicated for us to discuss. Suffice it to say that this system was arrived at by a sequence of simplifying assumptions and changes of variable rather typical of the process of modeling in applied mathematics.

Discussion Questions

A) As usual in our study of nonlinear systems, we will start by finding the equilibria of the system: the solutions of

f(x, y, z) = 0, g(x, y, z) = 0, h(x, y, z) = 0.

1 Students of mathematics might recognize a kindred spirit in Poincaré – he was amazingly good at seeing the “big picture” and identifying the key points about the problems he studied, but he was notoriously bad at developing complete, correct proofs of his assertions. His published papers are full of small mistakes and proofs that omit cases or have other gaps(!) If you have a mind as fertile and insightful as Poincaré’s, those failings can be accepted. To put that into perspective, though, Poincaré has been called the last person to understand all of the mathematics known in his day, and there will certainly never be another person like that again!

Where are the equilibrium points of the Lorenz system (LS)? (There are three in all; one is “obvious” two will require a bit of algebra to find.)

B) Next we will study the linearizations of the Lorenz system at the equilibrium points. We haven’t done this before for a 3-dimensional system, but the idea is exactly the same as in the 2-dimensional case.

1. First, compute the Jacobian matrix

J = ( \begin{matrix}\frac{∂f}{∂x} & \frac{∂f}{∂y} & \frac{∂f}{∂z} \\ \frac{∂g}{∂x} & \frac{∂g}{∂y} & \frac{∂g}{∂z} \\ \frac{∂h}{∂x} & \frac{∂h}{∂y} & \frac{∂h}{∂z} \end{matrix} )

of the system (LS). Then for each of your critical points, substitute in the coordinates and determine the Jacobian there.

2. There is one equilibrium point where the form of the Jacobian matrix is much simpler than the other two. (This should be immediately clear!) Find the eigenvalues and eigenvectors of the “simple one”. What does the phaseportait for the linearized system look like? (Describe in words or sketch.)

If you have some time left at this point, you can also look at the other equilibrium points, but the algebra of finding the eigenvalues and eigenvectors is much messier to do by hand.

Lab Questions

C) For the two “complicated” critical points, we can use Maple to compute the eigenvalues and eigenvectors like this. First, enter the command:

                               with(linalg);

to load the linear algebra package. Then you will need to enter the matrix for the linearized system. The format for entering a matrix looks like this. To enter

A = ( \begin{matrix}1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{matrix} )

for instance, use

                  A := matrix([[1,2,3],[4,5,6],[7,8,9]]);

(matrices are entered as lists of lists, row-wise). Then

                          evalf(eigenvectors(A));

will compute the eigenvalues and eigenvectors. See if you can interpret the output (Note: I is the imaginary unit i = \sqrt{−1} in Maple). See the online help page on the eigenvalues command if the output is too cryptic. What are the types of the other two critical points of the Lorenz system?

D) The linearizations do not tell the “whole story” about the Lorenz system, though. Next, let’s use Maple to plot its solutions and get a feeling for the global structure of the phase portrait. First enter:

                              with(DEtools);

1. Then define the Lorenz system as we did in the last lab (three equations now). Call the system LEqns. Then enter the following command to plot the solution with initial conditions x(0) = y(0) = z(0) = 1, for 0 ≤ t ≤ 50:

     DEplot3d(LEqns,[x(t),y(t),z(t)],t=0..50,[[x(0)=1,y(0)=1,z(0)=1]],
                 stepsize=.01,linecolor=red,thickness=1);

(plotting options chosen to make the structure of the solution more visible to start). The output here is a Maple 3D plot – you can rotate the axes and look at it from different viewpoints. Try to describe what the solution curve is doing as t increases. (You may want to change the interval of t-values, letting the endpoint gradually increase up to 50 to help you visualize this.) Also try a few different initial conditions and see if any different qualitative behaviors emerge.

2. Where are the three equilibrium points in this picture, and how do the types of the equilibria relate to the structure of the solutions you are observing?

E) Now, let’s focus on one coordinate function at a time.

1. We can plot x(t) versus t by using a different command from the DEtools package – DEplot (no “3d” this time) command as follows:

           DEplot(LEqns,[x(t),y(t),z(t)],t=0..20,scene=[t,x(t)],
    [[x(0)=1,y(0)=1,z(0)=1]], stepsize=.01,linecolor=red,thickness=1);

The scene option controls how the plot is generated. Here we are plotting x(t) versus the independent variable t.

2. Also plot y(t) versus t and z(t) versus t.

F) One of the most surprising things about the solutions of the Lorenz system is the following. Suppose we plot x(t) versus t for two solutions (one with x(0) = 1, y(0) = 1, z(0) = 1, the other with x(0) = 3, y(0) = 1, z(0) = 1) together, in different colors, using this kind of command (note the t-range):

           DEplot(LEqns,[x(t),y(t),z(t)],t=0..10,scene=[t,x(t)],
             [[x(0)=1,y(0)=1,z(0)=1],[x(0)=3,y(0)=1,z(0)=1]],
             stepsize=.01,linecolor=[red,black],thickness=1);

1. What happens if you increase the t range to 0 ≤ t ≤ 20 in the above? What does that say about the full solution (how does it relate to the full phase portrait)?

2. Now try changing the initial value for x in the second solution to x(0) = 1.1. What happens now?

3. What if x(0) = 1.01, 1.001, 1.0001, etc. (possibly after extending the t range)?

4. Experiment with the y and z coordinates of the solutions and see if the same sort of thing happens (possibly after extending the t range). You can also plot several solutions of the Lorenz system together in R^{3} using the DEplot3d command from question D above. That is also instructive!

5. Explain why the phrase “sensitive dependence on initial conditions” has been applied to this system.

G) (“Putting it all together”) Recall that the Lorenz system was supposed to be a model for some aspect of the atmosphere for weather forecasting. Supposing that its behavior reflects the properties of the real world (a fairly large assumption!), at least in part, what implications do you see? (Note: In weather forecasting, we would need to take actual measurements of conditions first, then plug that information into the model. Instruments for making those measurements have only limited precision and accuracy.) The solutions of the Lorenz system have sometimes been described as “chaotic”. Does that seem like a reasonable term?

If you want to look into some of the history of this system and other aspects of the “chaos theory,” that grew from the study of the Lorenz and other related systems, the book Chaos, Making a New Science by James Gleick is an enjoyable, general-audience introduction.

Assignment

Group write-ups due no later than 5:00 pm on Tuesday, December 7.