The Lorenz model Based in part on a Maple worksheet by H.J. L\303\274dde Modified by J. Reinhardt, 2007, 2009 We study the three coupled differential equations set up by E.N. Lorenz in 1963: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYsLUkmbWZyYWNHRiQ2KC1JI21vR0YkNi1RMCZEaWZmZXJlbnRpYWxEO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R0Y3LyUqc3ltbWV0cmljR0Y3LyUobGFyZ2VvcEdGNy8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZGLUYjNiZGLi1GLzYtUSJ+RidGMi9GNlEmZmFsc2VGJy9GOUZPL0Y7Rk8vRj1GTy9GP0ZPL0ZBRk8vRkNGT0ZERkctSSNtaUdGJDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvRjNRJ2l0YWxpY0YnRjIvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRl5vLyUpYmV2ZWxsZWRHRk9GSy1GVzYlUSJ4RidGWkZnbkZLLUYvNi1RIj1GJ0YyRk5GUEZRRlJGU0ZURlUvRkVRLDAuMjc3Nzc3OGVtRicvRkhGam9GSy1GVzYlUScmIzk2MztGJy9GZW5GT0YyRkstSShtZmVuY2VkR0YkNiQtRiM2Ji1GVzYlUSJ5RidGWkZnbi1GLzYtUSgmbWludXM7RidGMkZORlBGUUZSRlNGVEZVL0ZFUSwwLjIyMjIyMjJlbUYnL0ZIRlxxRmNvRjJGMkYy LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYyLUkmbWZyYWNHRiQ2KC1JI21vR0YkNi1RMCZEaWZmZXJlbnRpYWxEO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R0Y3LyUqc3ltbWV0cmljR0Y3LyUobGFyZ2VvcEdGNy8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZGLUYjNiZGLi1GLzYtUSJ+RidGMi9GNlEmZmFsc2VGJy9GOUZPL0Y7Rk8vRj1GTy9GP0ZPL0ZBRk8vRkNGT0ZERkctSSNtaUdGJDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvRjNRJ2l0YWxpY0YnRjIvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRl5vLyUpYmV2ZWxsZWRHRk9GSy1GVzYlUSJ5RidGWkZnbkZLLUYvNi1RIj1GJ0YyRk5GUEZRRlJGU0ZURlUvRkVRLDAuMjc3Nzc3OGVtRicvRkhGam9GSy1GVzYlUSJyRidGWkZnbkZLLUZXNiVRInhGJ0ZaRmduLUYvNi1RKCZtaW51cztGJ0YyRk5GUEZRRlJGU0ZURlUvRkVRLDAuMjIyMjIyMmVtRicvRkhGZnBGY29GYnBGX3BGSy1GVzYlUSJ6RidGWkZnbkYy LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYwLUkmbWZyYWNHRiQ2KC1JI21vR0YkNi1RMCZEaWZmZXJlbnRpYWxEO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R0Y3LyUqc3ltbWV0cmljR0Y3LyUobGFyZ2VvcEdGNy8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZGLUYjNiZGLi1GLzYtUSJ+RidGMi9GNlEmZmFsc2VGJy9GOUZPL0Y7Rk8vRj1GTy9GP0ZPL0ZBRk8vRkNGT0ZERkctSSNtaUdGJDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvRjNRJ2l0YWxpY0YnRjIvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRl5vLyUpYmV2ZWxsZWRHRk9GSy1GVzYlUSJ6RidGWkZnbkZLLUYvNi1RIj1GJ0YyRk5GUEZRRlJGU0ZURlUvRkVRLDAuMjc3Nzc3OGVtRicvRkhGam9GSy1GVzYlUSJ4RidGWkZnbkZLLUZXNiVRInlGJ0ZaRmduLUYvNi1RKCZtaW51cztGJ0YyRk5GUEZRRlJGU0ZURlUvRkVRLDAuMjIyMjIyMmVtRicvRkhGZnAtRlc2JVEiYkYnRlpGZ25GS0Zjb0Yy

