The Lorenz Attractor

- Marcus McGuff, Austin Community College

Manipulate[Refresh[If[var===Null,Clear[var];var["computePlotrange"]=plotrange;
var["highrect"]={}],None];


Refresh[(*Print["Eq pt solve"];*)

eqLorenz={x'[t]==s (-x[t]+y[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==-b z[t]+x[t] y[t]};

eqPts=Solve[eqLorenz/.{x'[t]->0,y'[t]->0,z'[t]->0},{x[t],y[t],z[t]}];

(*Make a function to find the eq pt nearest to any other pt (exclude origin)*)

If[Length[eqPts]>1,
nearestEqPt=Nearest[N[Select[({x[t],y[t],z[t]}/.eqPts),#≠{0.,0.,0.}&]]],nearestEqPt=Nearest[{{0,0,0}}]];

(*localLin[[j]] is lin matirx for pt k*)localLin=D[eqLorenz[[All,2]],{{x[t],y[t],z[t]}}]/.eqPts//N;

(*eigenSets[[k]] is eigset for pt k,eigenSets[[k,1]] is set of eigenvalues-eigenSets[[k,1,j]] is jth eigenvalue for pt k,eigenSets[[k,2]] is set of all eigenvectors,eigenSets[[k,2,j]] is jth eigenvector*)eigenSets=Map[Eigensystem,localLin],
TrackedSymbols:>{rPar,sPar,bPar}];

(* 
solveLorenzIC[{x0,y0,z0},currentPos,0];

solveLorenzIC[{xIC_,yIC_,zIC_},end_,pert_]:=First[{x[t],y[t],z[t]}/.Quiet[NDSolve[SetPrecision[{Dynamic[eqLorenz],Thread[{x[0],y[0],z[0]}==N[({xIC,yIC,zIC}+pert*RandomInteger[1,3])]]},Infinity],{x,y,z},{t,0,end},PrecisionGoal->ControlActive[4,8],WorkingPrecision->ControlActive[MachinePrecision,20]]]]

*)

Refresh[(*Print["DE solve"];*)
solLorenz=First[{x[t],y[t],z[t]}/.Quiet[NDSolve[SetPrecision[{eqLorenz,Thread[{x[0],y[0],z[0]}=={x0,y0,z0}]},Infinity],{x,y,z},{t,0,currentPos},PrecisionGoal->ControlActive[4,8],WorkingPrecision->ControlActive[MachinePrecision,20]]]];

endpt=Min[Flatten[Thread@DifferentialEquations`InterpolatingFunctionAnatomy`InterpolatingFunctionDomain[solLorenz[[All,0]]],1][[All,2]]],
TrackedSymbols:>{x0,y0,z0,eqLorenz,currentPos}];

With[{
optGr2D=Sequence[Axes->True,AxesOrigin->{0,0},AspectRatio->0.5,ImageSize->150,Ticks->None],plotOpt=Sequence[ColorFunction->(Hue[#1]&),PlotPoints->ControlActive[30,50],ColorFunctionScaling->False],
parPlotOpt=Sequence[PlotPoints->ControlActive[100,100],MaxRecursion->ControlActive[4,4],ColorFunctionScaling->False,ColorFunction->(Hue[#3]&)]},

Switch[ctrlPane,
2,
Dynamic@Grid[(*Print["Main loop 2"];*){{Graphics3D[{Dynamic[First[ParametricPlot3D[solLorenz,{t,0,endpt},PlotPoints->ControlActive[100,150],MaxRecursion->ControlActive[4,5],ColorFunction->(Hue[#4]&),ColorFunctionScaling->False]]],

{Red,PointSize[Large],Point[solLorenz/.t->Evaluate@Min[currentPos,endpt]]},

Dynamic[Table[{Blue,PointSize[Large],Tooltip[Point[{x[t],y[t],z[t]}/.eqPts[[k]]],Text@Style[Row[{"Eigenvalues: ",eigenSets[[k,1]]}],FontSize->12]]},{k,1,Length[eqPts]}]],

With[{curpt=(solLorenz/.t->Evaluate@Min[currentPos,endpt])},Tooltip[{Blue,Dashed,Arrow[{curpt,First[nearestEqPt[curpt]]}]},First[nearestEqPt[curpt]]]]},
ImageSize->400,AxesLabel->{"X","Y","Z"},Axes->True],

Graphics[{First[Plot[solLorenz[[1]],{t,0,endpt},plotOpt]],
{Red,PointSize[Large],Tooltip[Point[{t,solLorenz[[1]]}/.t->Evaluate@Min[currentPos,endpt]],{t,solLorenz[[1]]}/.t->Evaluate@Min[currentPos,endpt]]}},
AxesLabel->{"t","X"},Axes->True,optGr2D]},

{SpanFromAbove,
Graphics[{First[Plot[solLorenz[[2]],{t,0,endpt},plotOpt]],
{Red,PointSize[Large],Tooltip[Point[{t,solLorenz[[2]]}/.t->Evaluate@Min[currentPos,endpt]],{t,solLorenz[[2]]}/.t->Evaluate@Min[currentPos,endpt]]}},
AxesLabel->{"t","Y"},Axes->True,optGr2D]},

{SpanFromAbove,
Graphics[{First[Plot[solLorenz[[3]],{t,0,endpt},plotOpt]],
{Red,PointSize[Large],Tooltip[Point[{t,solLorenz[[3]]}/.t->Evaluate@Min[currentPos,endpt]],{t,solLorenz[[3]]}/.t->Evaluate@Min[currentPos,endpt]]}},
AxesLabel->{"t","Z"},Axes->True,optGr2D]}}],

1,
Dynamic@Grid[(*Print["Main loop 1"];*)Transpose@{{
Graphics[{First[ParametricPlot[Evaluate[solLorenz[[1;;2]]],{t,0,endpt},parPlotOpt]],

{Red,PointSize[Large],Tooltip[Point[solLorenz[[1;;2]]/.t->Min[currentPos,endpt]],solLorenz[[1;;2]]/.t->Min[currentPos,endpt]]},

Dynamic[
Table[{Blue,PointSize[Large],
Tooltip[Point[{x[t],y[t]}/.eqPts[[k]]],Text@Style[Row[{"Eigenvalues: ",eigenSets[[k,1]]}],FontSize->12]]},{k,1,3}]]


},
ImageSize->150,AxesLabel->{"X","Y"},Axes->True
],

Graphics[{First[Plot[solLorenz[[1]],{t,0,endpt},plotOpt]],
{Red,PointSize[Large],Tooltip[Point[{t,solLorenz[[1]]}/.t->Evaluate@Min[currentPos,endpt]],{t,solLorenz[[1]]}/.t->Evaluate@Min[currentPos,endpt]]}},
AxesLabel->{"t","X"},Axes->True,optGr2D]},

{Graphics[{First[ParametricPlot[Evaluate[solLorenz[[1;;3;;2]]],{t,0,endpt},parPlotOpt]],

{Red,PointSize[Large],Tooltip[Point[solLorenz[[1;;3;;2]]/.t->Min[currentPos,endpt]],solLorenz[[1;;3;;2]]/.t->Min[currentPos,endpt]]},

Dynamic[
Table[{Blue,PointSize[Large],
Tooltip[Point[{x[t],z[t]}/.eqPts[[k]]],Text@Style[Row[{"Eigenvalues: ",eigenSets[[k,1]]}],FontSize->12]]},{k,1,3}]]


},
ImageSize->150,AxesLabel->{"X","Z"},Axes->True
],

Graphics[{First[Plot[solLorenz[[2]],{t,0,endpt},plotOpt]],
{Red,PointSize[Large],Tooltip[Point[{t,solLorenz[[2]]}/.t->Evaluate@Min[currentPos,endpt]],{t,solLorenz[[2]]}/.t->Evaluate@Min[currentPos,endpt]]}},
AxesLabel->{"t","Y"},Axes->True,optGr2D]},

{Graphics[{First[ParametricPlot[Evaluate[solLorenz[[2;;3]]],{t,0,endpt},parPlotOpt]],

{Red,PointSize[Large],Tooltip[Point[solLorenz[[2;;3]]/.t->Min[currentPos,endpt]],solLorenz[[2;;3]]/.t->Min[currentPos,endpt]]},

Dynamic[
Table[{Blue,PointSize[Large],
Tooltip[Point[{y[t],z[t]}/.eqPts[[k]]],Text@Style[Row[{"Eigenvalues: ",eigenSets[[k,1]]}],FontSize->12]]},{k,1,3}]]


},
ImageSize->150,AxesLabel->{"Y","Z"},Axes->True
],

Graphics[{First[Plot[solLorenz[[3]],{t,0,endpt},plotOpt]],
{Red,PointSize[Large],Tooltip[Point[{t,solLorenz[[3]]}/.t->Evaluate@Min[currentPos,endpt]],{t,solLorenz[[3]]}/.t->Evaluate@Min[currentPos,endpt]]}},
AxesLabel->{"t","Z"},Axes->True,optGr2D]}
},Dividers->{False,{False,True,False}}],

3,Dynamic[(*Print["Main loop 3"];*)
With[{solPert1=First[{x[t],y[t],z[t]}/.Quiet[NDSolve[SetPrecision[{eqLorenz,Thread[{x[0],y[0],z[0]}=={x0+perturbation RandomInteger[1],y0+perturbation RandomInteger[1],z0+perturbation RandomInteger[1]}]},Infinity],{x,y,z},{t,0,currentPos},PrecisionGoal->ControlActive[4,8],WorkingPrecision->ControlActive[MachinePrecision,20]]]],solPert2=First[{x[t],y[t],z[t]}/.Quiet[NDSolve[SetPrecision[{eqLorenz,Thread[{x[0],y[0],z[0]}=={x0+2perturbation RandomInteger[1],y0+2perturbation RandomInteger[1],z0+2perturbation RandomInteger[1]}]},Infinity],{x,y,z},{t,0,currentPos},PrecisionGoal->ControlActive[4,8],WorkingPrecision->ControlActive[MachinePrecision,20]]]],
plotOpt2=Sequence[PlotPoints->ControlActive[30,30]],parPlotOpt2=Sequence[PlotPoints->ControlActive[80,80],MaxRecursion->ControlActive[4,4]]},
Grid[{{Graphics[{First@DensityPlot[Evaluate[-Sign[solLorenz[[2]]/.t->currentPos]x+2(1-(y+If[(solLorenz[[3]]/.t->currentPos)>200,1/3.,(solLorenz[[3]]/.t->currentPos)/600] Sin[π y]))],{x,-1.7,1.7},{y,0,1.1},ColorFunction->ColorData["TemperatureMap"],Frame->None,AspectRatio->1/3],

({Red,PointSize[Large],
Arrowheads[
     {{-Sign[solPert1[[1]]]0.04, 1}, {-Sign[solPert1[[1]]]0.04,0.5},
      {-Sign[solPert1[[1]]]0.04, 0}}],Arrow[{0.5Sin[solPert1[[1]]t-tLag]+xShift,0.5Cos[solPert1[[1]]t-tLag]+0.55}/.tLag->Range[0,4,0.2]/.xShift->-1.05//Transpose],

Green,
Arrowheads[
     {{-Sign[solLorenz[[1]]]0.04, 1}, {-Sign[solLorenz[[1]]]0.04,0.5},
      {-Sign[solLorenz[[1]]]0.04, 0}}],
Arrow[{0.5Sin[solLorenz[[1]]t-tLag]+xShift,0.5Cos[solLorenz[[1]]t-tLag]+0.55}/.tLag->Range[0,4,0.2]/.xShift->0//Transpose],

Blue,
Arrowheads[
     {{-Sign[solPert2[[1]]]0.04, 1}, {-Sign[solPert2[[1]]]0.04,0.5},
      {-Sign[solPert2[[1]]]0.04, 0}}],Arrow[{0.5Sin[solPert2[[1]]t-tLag]+xShift,0.5Cos[solPert2[[1]]t-tLag]+0.55}/.tLag->Range[0,4,0.2]/.xShift->1.05//Transpose],
Thick,Line[{{-1.7,0},{1.7,0}}],
Line[{{-1.7,1.1},{1.7,1.1}}]}/.t->currentPos)
},ImageSize->500],
SpanFromLeft,SpanFromLeft},

{Text@Row[{Style["Perturbed",Red],"        ",Style["Initial",Green],"        ",Style["Extremely Perturbed",Blue]}],
SpanFromLeft,SpanFromLeft},

Switch[whichView,
"Time plots",
{Graphics[Dynamic[First[Plot[solLorenz[[1]],{t,0,endpt},plotOpt2]]],AxesLabel->{"t","X"},Axes->True,optGr2D],

Graphics[Dynamic[First[Plot[solLorenz[[2]],{t,0,endpt},plotOpt2]]],AxesLabel->{"t","Y"},Axes->True,optGr2D],

Graphics[Dynamic[First[Plot[solLorenz[[3]],{t,0,endpt},plotOpt2]]],AxesLabel->{"t","Z"},Axes->True,optGr2D]
},

"Phase plots",
{Graphics[{Dynamic[First[ParametricPlot[Evaluate[solLorenz[[1;;2]]],{t,0,endpt},parPlotOpt2]]],

{Red,PointSize[Large],Point[solLorenz[[1;;2]]/.t->Min[currentPos,endpt]]},

Dynamic[
Table[{Blue,PointSize[Large],
Tooltip[Point[{x[t],y[t]}/.eqPts[[k]]],Text@Style[Row[{"Eigenvalues: ",eigenSets[[k,1]]}],FontSize->12]]},{k,1,3}]]
},
ImageSize->150,AxesLabel->{"X","Y"},Axes->True
],

Graphics[{Dynamic[First[ParametricPlot[Evaluate[solLorenz[[1;;3;;2]]],{t,0,endpt},parPlotOpt]]],

{Red,PointSize[Large],Point[solLorenz[[1;;3;;2]]/.t->Min[currentPos,endpt]]},

Dynamic[
Table[{Blue,PointSize[Large],
Tooltip[Point[{x[t],z[t]}/.eqPts[[k]]],Text@Style[Row[{"Eigenvalues: ",eigenSets[[k,1]]}],FontSize->12]]},{k,1,3}]]
},
ImageSize->150,AxesLabel->{"X","Z"},Axes->True
],

Graphics[{Dynamic[First[ParametricPlot[Evaluate[solLorenz[[2;;3]]],{t,0,endpt},parPlotOpt]]],

{Red,PointSize[Large],Point[solLorenz[[2;;3]]/.t->Min[currentPos,endpt]]},

Dynamic[
Table[{Blue,PointSize[Large],
Tooltip[Point[{y[t],z[t]}/.eqPts[[k]]],Text@Style[Row[{"Eigenvalues: ",eigenSets[[k,1]]}],FontSize->12]]},{k,1,3}]]
},
ImageSize->150,AxesLabel->{"Y","Z"},Axes->True
]
}]}

]]
],

4,
Dynamic@
With[{solPert3=First[{x[t],y[t],z[t]}/.Quiet[NDSolve[SetPrecision[{eqLorenz,Thread[{x[0],y[0],z[0]}=={x0+perturbation*pertX,y0+perturbation*pertY,z0+perturbation*pertZ}]},Infinity],{x,y,z},{t,0,currentPos},PrecisionGoal->ControlActive[4,6],WorkingPrecision->ControlActive[MachinePrecision,20]]]],
pertPt={x0+perturbation*pertX,y0+perturbation*pertY,z0+perturbation*pertZ},
plotOpt3=Sequence[PlotPoints->ControlActive[30,30],PlotStyle->{Green,Red}],parPlotOpt2=Sequence[PlotPoints->ControlActive[80,80],MaxRecursion->ControlActive[4,4]],
optGr2D2=Sequence[Axes->True,AxesOrigin->{0,0},AspectRatio->0.8,ImageSize->200,Ticks->None]},

Row[
{
Graphics3D[(*Print["3d loop"];*){First[ParametricPlot3D[{Tooltip[solLorenz,{x0,y0,z0}],Tooltip[solPert3,pertPt]},{t,0,endpt},PlotPoints->ControlActive[100,100],MaxRecursion->ControlActive[4,4],PlotStyle->{Green,Red}]],

{Green,PointSize[Large],Point[solLorenz/.t->Evaluate@Min[currentPos,endpt]],
Red,Point[solPert3/.t->Evaluate@Min[currentPos,endpt]]},

Dynamic[Table[{Blue,PointSize[Large],Tooltip[Point[{x[t],y[t],z[t]}/.eqPts[[k]]],Text@Style[Row[{"Eigenvalues: ",eigenSets[[k,1]]}],FontSize->12]]},{k,1,Length[eqPts]}]]
},
ImageSize->300,AxesLabel->{"X","Y","Z"},Axes->True],

Graphics[First[Plot[Evaluate[{Tooltip[solLorenz[[Switch[whichGraph,"X(t)",1,"Y(t)",2,"Z(t)",3]]],{x0,y0,z0}],Tooltip[solPert3[[Switch[whichGraph,"X(t)",1,"Y(t)",2,"Z(t)",3]]],Evaluate[{x0+perturbation*pertX,y0+perturbation*pertY,z0+perturbation*pertZ}]]}],{t,0,endpt},plotOpt3]],AxesLabel->{"t",Switch[whichGraph,"X(t)","X","Y(t)","Y","Z(t)","Z"]},Axes->True,optGr2D2]
}

]]
]],

Grid[{{"The Lorenz Equations: ","Parameters:","Initial conditions:"},{"  "<>ToString@TraditionalForm[X'[t]==σ (-X[t]+Y[t])]<>"  ",Control[{{s,10,"σ"},0.1,100,0.01,ImageSize->Small,Appearance->"Labeled"}],Control[{{x0,5,"   (Subscript[X, 0])"},1,200,0.01,ImageSize->Small,Appearance->"Labeled"}]},{"  "<>ToString@TraditionalForm[Y'[t]==r*X[t]-Y[t]-X[t] Z[t]]<>"  ",Control[{{r,28,"r"},1,100,0.01,ImageSize->Small,Appearance->"Labeled"}],Control[{{y0,5,"   (Subscript[Y, 0])"},0,200,0.01,ImageSize->Small,Appearance->"Labeled"}]},{"  "<>ToString@TraditionalForm[Z'[t]==-b Z[t]+X[t] Y[t]]<>"  ",Control[{{b,8/3,"b"},0,10,0.01,ImageSize->Small,Appearance->"Labeled"}],Control[{{z0,5,"   (Subscript[Z, 0])"},1,200,0.01,ImageSize->Small,Appearance->"Labeled"}]}},Alignment->Left,Dividers->{False,{False,True,False,False,False}},ItemStyle->{Directive[Black,12,FontFamily->"Times"],None}],


Row[{Control[{{currentPos,4,"Time: "},1,40,0.01,ControlType->Animator,AnimationRepetitions->1,DefaultDuration->180,AnimationDirection->Forward,AnimationRunning->False}](*,
"Time = ",Dynamic@currentPos*)},ItemStyle->{Directive[Black,12,FontFamily->"Times"],None}],

{{ctrlPane,1,""},{
1->"Projections into coordinate planes",
2->"Lorenz attractor in 3 dimensions",
3->"Lorenz Simulation",
4->"Perturbations to initial conditions"}},

PaneSelector[
{
3->Row[{Control[{{whichView,"Time plots"},{"Time plots","Phase plots"}}],
Control[{{perturbation,0.1,"Perturbation size"},0,5,0.005}]}],
4-> Row[{Control[{{whichGraph,"X(t)"},{"X(t)","Y(t)","Z(t)"}}],
Control[{{perturbation,0.1,"Perturbation size"},0,5,0.005,ImageSize->Small,Appearance->"Labeled"}],"Direction:",
Control[{{pertX,1,"X"},{0,1},ControlType->Checkbox}],
Control[{{pertY,0,"Y"},{0,1},ControlType->Checkbox}],
Control[{{pertZ,0,"Z"},{0,1},ControlType->Checkbox}]}]
},
Dynamic[ctrlPane]],

(*Used for built-in scrolling/zoom-can reference directly or use var["computePlotrange"] for better performance during drags*)

{{plotrange,{{-5,5},{-5,5}}},ControlType->None},
{{eqPts,{}},ControlType->None},
{{localLin,{}},ControlType->None},
{eqLorenz,ControlType->None},
{{eigenSets,{}},ControlType->None},
{{solLorenz,{}},ControlType->None},
{nearestEqPt,ControlType->None},
{endpt,ControlType->None},

(*temporary holding variable*)
{var,ControlType->None},

Initialization:>(mouseMode[]:=Switch[CurrentValue/@{"ShiftKey","CommandKey"},{False,False},"special",{True,False},"pan",{False,True},"zoomIn",{True,True},"zoomOut"];

readMouse[startOrEnd_]:=(var[startOrEnd]=MousePosition["Graphics"];
var[startOrEnd<>"Scaled"]=MousePosition["GraphicsScaled"]);


Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]),
Alignment->Center,
SynchronousInitialization->False(*,SynchronousUpdating->False*)]