Non-singularly perturbed systems: Lorenz Model (Ch 11, page 197)

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

<<VectorFieldPlots`;

System

Parameters

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

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

Over-Acceleration (jerk)

Γ=(J . γ +Dt[J] . V)/.{Dt[x]->f[x,y,z],Dt[y]->g[x,y,z],Dt[z]->h[x,y,z],Dt[parameters]->0}

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]

Flow Curvature Manifold

ϕ=Factor[Γ  . Cross[V,γ]]

flowmanifold=ContourPlot3D[ϕ==0,{x,-25,25},{y,-30,30},{z,0,50},PlotPoints->{10,10,10},BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"},AspectRatio->1];

Show[phaseplot,fixedpoint1,fixedpoint2,fixedpoint3,flowmanifold,BoxRatios->{1,1,1},ImageSize->{600,600},Boxed->False,Axes->False]

Flow Curvature Invariance

Normal vector

η={D[ϕ,x],D[ϕ,y],D[ϕ,z]}

Invariance

Ψ=Factor[(V . η)]

invariantpart=ContourPlot3D[ϕ==0,{x,-25,25},{y,-30,30},{z,0,50},PlotPoints->{10,10,10},BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"},AspectRatio->1];

Show[phaseplot,fixedpoint1,fixedpoint2,fixedpoint3,invariantpart,BoxRatios->{1,1,1},ImageSize->{600,600},Boxed->False,Axes->False]