In this section we study the fixed points of the Lorenz system and their stability No numerical calculations are involved. Instead MAPLE is used to perform the necessary algebra analytically. restart: with(plots): with(linalg): Definition of the set of equations for the fixed points: lorenzfp:={sigma*(y-x)=0,r*x-y-x*z=0,x*y-b*z=0}: Solution of the systems of algebraic fixed-point equations: solcp:=solve(lorenzfp,{x,y,z}); The second solution is not unique. In such cases the command "allvalues()" gives the full set of solutions solcp2 := allvalues(solcp[2]); These solutions only exist (i.e. are real) for r \342\211\245 1. Stability analysis The Jacobian reads M:=matrix(3,3,[[-sigma,sigma,0],[r-z,-1,-x],[y,x,-b]]); The general characteristic polynomial reads CP:=det(M-lambda*diag(1,1,1)); We sort the polynomial according to the powers of lambda: CP:=collect(%,lambda); Now we assign the x,y,z-values of the first fixed point and evaluate the eigenvalues: assign(solcp[1]); x,y,z; eigenv1:=solve(CP,lambda); We find that the eigenvalue LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JlEnJiM5NTU7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYnLUkjbW5HRiQ2JFEiMUYnRjgvRjNGNy8lK2ZvcmVncm91bmRHUSxbMjAwLDAsMjAwXUYnLyUscGxhY2Vob2xkZXJHRjcvRjlRJ2l0YWxpY0YnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGOA== changes its sign at r=1: solve(eigenv1[1],r); Therefore fixed point #1 becomes instable for r>1. Now we inspect the second fixed point. The characteristic polynomial looks only slightly more complicated: x:='x':y:='y':z:='z': assign(solcp2[1]); x,y,z; CP; We simplify and sort the polynomial. It looks rather harmless: simplify(CP); charpoly2:=collect(%,lambda); An explicit solution of the characteristic equation can be found, but it is of shocking complexity: eigenv2:=solve(charpoly2,lambda); Not much can be seen from this general expression. Here the Hurwitz criterion comes to the rescue. We first extract the coefficients of the characteristic polynomial. MAPLE offers the command "coeff( )" for this task. (The minus sign is introduced to make a0 positive.) a3:=coeff(-charpoly2,lambda,3); a2:=coeff(-charpoly2,lambda,2); a1:=coeff(-charpoly2,lambda,1); a0:=coeff(-charpoly2,lambda,0); H2:=det(matrix(2,2,[[a1,a0],[a3,a2]])); The third Hurwitz determinant H3 obviously agrees with H2: H3:= det(matrix(3,3,[[a1,a0,0],[a3,a2,a1],[0,0,a3]])); The critical value of r is defined as the point where H2 (and H3) changes its sign rc:=solve(H2,r); For the 'standard parameters' this leads to the critical value sigma:=10:b:=8/3: evalf(rc); We can plot the three exact eigenvalues as a function of r. The dependence is quite complicated. However, it is clearly seen that the real part of all three eigenvalues (blue curves) is negative below for r<rc. eigen1:=unapply(eigenv2[1],r): plot([Re(eigen1(r)),Im(eigen1(r))],r=1..30,color=[blue,red]); eigen2:=unapply(eigenv2[2],r): plot([Re(eigen2(r)),Im(eigen2(r))],r=1..30,color=[blue,red]); eigen3:=unapply(eigenv2[3],r): plot([Re(eigen3(r)),Im(eigen3(r))],r=1..30,color=[blue,red]); Further study shows, that at r=13.926 an unstable limit cycle emerges which at LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEickYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUYvNiVRImNGJ0YyRjVGMi8lK2ZvcmVncm91bmRHUSxbMjAwLDAsMjAwXUYnLyUscGxhY2Vob2xkZXJHRjRGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnL0Y2USdub3JtYWxGJw===24.73... joins with the stable fixed points in a subcritical Hopf bifurcation, resulting in a strange attractor. For more details see, e.g. Strogatz, Nonlinear Dynamics and Chaos, Capt. 9.

