Lorenz Model

Clear["Global`*"]
Off[General::spell1]

System

Parameters

σ=10;
r=28;
β=(8/3);

{(Subscript[x, min]),(Subscript[x, max])}={-15,15};
{(Subscript[y, min]),(Subscript[y, max])}={-25,25};
{(Subscript[z, min]),(Subscript[z, max])}={0,50};

f[x_,y_,z_]:=σ(y- x)

g[x_,y_,z_]:=-x z+r x-y

h[x_,y_,z_]:=x y-β z

Vector field

V={f[x,y,z],g[x,y,z],h[x,y,z]}

Fixed Points

fp=Solve[{f[(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])]==0,g[(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])]==0,h[(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])]==0},{(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])}]

fp1=Sort[fp[[1]]];fp2=Sort[fp[[2]]];fp3=Sort[fp[[3]]];

{(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])}={Extract[fp1,{1}][[2]],Extract[fp1,{2}][[2]],Extract[fp1,{3}][[2]]};
fixedpoint1=Show[Graphics3D[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])}]},PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])},{(Subscript[z, min]),(Subscript[z, max])}}]];

{(Subscript[x, 2]),(Subscript[y, 2]),(Subscript[z, 2])}={Extract[fp2,{1}][[2]],Extract[fp2,{2}][[2]],Extract[fp2,{3}][[2]]};
fixedpoint2=Show[Graphics3D[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 2]),(Subscript[y, 2]),(Subscript[z, 2])}]},PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])},{(Subscript[z, min]),(Subscript[z, max])}}]];

{(Subscript[x, 3]),(Subscript[y, 3]),(Subscript[z, 3])}={Extract[fp3,{1}][[2]],Extract[fp3,{2}][[2]],Extract[fp3,{3}][[2]]};
fixedpoint3=Show[Graphics3D[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 3]),(Subscript[y, 3]),(Subscript[z, 3])}]},PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])},{(Subscript[z, min]),(Subscript[z, max])}}]];

Jacobian

a11=Simplify[D[f[x,y,z],x]]

a12=Simplify[D[f[x,y,z],y]]

a13=Simplify[D[f[x,y,z],z]]

a21=Simplify[D[g[x,y,z],x]]

a22=Simplify[D[g[x,y,z],y]]

a23=Simplify[D[g[x,y,z],z]]

a31=Simplify[D[h[x,y,z],x]]

a32=Simplify[D[h[x,y,z],y]]

a33=Simplify[D[h[x,y,z],z]]

J={{a11,a12,a13},{a21,a22,a23},{a31,a32,a33}}

MatrixForm[J]

Acceleration

γ=J . V

Fixed Points Stability - Eigenvalues

Eigenvalues[J]/.{x-> (Subscript[x, 1]),y->(Subscript[y, 1]),z->(Subscript[z, 1])}//N

Eigenvalues[J]/.{x-> (Subscript[x, 2]),y->(Subscript[y, 2]),z->(Subscript[z, 2])}//N

Eigenvalues[J]/.{x-> (Subscript[x, 3]),y->(Subscript[y, 3]),z->(Subscript[z, 3])}//N

Phase Portrait

solution=NDSolve[{x'[t]==f[x[t],y[t],z[t]],x[0]==1,

y'[t]==g[x[t],y[t],z[t]],y[0]==1,

z'[t]==h[x[t],y[t],z[t]],z[0]==1},{x,y,z},

{t,0,100},Method->"ExplicitRungeKutta",MaxSteps->Infinity,StartingStepSize->0.001];phaseplot=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solution],{t,50,100},PlotPoints->2000,PlotRange->All,PlotStyle->{Red,Thickness[0.0007]},ImageSize->{400,400},AspectRatio->1,AxesLabel->{"X","Y","Z"}];

Show[phaseplot,fixedpoint1,fixedpoint2,fixedpoint3]

plotX=Plot[Evaluate[x[t]/.solution],{t,0,60},PlotRange->{(Subscript[x, min]),(Subscript[x, max])},PlotStyle->{RGBColor[0,0,1]}]

plotY=Plot[Evaluate[y[t]/.solution],{t,0,60},PlotRange->{(Subscript[y, min]),(Subscript[y, max])},PlotStyle->{RGBColor[0,1,0]}]

plotZ=Plot[Evaluate[z[t]/.solution],{t,0,60},PlotRange->{(Subscript[z, min]),(Subscript[z, max])},PlotStyle->{RGBColor[1,0,0]}]

Show[plotX,plotY]

Curvature equations (relations)

Cross[V,γ]

