Chaotic systems

The Lorenz attractor

Clear[tmax]
Tmax=100;
x0=1.1;y0=1;z0=1;

d=3;
r=26.5;
b=1;

Manipulate[
soln=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,tmax},MaxSteps->Infinity];
px=Plot[Evaluate[x[t]/.soln],{t,0,tmax},AxesLabel-> {t,x[t]},PlotRange-> All,ImageSize->280,LabelStyle->20];
py=Plot[Evaluate[y[t]/.soln],{t,0,tmax},AxesLabel-> {t,y[t]},PlotRange-> All,ImageSize->280,LabelStyle->20];
pz=Plot[Evaluate[z[t]/.soln],{t,0,tmax},AxesLabel-> {t,z[t]},PlotRange-> All,ImageSize->280,LabelStyle->20];
Row[{px,py,pz}],
{{tmax,0.25 Tmax,"Integration time, t"},0,Tmax,Appearance->"Labeled"},SaveDefinitions->True,TrackedSymbols:>{tmax}]

soln2=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,Tmax},MaxSteps->Infinity];
diffx=d (y-x);diffy=r x-y-x z;diffz=x y-b z;
fixedpoints=Quiet[Solve[{diffx==0,diffy==0,diffz==0},{x,y,z}]];
Print[Style["Coordinates of fixed points",18,Blue]]
coordinates={x,y,z}/.fixedpoints
Gfp=Graphics3D[{PointSize[0.03],Red,Point[coordinates]}];
pplot=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.soln2],{t,0,Tmax},PlotPoints->1000,AxesLabel->{x[t],y[t],z[t]},PlotStyle->Thin,LabelStyle->20];
Show[pplot,Gfp]

The Lorenz attractor, two initial conditions

Clear[tmax]
Tmax=40;
x0=1.1;y0=1;z0=1;
dx0=0.0001;dy0=0.0;dz0=0.0;

d=3;
r=26.5;
b=1;

Manipulate[
sol1=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,tmax},MaxSteps->Infinity];
sol2=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0+dx0,y[0]==y0+dy0,z[0]==z0+dz0},{x,y,z},{t,0,tmax},MaxSteps->Infinity];
px1=Plot[{Evaluate[x[t]/.sol1],Evaluate[x[t]/.sol2]},{t,0,tmax},AxesLabel-> {t,x[t]},PlotRange-> All,ImageSize->270,LabelStyle->20];
py2=Plot[{Evaluate[y[t]/.sol1],Evaluate[y[t]/.sol2]},{t,0,tmax},AxesLabel-> {t,y[t]},PlotRange-> All,ImageSize->270,LabelStyle->20];
pz3=Plot[{Evaluate[z[t]/.sol1],Evaluate[z[t]/.sol2]},{t,0,tmax},AxesLabel-> {t,z[t]},PlotRange-> All,ImageSize->270,LabelStyle->20];
Row[{px1,py2,pz3}],
{{tmax, 0.25 Tmax,"Integration time, t"},0,Tmax,Appearance->"Labeled"},TrackedSymbols:>{tmax},SaveDefinitions->True]

sol12=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,Tmax},MaxSteps->Infinity];
sol22=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0+dx0,y[0]==y0+dy0,z[0]==z0+dz0},{x,y,z},{t,0,Tmax},MaxSteps->Infinity];
diffx=d (y-x);diffy=r x-y-x z;diffz=x y-b z;
fixedpoints=Quiet[Solve[{diffx==0,diffy==0,diffz==0},{x,y,z}]];
Print[Style["Coordinates of fixed points",18,Blue]]
coordinates={x,y,z}/.fixedpoints
Gfp=Graphics3D[{PointSize[0.03],Red,Point[coordinates]}];
pplot1=ParametricPlot3D[{Evaluate[{x[t],y[t],z[t]}/.sol12],Evaluate[{x[t],y[t],z[t]}/.sol22]},
{t,0,Tmax},PlotPoints->2000,AxesLabel->{x[t],y[t],z[t]},PlotStyle->Thin,LabelStyle->20];
Show[pplot1,Gfp]

Animation of the Lorenz attractor in 2-D projection

Clear[tmax]
Tmax=25;
x0=1.1;y0=1;z0=7;
xlimit=12;zlimit=42;

d=3;
r=26.5;
b=1;

Animate[
intsol1=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,tmax},MaxSteps->Infinity];
intsol2=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0+dx0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,tmax},MaxSteps->Infinity];
pxz1=ParametricPlot[Evaluate[{x[t],z[t]}/.intsol1],{t,0,tmax},AxesLabel-> {x[t],z[t]},PlotRange->{{-xlimit,xlimit},{5,zlimit}},AspectRatio->1,PlotPoints->450,PlotStyle->{Thickness[0.0025]},ImageSize->530,LabelStyle->20,PerformanceGoal->"Quality"];
pxz2=ParametricPlot[Evaluate[{x[t],z[t]}/.intsol2],{t,0,tmax},AxesLabel-> {x[t],z[t]},PlotRange-> {{-xlimit,xlimit},{5,zlimit}},AspectRatio->1,PlotPoints->450,PlotStyle->{Thickness[0.0015],Red},ImageSize->530,LabelStyle->20,PerformanceGoal->"Quality"];
Show[pxz1,pxz2,Graphics[{PointSize[0.027],Blue,Point[Dynamic[Evaluate[{x[tmax],z[tmax]}/.intsol1]]]}],Graphics[{PointSize[0.027],Red,Point[Dynamic[Evaluate[{x[tmax],z[tmax]}/.intsol2]]]}]],
{tmax,0.01,Tmax},AnimationRunning->False,AnimationRate->0.6,SaveDefinitions->True]

Animation of the Lorenz attractor in 3-D, two initial conditions

Clear[dt]
Tmax=25;
x0=1.1;y0=1;z0=7;
xlimit=12;ylimit=20;zlimit=42;

d=3;
r=26.5;
b=1;

Animate[
intsol1s=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,tmax},MaxSteps->Infinity];
intsol2s=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0+dx0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,tmax},MaxSteps->Infinity];
pxyz1=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.intsol1s],{t,0,tmax},AxesLabel-> {x[t],y[t],z[t]},PlotRange->{{-xlimit,xlimit},{-ylimit,ylimit},{5,zlimit}},BoxRatios->{1,1,1},PlotPoints->450,PlotStyle->{Thickness[0.0029]},ImageSize->530,LabelStyle->20,PerformanceGoal->"Quality"];
pxyz2=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.intsol2s],{t,0,tmax},AxesLabel-> {x[t],y[t],z[t]},PlotRange-> {{-xlimit,xlimit},{-ylimit,ylimit},{5,zlimit}},BoxRatios->{1,1,1},PlotPoints->450,PlotStyle->{Thickness[0.0029],Red},ImageSize->530,LabelStyle->20,PerformanceGoal->"Quality"];
Show[pxyz1,pxyz2,Graphics3D[{PointSize[0.027],Blue,Point[Dynamic[Evaluate[{x[tmax],y[tmax],z[tmax]}/.intsol1s]]]}],Graphics3D[{PointSize[0.027],Red,Point[Dynamic[Evaluate[{x[tmax],y[tmax],z[tmax]}/.intsol2s]]]}]],
{tmax,0.01,Tmax},AnimationRunning->False,AnimationRate->0.6,SaveDefinitions->True]

Creation of this notebook was supported by IT Academy program of Information Technology Foundation for Education (HITSA).