We integrate the Lorenz equations for the standard choice of parameters s=10 and b=8/3 and the set r ={0.8, 1.2, 3., 10., 22., 28}. We prepare one-dimensional plots x(t) and for the largest values of r also three-dimensional plots (x,y,z). restart; with(plots): lorenzeqs:=diff(x(t),t)=sigma*(y(t)-x(t)),diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t),diff(z(t),t)=-b*z(t)+x(t)*y(t): vars:={x(t),y(t),z(t)}: sigma:=10: b:=8/3: rvalue:=[0.8,1.2,3.,10.,22.,28.]: ics:=x(0)=1,y(0)=1,z(0)=1: #initial conditions mycolor:=['red','green','blue','red','green','blue']: for i from 1 to 6 do r:=rvalue[i]: lorenzsol:= dsolve({lorenzeqs,ics},vars,type=numeric, method=classical[rk4],output=listprocedure): lorenzx[i]:=odeplot(lorenzsol,[t,x(t)],0..100,numpoints=500, color=mycolor[i]): if i=5 then lorenzxyz[1]:=odeplot(lorenzsol,[x(t),y(t),z(t)],0..100,numpoints=5000): fi: if i=6 then lorenzxyz[2]:=odeplot(lorenzsol,[x(t),y(t),z(t)],0..100, numpoints=5000): fi: od: We display the trajectory x(t). For r<1 the fixed point (0,0,0) is a stable attractor. For r>1 it becomes instable and is replaced by the two attractors (LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYyLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JJm1zcXJ0R0YkNiMtRiM2KS1GLDYjUSFGJy1GIzYoLUYsNiVRImJGJ0YvRjItRjY2LVExJkludmlzaWJsZVRpbWVzO0YnRjlGO0Y+RkBGQkZERkZGSC9GS1EmMC4wZW1GJy9GTkZqbi1GLDYlUSJyRidGL0YyLyUrZm9yZWdyb3VuZEdRKlswLDAsMjU1XUYnLyUwZm9udF9zdHlsZV9uYW1lR1EqMkR+T3V0cHV0RidGOS1GNjYtUSgmbWludXM7RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjIyMjIyMjJlbUYnL0ZORmlvRllGX29GYm9GOS1GNjYtUSIsRidGOUY7L0Y/RjFGQEZCRkRGRkZIRmluL0ZOUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGNUZPRltwLUYsNiVRInpGJ0YvRjJGNUZcb0Zlby1JI21uR0YkNiRRIjFGJ0Y5Rl9vRmJvRjk=) and (LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzY0LUkjbWlHRiQ2JlEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUwZm9udF9zdHlsZV9uYW1lR1EnTm9ybWFsRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LlEiPUYnRjIvRjZRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZALyUpc3RyZXRjaHlHRkAvJSpzeW1tZXRyaWNHRkAvJShsYXJnZW9wR0ZALyUubW92YWJsZWxpbWl0c0dGQC8lJ2FjY2VudEdGQC8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRk8tRjk2LlEqJnVtaW51czA7RidGMkY8Rj5GQUZDRkVGR0ZJRksvRk5RLDAuMjIyMjIyMmVtRicvRlFGVi1JJm1zcXJ0R0YkNiMtRiM2KS1GLDYjUSFGJy1GIzYoLUYsNiVRImJGJ0YvRjUtRjk2LVExJkludmlzaWJsZVRpbWVzO0YnRjxGPkZBRkNGRUZHRklGSy9GTlEmMC4wZW1GJy9GUUZjby1GLDYlUSJyRidGL0Y1LyUrZm9yZWdyb3VuZEdRKlswLDAsMjU1XUYnL0YzUSoyRH5PdXRwdXRGJ0Y8LUY5Ni1RKCZtaW51cztGJ0Y8Rj5GQUZDRkVGR0ZJRktGVUZXRlxvRmhvRltwRjwtRjk2LlEiLEYnRjJGPEY+L0ZCRjFGQ0ZFRkdGSUZLRmJvL0ZRUSwwLjMzMzMzMzNlbUYnLUYsNiZRInlGJ0YvRjJGNUY4RlJGWC1GOTYtRmJwRjxGPkZjcEZDRkVGR0ZJRktGYm9GZHAtRiw2JVEiekYnRi9GNS1GOTYtRjtGPEY+RkFGQ0ZFRkdGSUZLRk1GUEZlb0ZdcC1JI21uR0YkNiRRIjFGJ0Y8RmhvRltwRjw=). display(seq(lorenzx[i],i=1..3),view=[0..20,0..3],axes=boxed,labels=[t,x], title=`Trajectories for r=0.8,1.2,3`); With increasing r the transient time it takes to reach the stable state is growing. At r =LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbW9HRiQ2LVEifkYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R0Y0LyUqc3ltbWV0cmljR0Y0LyUobGFyZ2VvcEdGNC8lLm1vdmFibGVsaW1pdHNHRjQvJSdhY2NlbnRHRjQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZDLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEickYnLyUnaXRhbGljR1EldHJ1ZUYnL0YwUSdpdGFsaWNGJy1GIzYnLUZKNiVRImNGJ0ZNRlBGTS8lK2ZvcmVncm91bmRHUSxbMjAwLDAsMjAwXUYnLyUscGxhY2Vob2xkZXJHRk9GUC8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRi8==24.73... a bifurcation occurs and the fixed points become unstable. The trajectories approach the chaotic Lorenz attractor and display a highly irregular behaviour. display(seq(lorenzx[i],i=4..6),view=[0..100,-20..20],axes=boxed,labels=[t,x], title=`Trajectories for r=10,22,28`); For r < LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEickYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUYvNiVRImNGJ0YyRjVGMi8lK2ZvcmVncm91bmRHUSxbMjAwLDAsMjAwXUYnLyUscGxhY2Vob2xkZXJHRjRGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnL0Y2USdub3JtYWxGJw== the approach to one of the two spiral attractors is nicely seen in the 3-dimensional phase space plot: display(lorenzxyz[1],axes=boxed,labels=[x,y,z],orientation=[100,82], title=`The spiral fixed point attractor for r=20`,thickness=1,style=line); The shape of the Lorenz attractor (in the chaotic region r > LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEickYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUYvNiVRImNGJ0YyRjVGMi8lK2ZvcmVncm91bmRHUSxbMjAwLDAsMjAwXUYnLyUscGxhY2Vob2xkZXJHRjRGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnL0Y2USdub3JtYWxGJw==) is seen below. It consists of two connected leaves with the unstable fixed points in the respective centers. The trajectories unpredictably jump between the two leaves, circling aound the centers repeatedly. display(lorenzxyz[2],axes=boxed,labels=[x,y,z],orientation=[100,82], title=`The strange attractor for r=28`,thickness=1,style=line);