A=Factor[784 (x^3)-((308 (x^2) y)/3)+((38 x (y^2))/3)+(x^3) (y^2)-10 (y^3)-((5600 x z)/9)-56 (x^3) z+((6680 y z)/9)-(5/3) (x^2) y z+((88 x (z^2))/3)+(x^3) (z^2)-((80 y (z^2))/3)];

B=Factor[280 (x^3)-((110 (x^2) y)/3)+((380 x (y^2))/3)-100 (y^3)-((8480 x z)/9)-10 (x^3) z+((2000 y z)/9)+((80 x (z^2))/3)];

G=-7560 (x^2)-2430 x y+10 (x^3) y+2700 (y^2)-10 (x^2) (y^2)+((1570 (x^2) z)/3)+((350 x y z)/3)-100 (y^2) z-10 (x^2) (z^2);

Wilkinson method

Normal vector to S1

η1={D[A,x],D[A,y],D[A,z]}

Normal vector to S2

η2={D[G,x],D[G,y],D[G,z]}

Cross[η1,η2]

Expand[{-1512000 (x^2)-((550760 (x^4))/3)+560 (x^6)+5163600 x y+((839290 (x^3) y)/3)-((170 (x^5) y)/3)-4008000 (y^2)+((64100 (x^2) (y^2))/3)+200 (x^4) (y^2)-((18100 x (y^3))/3)-200 (x^3) (y^3)+3000 (y^4)+((5432240 (x^2) z)/9)+((107890 (x^4) z)/9)-20 (x^6) z-((5137600 x y z)/9)-20800 (x^3) y z+((3260000 (y^2) z)/9)-(1900/3) (x^2) (y^2) z-((320800 (x^2) (z^2))/9)-200 (x^4) (z^2)+(133600/9) x y (z^2)+400 (x^3) y (z^2)-8000 (y^2) (z^2)+((1600 (x^2) (z^3))/3),9408000 (x^2)-384160 (x^4)-9710400 x y-((218960 (x^3) y)/9)-1680 (x^5) y-1803600 (y^2)+((2621590 (x^2) (y^2))/9)-500 (x^4) (y^2)-((331700 x (y^3))/9)-((950 (x^3) (y^3))/3)+((3800 (y^4))/3)+300 (x^2) (y^4)-((32742080 (x^2) z)/27)+((138320 (x^4) z)/3)+((38898880 x y z)/27)+5860 (x^3) y z+60 (x^5) y z+((4157200 (y^2) z)/27)-(172130/9) (x^2) (y^2) z+20 (x^4) (y^2) z+(2200/3) x (y^3) z+((138160 (x^2) (z^2))/3)-((5150 (x^4) (z^2))/3)-(605200/9) x y (z^2)-150 (x^3) y (z^2)-((29600 (y^2) (z^2))/9)+300 (x^2) (y^2) (z^2)-((1760 (x^2) (z^3))/3)+20 (x^4) (z^3)+(3200/3) x y (z^3),-7267680 (x^3)+23520 (x^5)+13333320 (x^2) y-((47320 (x^4) y)/3)-1531620 x (y^2)-1010 (x^3) (y^2)-30 (x^5) (y^2)-4500 (y^3)+((52060 (x^2) (y^3))/3)-20 (x^4) (y^3)-600 x (y^4)+12734400 x z+((6828080 (x^3) z)/9)-1680 (x^5) z-1556400 y z-(12796790/9) (x^2) y z+(3850/3) (x^4) y z+(203500/3) x (y^2) z+150 (x^3) (y^2) z+((2900 (y^3) z)/3)-600 (x^2) (y^3) z-((35746160 x (z^2))/27)-((242150 (x^3) (z^2))/9)+30 (x^5) (z^2)+((3549200 y (z^2))/27)+(452930/9) (x^2) y (z^2)-20 (x^4) y (z^2)-(1400/3) x (y^2) (z^2)+((415600 x (z^3))/9)+((950 (x^3) (z^3))/3)-((24800 y (z^3))/9)-600 (x^2) y (z^3)-((1600 x (z^4))/3)}/(1×(10^6)) ]//N

Associated dynamical system

