ME 406The Lorenz EquationssysidMathematica 6.0.3, DynPac 11.01, 1ê13ê2009intreset;plotreset; imsize = 250;‡1. IntroductionThis notebook contains all of the material given in class on the Lorenz equations, and it constitutessection 2.5 of the class notes. The Lorenz equations are given by(1)x° = s(y - x) , y° = rx - y -xz , z° = xy - bz .These equations contain three parameters: s, r and b. In what follows, we will always assume that theseparameters are positive. In all of the numerical calculations below, we take s = 10.0, and b = 8/3. These are themost used values in the study of the Lorenz equations. We will vary the parameter r over a wide range, andstudy how the solutions depend on r.These equations, which are simple in appearance, have solutions with extraordinary properties. Theywere first studied in the 1960's by the M.I.T. meteorologist Edward Lorenz. He developed the equations as amodel for the modal amplitudes in a nonlinear thermal convection problem. Lorenz recognized that the solutionsof the equations can exhibit an unusual form of behavior which we now call chaos. It took time for others torealize exactly what Lorenz had discovered. Lorenz has told the story of the discovery in his book The Essenceof Chaos, University of Washington Press, 1993. For a very readable and basic treatment of the equations, seeChapter 9 of Nonlinear Dynamics and Chaos, S.H. Strogatz, Addison-Wesley, 1994. For a general workcontaining a more advanced treatment, see Nonlinear Oscillations, Dynamical Systems, and Bifurcations ofVector Fields, J. Guckenheimer and P. Holmes, Springer-Verlag, 1983. For a book-length treatment containingmany detailed results, see The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors, C. Sparrow,Springer-Verlag, 1982. In our study here, we will continue our primarily experimental approach and use thecomputer to learn about the system. We begin by defining the equations for DynPac.setstate@8x, y, z<D; setparm@8s, r, b<D; sysname = "Lorenz";slopevec = 8s * Hy - xL, r * x - y - x * z, x * y - b * z<;‡2. Basic Propertiesü2.1 SymmetryThe equations are invariant to the transformation (2)(x, y, z) ö (-x, -y, z) .Thus if x(t), y(t), z(t) is a solution, so is -x(t), -y(t), z(t). As we shall see, this symmetry shows up in a number ofways, including the location of the equilibrium points of the system.ü2.2 The z-Axis is InvariantWe see from equations (1) that if x(0) = 0 and y(0) = 0, then x and y remain zero for all t. Thus the z-axisis an orbit, on which(3)z° = -bz , hence z(t) = z(0)e-bt, for x, y = 0. Thus the z-axis is always a part of the stable manifold for the equilibrium at the origin.ü2.3 The System is DissipativeThe divergence of the slopevector isdiv@slopevecD-1 - b - sand this is always negative. As we saw earlier in class, for any given volume V of phase points moving with theflow, we have (4)dVdt= V div[slopevec] = -(1 + b + s)V, hence V(t) = V(0)e-H1+b+sL t . With our canonical values of 10 for s and 8/3 for b, this is V(t) = V(0)e-13.67 t, so that volumes of initial pointsare reduced by a factor of e in a time 1 ê 13.670.0731529ü2.4 The Solutions are BoundedIt is not hard to prove that the solutions of the Lorenz equations are bounded. Consider first the casewhen r < 1. We examine a potential Liapunov function V = x2+ s * Iy2+ z2M;2 loreq.nbThe orbital derivative isSimplify@orbdt@VDD-2 Ix2- H1 + rL x y + y2+ b z2M sThe z-term is negative. We get for the x and y termsx2- H1 + rLx y + y2= Bx -12H1 + rL yF2+14 H1 - rL H3 + rL y2,and this quantity is postive except for x = y = 0. Thus V° < 0, and this shows that the origin is a global attractorfor the system when r < 1.When r > 1, we have to work a little harder to show the solutions are bounded. We start with a newLiapunov functionV = r * x2+ s * y2+ s * Hz - 2 * rL2;Level surfaces V = V0 are ellipsoids with center at (0, 0, 2r), and semi-axes (V0êr, V0ês, V0ês). Theorbital derivative isVdot = Simplify@orbdt@VDD-2 Iy2+ b z2+ r Ix2- 2 b zMM sWe work on this a bit.Vdotmod = Vdot ë I2 s b r2M-y2+ b z2+ r Ix2- 2 b zMb r2We can write this as(5)- x2br- y2b r2- Hz - rL2r2 + 1The expression in equation (5) is zero on the ellipsoid given by (6) x2br+ y2b r2+ Hz - rL2r2 = 1 ,and is positive inside that ellipsoid and negative outside that ellipsoid. By choosing V0 large enough, we canmake the level surface V = V0 lie entirely outside the ellipsoid (6). Any orbit starting outside that level surfacewill cross to the inside, and any orbit starting inside that level surface will remain inside. Thus all of the solu-tions are bounded, subject only to the assumption we made at the beginning that s, r, and b are positive. Unlikethe case r < 1, we cannot in this case draw any conclusions about a global attractor or even a stable equilibrium. loreq.nb 3ü2.5 EquilibriaWe use findpolyeq to look for equilibria.eqpoints = findpolyeq:80, 0, 0<, :- -b + b r , - -b + b r , -1 + r>,: -b + b r , -b + b r , -1 + r>>We see that the origin is an equilibrium for any values of the parameters. The other two equilibria are real if andonly if r ¥ 1. In fact for r = 1, the three equilibria coincide at {0,0,0}, and we therefore have a pitchfork bifurca-tion at r = 1. We will look at this in more detail later. We name the three equilibria for later reference.C0 = eqpoints[[1]]80, 0, 0<C1 = eqpoints[[2]]:- -b + b r , - -b + b r , -1 + r>C2 = eqpoints[[3]]: -b + b r , -b + b r , -1 + r>Note that the x and y components of C1 and C2 are the negatives of one another, consistent with the symmetrydiscussed earlier. In addition, for each equilibrium we have x = y, which is a direct consequence of x° = s(y - x)at equilibrium.‡3. Bifurcations with Changing rWe now begin the major task in this notebook, namely the study of the solutions as the parameter r ischanged.ü3.1 The Range 0 < r < 1.We look at the eigenvalues at the origin, which is the only equilibrium in this range of r-values.eigval@C0D:-b,12-1 - s - 1 - 2 s + 4 r s + s2,12-1 - s + 1 - 2 s + 4 r s + s2>The quantity under the square root sign can be written as (1 - sL2+4rs, which shows that the square root is real.The first and second eigenvalues are obviously negative. The condition for the third eigenvalue to be negative is(1 - 2s + 4rs