We repeat the calculation with a manual imput of r.Regions of chaotic and regular trajectories are intermixed. Suggested values, e.g., r = 100, 120, 150, 170, 200, 240 restart; with(plots): lorenzeqs:=diff(x(t),t)=sigma*(y(t)-x(t)),diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t),diff(z(t),t)=-b*z(t)+x(t)*y(t): vars:={x(t),y(t),z(t)}: sigma:=10:b:=8/3: ics:=x(0)=10,y(0)=10,z(0)=10: #initial conditions tstart:=30: #don't display the initial period r:=100; lorenzsol:= dsolve({lorenzeqs,ics},vars,type=numeric,relerr=1e-12,maxfun=200000,output=listprocedure): lorenzx:=odeplot(lorenzsol,[t,x(t)],tstart..50,numpoints=1000,color=blue): lorenzxyz:=odeplot(lorenzsol,[x(t),y(t),z(t)],tstart..100,numpoints=10000): Note that an initial transient time period ( t<tstart) is not displayed. display(lorenzx,axes=boxed,labels=[t,x]); display(lorenzxyz,axes=boxed,labels=[x,y,z],thickness=1,style=line); Period doubling bifurcations With varying parameter r there is a complicated sequence of regular (periodic) and chaotic (aperiodic) attractors. For the periodic trajectories, sequences of period doublings can be identified. Actually they are found for decreasing values of r. For examples we can observe the cascade: r = 99.980 period 2 r = 99.630 period 4 r = 99.548 period 8 Transient chaos Even when a trajectory finally reaches a stable fixed point it can show transient chaotic behaviour before it finally settles down to the stable state. As an example look at r = 22.4 tstart = 0 tmax = 200

Instead of going through a bifurcation, a system also can become chaotic through the phenomenon of intermittency. The 'intermittency route to chaos' is observed in the Lorenz system in the vicinity of r=166. Starting from a stable periodic motion an increase of r leads to sudden 'outbursts'. The time interval between these outbursts shrinks until the trajectory becomes completely irregular. restart; with(plots): lorenzeqs:=diff(x(t),t)=sigma*(y(t)-x(t)),diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t),diff(z(t),t)=-b*z(t)+x(t)*y(t): vars:={x(t),y(t),z(t)}: sigma:=10:b:=8/3: rvalue:=[166.,166.2,166.4]: ics:=x(0)=10,y(0)=10,z(0)=10: #initial conditions for i from 1 to 3 do r:=rvalue[i]: lorenzsol:= dsolve({lorenzeqs,ics},vars,type=numeric,method=classical[rk4],output=listprocedure): lorenzz[i]:=odeplot(lorenzsol,[t,z(t)],0..100,numpoints=10000,color=blue): od: for i from 1 to 3 do display(lorenzz[i],view=[30..100,60..250], axes=boxed,labels=[t,z], title=`Trajectories with intermittency`); od;

