Orbits of the Lorenz Equation

We will compute some orbits of the famous Lorenz equation.  This defines a flow in R^3 with a strange attractor.  In other words, there is a very complicated invariant set which attracts and forms the omega limit set for most initial conditions.

lorenzeq = {x'[t]==10*(y[t]-x[t]),y'[t]==x[t]*(28-z[t])-y[t],z'[t]==x[t]*y[t]-8*z[t]/3}

Here is a function which computes an orbit and draws it using parametric plot.

MakeLorenzPicture[t1_,x0_,y0_,z0_]:= Module[{sol,ppl,x1,y1,z1},
sol=NDSolveValue[lorenzeq~Join~{x[0]==x0,y[0]==y0,z[0]==z0},{x[t],y[t],z[t]},{t,0,Max[t1,0.01]}];
ppl=ParametricPlot3D[sol,{t,0,Max[t1,0.01]},PlotStyle->{Thin,Blue},PlotRange->{{-30,30},{-30,30},{0,50}}];
{x1,y1,z1}=sol/.t->Max[t1,0.01];
Show[ppl,Graphics3D[{PointSize[0.02],Green,Point[{x0,y0,z0}],Black,Point[{x1,y1,z1}]}]]
]

The forward orbit of a typical initial points shows a very complicated structure.

MakeLorenzPicture[30,0,10,50]

To see the omega limit set we can just plot the tail of the orbit.  In the following function we only plot the part of the orbit from t1 - tail to t1.

MakeLorenzPictureTail[t1_,x0_,y0_,z0_,tail_]:= Module[{sol,ppl,x1,y1,z1,t1alt,t0},
t1alt = Max[t1,0.01];
t0=Max[0.001,t1alt-tail];
sol=NDSolveValue[lorenzeq~Join~{x[0]==x0,y[0]==y0,z[0]==z0},{x[t],y[t],z[t]},{t,0,t1alt},MaxStepSize->0.2];
ppl=ParametricPlot3D[sol,{t,t0,t1alt},PlotStyle->{Thin,Blue},PlotRange->{{-30,30},{-30,30},{0,50}},MaxRecursion->10];
{x1,y1,z1}=sol/.t->t1alt;
Show[ppl,Graphics3D[{PointSize[0.02],Black,Point[{x1,y1,z1}]}]]
]

Here we skip the first 50 time units to give the orbit time to get attracted and then plot the next 50.

attractorpic =MakeLorenzPictureTail[250,0,10,50,50]

By plotting with a short tail you can visualize the motion.

Manipulate[MakeLorenzPictureTail[t1,x,y,z,tail],{t1,0.01,100},{{x,0.1},-10,10},{y,-10,10},{{z,10},0,40},{{tail,5},0.01,10}]

Poincare Map

We will use the plane z=27 as a Poincare section.  Every time an orbit crosses the plane from up to down, we get a new point of the section.

sigma = Graphics3D[Polygon[{{-7,7,27},{7,7,27},{7,-7,27},{-7,-7,27}}]];
Show[attractorpic,sigma]

The next function computes the next intersection of a solution with the Poincare section.

NextPoincarePoint[{x0_,y0_,z0_}]:=
(Reap@NDSolve[lorenzeq~Join~{x[0]==x0,y[0]==y0,z[0]==z0,WhenEvent[z[t]<27&&t>0.1,Sow[{x[t],y[t],z[t]}];"StopIntegration"]},{x[t],y[t],z[t]},{t,0,10}])[[2,1,1]]

Make a list of a lot of Poincare points.

start = NextPoincarePoint[{0,10,50}];
poincarelist= NestList[NextPoincarePoint,start,1000];
poincareplot = ListPlot[Map[Part[#,{1,2}]&,poincarelist],AspectRatio->1,AxesLabel->{"x","y"}]

Thus the attractor intersects the section in what looks like two curves.  Actually the two sets are not curves but have many layers.  However, thinking of curves gives a good approximation. Now roughly speaking, the Poincare map smashes these two curves together to make one curve, then tears the new curve in half and stretches the halves onto the original two curves again.  All of this can be described approximately by a one-dimensional discrete dynamical system.  One variable which can describe the map is u = .36 x+y (found by trial and error).  The level lines of this are shown below.  Points on the same contour line will end up getting approximately smashed together after one iteration of the Poincare map.  The level curve u = 0 is approximately the stable manifold of (0,0,0) and this is where the map tears things in half.  The points with u>0 go around one side of the attractor and then return while points with u<0 go around the other side and then return.

ContourPlot[.36*x+y,{x,-6,6},{y,-6,6},ContourShading->False,Contours->20];
Show[poincareplot,%]

To understand the dynamics we can make a delay plot showing u(n) versus u(n+1).   Iteration of a map like this defines a “chaotic” discrete dynamical system.  The u coordinate is repeatedly stretched and cut.  After a while, it is hard to predict where you will be, or even whether u>0 or u<0.

u[{x_,y_,z_}]:= 0.36*x+y;
ulist= Map[u,poincarelist];
delaylist={Drop[ulist,-1],Drop[ulist,1]}//Transpose;
ListPlot[delaylist,AspectRatio->1,AxesLabel->{"u(n)","u(n+1)"}]