\magnification=\magstep1
\def\p{\partial}
\centerline{Mathematics 304, section 1 -- Ordinary Differential Equations}
\centerline{Discussion/Lab 4 -- The Lorenz Equations}
\centerline{December 1-3, 2004}
\bigskip
\noindent
{\it Background}
\bigskip
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\'e in the
1880's on ODE problems in {\it celestial mechanics} already hinted
at what could happen in more dimensions.  However, Poincar\'e'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\'e's prodigious
intuition\footnote{${}^1$}{Students of mathematics
might recognize a kindred spirit in Poincar\'e -- 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\'e's, those
failings can be accepted.  To put that into perspective, though, Poincar\'e 
has been called the last person to understand {\it all}
of the mathematics known in his day, and there will
certainly never be another person like that again!})  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:
$$\cases{x' = f(x,y,z) = -10x + 10y &\cr
         y' = g(x,y,z) =  28x - y - xz & \cr
         z' = h(x,y,z) = {-8\over 3}z + xy&\cr}\leqno(LS)$$
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 {\it 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.

\bigskip
\noindent
{\it Discussion Questions}
\bigskip
\noindent
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,\quad g(x,y,z) = 0, \quad h(x,y,z) = 0.$$
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.)
\bigskip
\noindent
B)  Next we will study the {\it 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.
\item{1.}  First, compute the Jacobian matrix 
$$J = \pmatrix{ {\p f \over \p x} & {\p f \over \p y} & {\p f \over \p z} \cr
{\p g \over \p x} & {\p g \over \p y} & {\p g \over \p z} \cr 
{\p h \over \p x} & {\p h \over \p y} & {\p h \over \p z} \cr}$$
of the system (LS).  Then for each of your critical points,
substitute in the coordinates and determine the Jacobian there.
\item{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.)
\bigskip
\item{}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. 
\bigskip
\noindent
{\it Lab Questions}
\bigskip
\noindent
C)  For the two ``complicated'' critical points, we can use 
Maple to compute the eigenvalues and eigenvectors like this.
First, enter the command:
\bigskip
\centerline{\tt with(linalg);}
\bigskip
\noindent
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 = \pmatrix{1 & 2 & 3\cr
           4 & 5 & 6\cr
           7 & 8 & 9\cr}$$
for instance, use 
\bigskip
\centerline{\tt A := matrix([[1,2,3],[4,5,6],[7,8,9]]);}
\bigskip
\noindent
(matrices are entered as lists of lists, row-wise).  Then 
\bigskip
\centerline{\tt evalf(eigenvectors(A));}
\bigskip
\noindent
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 {\tt eigenvalues} command if the output is too cryptic.
What are the types of the other two critical points of the 
Lorenz system?
\bigskip
\noindent
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:
\bigskip
\centerline{\tt with(DEtools);}
\bigskip
\noindent
\item{1.}  Then define the Lorenz system as we did in the last lab (three 
equations now).  Call the system {\tt LEqns}.  Then enter
the following command to plot the solution with initial
conditions $x(0) = y(0) = z(0) = 1$, for $0 \le t \le 50$:
\bigskip
\noindent
\centerline{\tt DEplot3d(LEqns,[x(t),y(t),z(t)],t=0..50,[[x(0)=1,y(0)=1,z(0)=1]],}
\centerline{\tt stepsize=.01,linecolor=red,thickness=1);}
\bigskip
\item{}
(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.  
\item{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?
\bigskip
\noindent
E)  Now, let's focus on one coordinate function at a time.
\bigskip
\item{1.} We can plot $x(t)$ versus $t$ by using a different command
from the {\tt DEtools} package -- {\tt DEplot} (no ``3d'' this time)
command as follows:
\bigskip
\noindent
\centerline{\tt DEplot(LEqns,[x(t),y(t),z(t)],t=0..20,scene=[t,x(t)],}
\centerline{\tt [[x(0)=1,y(0)=1,z(0)=1]], stepsize=.01,linecolor=red,thickness=1);}
\bigskip
\item{}
The {\tt scene} option controls how the plot is generated.
Here we are plotting $x(t)$ versus the independent variable $t$.
\item{2.} Also plot $y(t)$ versus $t$ and $z(t)$ versus $t$.
\bigskip
\noindent
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):
\bigskip
\noindent
\centerline{\tt DEplot(LEqns,[x(t),y(t),z(t)],t=0..10,scene=[t,x(t)],}
\centerline{\tt [[x(0)=1,y(0)=1,z(0)=1],[x(0)=3,y(0)=1,z(0)=1]],}
\centerline{\tt stepsize=.01,linecolor=[red,black],thickness=1);}
\bigskip
\noindent
\item{1.} What happens if you increase the $t$ range to 
$0 \le t \le 20$ in the above?  What does that say about the
full solution (how does it relate to the full phase portrait)?
\item{2.}  Now try changing the initial value for $x$ in 
the second solution to $x(0) = 1.1$.  What happens now?
\item{3.}  What if $x(0) = 1.01, 1.001, 1.0001$, etc. (possibly
after extending the $t$ range)?
\item{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 ${\bf R}^3$ using 
the {\tt DEplot3d} command from question D above.  That is also
instructive!
\item{5.}  Explain why the phrase ``sensitive dependence on 
initial conditions'' has been applied to this system.  
\bigskip
\noindent
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?  
\bigskip
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 {\it Chaos,
Making a New Science} by James Gleick is an enjoyable,
general-audience introduction.
\bigskip
\noindent
{\it Assignment}
\bigskip
Group write-ups due no later than 5:00 pm on Tuesday, December 7.




\bye