The Lorenz map (or return map) reduces the full dynamics of the Lorenz system to a done-dimensional discrete map z(n) -> z(n+1) where z(n) is the nth maximum of the trajectory z(t). It turns out that this map to a good approximation can be described by a simple continuous (but not everywhere differentiable) function: z(n+1) = f(z(n))! Thus part of the dynamics of the Lorenz system can be understood in terms of a discrete one-dimensional iteration, without the need to solve the differential equations! restart; with(plots): lorenzeqs:=diff(x(t),t)=sigma*(y(t)-x(t)),diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t),diff(z(t),t)=-b*z(t)+x(t)*y(t): vars:={x(t),y(t),z(t)}: ics:=x(0)=10,y(0)=10,z(0)=10: sigma:=10: b:=8/3: r:=28: tmax:=300: dt:=0.01: numb:=0: lorenzsol:=dsolve({lorenzeqs,ics},vars,numeric,method=rkf45,relerr=1e-12,output=listprocedure): assign(lorenzsol); Z:=z(t): imax := tmax/dt: tstart:=10: for i from 1 to imax do # generate a list of z(t_i) values zarray[i]:= Z(tstart+i*dt); od: for i from 1 to imax-2 do # search for the maxima and put them into the list pts if (zarray[i+1]-zarray[i]>0 and zarray[i+2]-zarray[i+1]<0) then numb:=numb+1: pts[numb]:=zarray[i+1]: fi; od: The points of the return map fall on a smooth curve (with very small deviations) with a sharp peak. plot([seq([pts[j],pts[j+1]],j=0..(numb-1))],style=point,symbol=POINT,color=blue,view=[30..50,30..50],title=`Lorenz map for r=28`,axes=boxed,labels=[`z(n)`,`z(n+1)`]); JSFH

