Example 8.3: Lorenz Model (Ch 8, page 147)

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 . γ +0Dt[J] . V)/.{Dt[x]->f[x,y,z],Dt[y]->g[x,y,z],Dt[z]->h[x,y,z],Dt[parameters]->0}

Flow Curvature Manifold

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

Osculating Plane

X={6 (Sqrt[2])+x,6 (Sqrt[2])+y,-27+z}

Π=Factor[X  . Cross[V,γ]]

X={x,y,z}

Π=Factor[X  . Cross[V,γ]]

X={-6 (Sqrt[2])+x,-6 (Sqrt[2])+y,-27+z}

Π=Factor[X  . Cross[V,γ]]

Hessian Theorem 8.3

a11=Simplify[D[ϕ,x,x]]

a12=Simplify[D[ϕ,x,y]]

a13=Simplify[D[ϕ,x,z]]

a21=Simplify[D[ϕ,y,x]]

a22=Simplify[D[ϕ,y,y]]

a23=Simplify[D[ϕ,y,z]]

a31=Simplify[D[ϕ,z,x]]

a32=Simplify[D[ϕ,z,y]]

a33=Simplify[D[ϕ,z,z]]

Matrice hessienne

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

MatrixForm[H]

K=Det[H]

Fixed point (SubMinus[I])

Factor[(K/ϕ)/.{x->-6 (Sqrt[2])+ϵ,y->-6 (Sqrt[2])+ϵ,z->27+ϵ}]

-(200 (-111480393626025984-96881770651189248 (Sqrt[2])+167174835583057920 ϵ+97002513644912640 (Sqrt[2]) ϵ-64872746025615360 (ϵ^2)-57834947731783680 (Sqrt[2]) (ϵ^2)+23499502759411712 (ϵ^3)+11748868500111360 (Sqrt[2]) (ϵ^3)-2594672117624832 (ϵ^4)-3350475545495552 (Sqrt[2]) (ϵ^4)+732614203038720 (ϵ^5)+177596576389120 (Sqrt[2]) (ϵ^5)-13100574501664 (ϵ^6)-63029140567360 (Sqrt[2]) (ϵ^6)+8445770487440 (ϵ^7)-17192307432 (Sqrt[2]) (ϵ^7)+59952266206 (ϵ^8)-431363143656 (Sqrt[2]) (ϵ^8)+32826519501 (ϵ^9)-3362869296 (Sqrt[2]) (ϵ^9)+192088755 (ϵ^10)-904607712 (Sqrt[2]) (ϵ^10)+34425996 (ϵ^11)-2945376 (Sqrt[2]) (ϵ^11)+38304 (ϵ^12)-408240 (Sqrt[2]) (ϵ^12)+4536 (ϵ^13)))/(9 (-16128-14016 (Sqrt[2])+3680 ϵ+1344 (Sqrt[2]) ϵ-71 (ϵ^2)-180 (Sqrt[2]) (ϵ^2)+6 (ϵ^3)))/.{ϵ->0}//N

Fixed point O (0, 0, 0)

FullSimplify[(K/ϕ)/.{x->0+ϵ,y->0+ϵ,z->0+ϵ}]

-(200 (-8154695373037056000+ϵ (1949339191261824000+ϵ (76905877814784000+ϵ (-193276265731787520+ϵ (-10488489197406912+ϵ (-11690309960110416+ϵ (4031622516163440+ϵ (-519951015582240+ϵ (34591993991084+ϵ (-1158617462428+ϵ (13395439528+3 ϵ (-37168492+ϵ (13169361+8 ϵ (-26459+21 ϵ (-167+9 ϵ))))))))))))))))/(9 (1704240+ϵ (203130+ϵ (-46314+ϵ (4868+ϵ (-233+6 ϵ))))))/.{ϵ->0}//N

Fixed point (SubPlus[I])

Factor[(K/ϕ)/.{x->6 (Sqrt[2])+ϵ,y->6 (Sqrt[2])+ϵ,z->27+ϵ}]

-(200 (-111480393626025984+96881770651189248 (Sqrt[2])+167174835583057920 ϵ-97002513644912640 (Sqrt[2]) ϵ-64872746025615360 (ϵ^2)+57834947731783680 (Sqrt[2]) (ϵ^2)+23499502759411712 (ϵ^3)-11748868500111360 (Sqrt[2]) (ϵ^3)-2594672117624832 (ϵ^4)+3350475545495552 (Sqrt[2]) (ϵ^4)+732614203038720 (ϵ^5)-177596576389120 (Sqrt[2]) (ϵ^5)-13100574501664 (ϵ^6)+63029140567360 (Sqrt[2]) (ϵ^6)+8445770487440 (ϵ^7)+17192307432 (Sqrt[2]) (ϵ^7)+59952266206 (ϵ^8)+431363143656 (Sqrt[2]) (ϵ^8)+32826519501 (ϵ^9)+3362869296 (Sqrt[2]) (ϵ^9)+192088755 (ϵ^10)+904607712 (Sqrt[2]) (ϵ^10)+34425996 (ϵ^11)+2945376 (Sqrt[2]) (ϵ^11)+38304 (ϵ^12)+408240 (Sqrt[2]) (ϵ^12)+4536 (ϵ^13)))/(9 (-16128+14016 (Sqrt[2])+3680 ϵ-1344 (Sqrt[2]) ϵ-71 (ϵ^2)+180 (Sqrt[2]) (ϵ^2)+6 (ϵ^3)))/.{ϵ->0}//N

Corollary 8.2

Fixed point (SubMinus[I])

FullSimplify[(ϕ/Π)/.{x->-6 (Sqrt[2])+ϵ,y->-6 (Sqrt[2])+ϵ,z->27+ϵ}]

((10 (6 (Sqrt[2])-ϵ) (192 (84+73 (Sqrt[2]))+ϵ (-32 (115+42 (Sqrt[2]))+(71+180 (Sqrt[2])-6 ϵ) ϵ)))/(-16 (73+42 (Sqrt[2]))+(71+108 (Sqrt[2])-6 ϵ) ϵ))/.{ϵ->0}//N

Fixed point (SubPlus[I])

FullSimplify[(ϕ/Π)/.{x->6 (Sqrt[2])+ϵ,y->6 (Sqrt[2])+ϵ,z->27+ϵ}]

-((10 (6 (Sqrt[2])+ϵ) (192 (-84+73 (Sqrt[2]))+ϵ (3680-1344 (Sqrt[2])+ϵ (-71+180 (Sqrt[2])+6 ϵ))))/(1168-672 (Sqrt[2])+ϵ (-71+108 (Sqrt[2])+6 ϵ)))/.{ϵ->0}//N