f[x_,y_,z_]:=-1.512` (x^2)-0.18358666666666668` (x^4)+0.00056` (x^6)+5.1636` x y+0.2797633333333333` (x^3) y-0.000056666666666666664` (x^5) y-4.008` (y^2)+0.021366666666666666` (x^2) (y^2)+0.0002` (x^4) (y^2)-0.006033333333333333` x (y^3)-0.0002` (x^3) (y^3)+0.003` (y^4)+0.6035822222222222` (x^2) z+0.011987777777777778` (x^4) z-0.00002` (x^6) z-0.5708444444444445` x y z-0.0208` (x^3) y z+0.3622222222222222` (y^2) z-0.0006333333333333333` (x^2) (y^2) z-0.035644444444444444` (x^2) (z^2)-0.0002` (x^4) (z^2)+0.014844444444444445` x y (z^2)+0.0004` (x^3) y (z^2)-0.008` (y^2) (z^2)+0.0005333333333333334` (x^2) (z^3)

g[x_,y_,z_]:=9.408` (x^2)-0.38416` (x^4)-9.7104` x y-0.02432888888888889` (x^3) y-0.00168` (x^5) y-1.8036` (y^2)+0.29128777777777776` (x^2) (y^2)-0.0005` (x^4) (y^2)-0.03685555555555556` x (y^3)-0.00031666666666666665` (x^3) (y^3)+0.0012666666666666666` (y^4)+0.0003` (x^2) (y^4)-1.2126696296296295` (x^2) z+0.046106666666666664` (x^4) z+1.4406992592592593` x y z+0.00586` (x^3) y z+0.00006` (x^5) y z+0.15397037037037037` (y^2) z-0.019125555555555555` (x^2) (y^2) z+0.00002` (x^4) (y^2) z+0.0007333333333333333` x (y^3) z+0.046053333333333335` (x^2) (z^2)-0.0017166666666666667` (x^4) (z^2)-0.06724444444444444` x y (z^2)-0.00015` (x^3) y (z^2)-0.003288888888888889` (y^2) (z^2)+0.0003` (x^2) (y^2) (z^2)-0.0005866666666666667` (x^2) (z^3)+0.00002` (x^4) (z^3)+0.0010666666666666667` x y (z^3)

h[x_,y_,z_]:=-7.26768` (x^3)+0.02352` (x^5)+13.33332` (x^2) y-0.015773333333333334` (x^4) y-1.53162` x (y^2)-0.00101` (x^3) (y^2)-0.00003` (x^5) (y^2)-0.0045` (y^3)+0.01735333333333333` (x^2) (y^3)-0.00002` (x^4) (y^3)-0.0006` x (y^4)+12.7344` x z+0.7586755555555555` (x^3) z-0.00168` (x^5) z-1.5564` y z-1.4218655555555555` (x^2) y z+0.0012833333333333334` (x^4) y z+0.06783333333333333` x (y^2) z+0.00015` (x^3) (y^2) z+0.0009666666666666667` (y^3) z-0.0006` (x^2) (y^3) z-1.323931851851852` x (z^2)-0.026905555555555554` (x^3) (z^2)+0.00003` (x^5) (z^2)+0.13145185185185185` y (z^2)+0.050325555555555554` (x^2) y (z^2)-0.00002` (x^4) y (z^2)-0.00046666666666666666` x (y^2) (z^2)+0.046177777777777776` x (z^3)+0.00031666666666666665` (x^3) (z^3)-0.0027555555555555554` y (z^3)-0.0006` (x^2) y (z^3)-0.0005333333333333334` x (z^4)

solution=NDSolve[{x'[t]==f[x[t],y[t],z[t]],x[0]==.001,

y'[t]==g[x[t],y[t],z[t]],y[0]==.01,

z'[t]==h[x[t],y[t],z[t]],z[0]==1},{x,y,z},

{t,0,40},Method->"ExplicitRungeKutta",MaxSteps->Infinity,StartingStepSize->0.001];phaseplotBobD=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solution],{t,20,40},PlotPoints->2000,PlotRange->All,PlotStyle->{Blue,Thickness[0.007]},ImageSize->{400,400},BoxRatios->{1,1,1},AxesLabel->{"X","Y","Z"}]

solution=NDSolve[{x'[t]==-f[x[t],y[t],z[t]],x[0]==-.001,

y'[t]==-g[x[t],y[t],z[t]],y[0]==-.01,

z'[t]==-h[x[t],y[t],z[t]],z[0]==.00001},{x,y,z},

{t,0,40},Method->"ExplicitRungeKutta",MaxSteps->Infinity,StartingStepSize->0.001];phaseplotBobG=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solution],{t,0,40},PlotPoints->2000,PlotRange->All,PlotStyle->{Blue,Thickness[0.007]},ImageSize->{400,400},BoxRatios->{1,1,1},AxesLabel->{"X","Y","Z"}]

Show[phaseplot,phaseplotBobG,phaseplotBobD,fixedpoint1,fixedpoint2,fixedpoint3]