To find the largest Lyapunov exponent we solve a 6-dimensional system of differential equations: We simultaneously integrate the reference trajectory LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUYvNiVRInJGJ0YyRjVGMi8lK2ZvcmVncm91bmRHUSxbMjAwLDAsMjAwXUYnLyUscGxhY2Vob2xkZXJHRjRGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnL0Y2USdub3JtYWxGJw==(t) solving the nonlinear equation of motion and the vector of perturbations \316\267(t) solving the linearized system. Separately calculating the reference trajectory (3-dimensional nonlinear diff. eq.), storing it, and then solving for the perturbations (3-dimensional linear diff. eq.) using this trajectory in principle should be more efficient but is less convenient to implement. restart; with(plots): lorenzeqs:= diff(xr(t),t)=sigma*(yr(t)-xr(t)), diff(yr(t),t)=r*xr(t)-yr(t)-xr(t)*zr(t), diff(zr(t),t)=-b*zr(t)+xr(t)*yr(t): lineqs:= diff(eta1(t),t)=-sigma*eta1(t)+sigma*eta2(t), diff(eta2(t),t)=(r-zr(t))*eta1(t)-eta2(t)-xr(t)*eta3(t), diff(eta3(t),t)=yr(t)*eta1(t)+xr(t)*eta2(t)-b*eta3(t): vars:= {xr(t),yr(t),zr(t),eta1(t),eta2(t),eta3(t)}: ics:= xr(0)=10,yr(0)=10,zr(0)=10, # initial conditions (nonlinear) eta1(0)=1e-10,eta2(0)=1e-10,eta3(0)=1e-10: # initial conditions (linear) sigma:=10: b:=8/3: r:=28: sol:=dsolve({lorenzeqs,lineqs,ics},vars,numeric,relerr=1e-12,maxfun=500000,output=listprocedure): e1:=subs(sol,eta1(t)): e2:=subs(sol,eta2(t)): e3:=subs(sol,eta3(t)): normeta:= t -> sqrt(e1(t)*e1(t)+e2(t)*e2(t)+e3(t)*e3(t)): ploteta:=logplot(normeta(t),t=0..50): display(ploteta,axes=boxed,labels=[t,`|eta|`], title= `Linear solution`); Approximate evaluation of the largest Lyapunov exponent LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JlEnJiM5NjM7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYoLUkjbW5HRiQ2JFEiMUYnRjgtSSNtb0dGJDYtUSJ+RidGOC8lJmZlbmNlR0Y0LyUqc2VwYXJhdG9yR0Y0LyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlUvRjNGNy8lK2ZvcmVncm91bmRHUSxbMjAwLDAsMjAwXUYnLyUscGxhY2Vob2xkZXJHRjcvRjlRJ2l0YWxpY0YnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGOA=== LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUknbXVuZGVyR0YkNiUtSSNtb0dGJDYtUSRsaW1GJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJnVuc2V0RicvJSpzZXBhcmF0b3JHRjcvJSlzdHJldGNoeUdGNy8lKnN5bW1ldHJpY0dGNy8lKGxhcmdlb3BHRjcvJS5tb3ZhYmxlbGltaXRzR1EldHJ1ZUYnLyUnYWNjZW50R0Y3LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMTY2NjY2N2VtRictRiM2Ji1JI21pR0YkNiVRInRGJy8lJ2l0YWxpY0dGQi9GM1EnaXRhbGljRictRi82LVEtJnJpZ2h0YXJyb3c7RidGMi9GNlEmZmFsc2VGJy9GOUZZL0Y7RkIvRj1GWS9GP0ZZL0ZBRlkvRkRGWS9GRlEsMC4yNzc3Nzc4ZW1GJy9GSUZbby1GLzYtUSgmaW5maW47RidGMkZYRlovRjtGWUZmbkZnbkZobkZpbkZFL0ZJRkdGMi8lLGFjY2VudHVuZGVyR0ZZLUZONiNRIUYnLUYjNiVGZG8tSSZtZnJhY0dGJDYoLUYjNigtSSNtbkdGJDYkUSIxRidGMkZRLyUrZm9yZWdyb3VuZEdRK1swLDE2MCw4MF1GJy8lLHBsYWNlaG9sZGVyR0ZCLyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR0ZCRlMtRiM2J0ZNRlEvRmNwUSxbMjAwLDAsMjAwXUYnRmVwRlMvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmJxLyUpYmV2ZWxsZWRHRllGMkZkb0Yy LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictSSNtb0dGJDYtUSNsbkYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGOC8lKXN0cmV0Y2h5R0Y4LyUqc3ltbWV0cmljR0Y4LyUobGFyZ2VvcEdGOC8lLm1vdmFibGVsaW1pdHNHRjgvJSdhY2NlbnRHRjgvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZHLUYwNi1RMCZBcHBseUZ1bmN0aW9uO0YnRjMvRjdRJmZhbHNlRicvRjpGTi9GPEZOL0Y+Rk4vRkBGTi9GQkZOL0ZERk5GRUZILUkobWZlbmNlZEdGJDYkLUYjNiZGKy1GIzYlRistSSZtZnJhY0dGJDYoLUYjNiUtRlY2Ji1GIzYmLUYsNiVRJyYjOTUxO0YnLyUnaXRhbGljR0ZORjMtRlY2JC1GIzYlLUYsNiVRInRGJy9GY29RJXRydWVGJy9GNFEnaXRhbGljRidGW3BGXXBGM0ZbcEZdcEYzLyUlb3BlbkdRInxnckYnLyUmY2xvc2VHRmFwRltwRl1wLUYjNictRlY2Ji1GIzYoRl9vLUZWNiQtRiM2Jy1JI21uR0YkNiRRIjBGJ0YzRltwLyUrZm9yZWdyb3VuZEdRLFsyMDAsMCwyMDBdRicvJSxwbGFjZWhvbGRlckdGXHBGXXBGM0ZbcEZicUZlcUZdcEYzRl9wRmJwRltwRmJxRmVxRl1wLyUubGluZXRoaWNrbmVzc0dRIjFGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zcci8lKWJldmVsbGVkR0ZORjNGK0YzRjNGMw== lyapunov:= t-> log(normeta(t)/normeta(0))/t: lyapunov(50); lyapunov(100); lyapunov(200); Remark: In parameter regions leading to non-chaotic trajectories, both 'nontrivial' Lyapunov exponents are negative. Our program then finds the value LUklbXN1Ykc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JlEoJnNpZ21hO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lNnNlbGVjdGlvbi1wbGFjZWhvbGRlckdRJXRydWVGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSVtcm93R0YkNiYtSSNtbkdGJDYkUSIxRidGNS1JI21vR0YkNi1RIn5GJ0Y1LyUmZmVuY2VHRjEvJSpzZXBhcmF0b3JHRjEvJSlzdHJldGNoeUdGMS8lKnN5bW1ldHJpY0dGMS8lKGxhcmdlb3BHRjEvJS5tb3ZhYmxlbGltaXRzR0YxLyUnYWNjZW50R0YxLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGUy8lMGZvbnRfc3R5bGVfbmFtZUdRJVRleHRGJ0Y1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRic= = 0 which in this situation is the largest Lyapunov exponent. Try out the value r = 100. LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=