Pikovskii-Rabinovich-Trakhtengerts model (Ch 14, page 219)

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

<<VectorFieldPlots`;

System

Parameters

k=6; (Subscript[ν, 1])=1; (Subscript[ν, 2])=4;(Subscript[ν, 3])=1;

{(Subscript[x, min]),(Subscript[x, max])}={-5,5};
{(Subscript[y, min]),(Subscript[y, max])}={-12,12};
{(Subscript[z, min]),(Subscript[z, max])}={0,10};

f[x_,y_,z_]:=k y - (Subscript[ν, 1])x -y z

g[x_,y_,z_]:=k x - (Subscript[ν, 2])y + x z

h[x_,y_,z_]:= x y - (Subscript[ν, 3])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])}]//N

fp1=Sort[fp[[1]]];fp2=Sort[fp[[4]]];fp3=Sort[fp[[5]]];

{(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,500},Method->"ExplicitRungeKutta",MaxSteps->Infinity,StartingStepSize->0.001];phaseplot=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solution],{t,400,500},PlotPoints->2000,PlotRange->All,PlotStyle->{Red,Thickness[0.001]},ImageSize->{400,400},BoxRatios->{1,1,1},AxesLabel->{"X","Y","Z"}]

Show[phaseplot,fixedpoint1,fixedpoint2,fixedpoint3]

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

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

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

Show[plotX,plotY]

Flow Curvature Manifold

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

flowmanifold=ContourPlot3D[ϕ==0,{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},{z,(Subscript[z, min]),(Subscript[z, max])},PlotPoints->{50,50,50},BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"}]

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

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]

MatrixForm[H]

Fixed point O (0, 0, 0)

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

FullSimplify[Factor[((ϵ^3) (-4.491072818385036`*^13+ϵ (5.024020797848505`*^13+ϵ (1.973695610880331`*^14+ϵ (1.6031504439706162`*^14+ϵ (-1.505602846310482`*^14+ϵ (-3.648009919672299`*^14+ϵ (1.5778817935871469`*^13+ϵ (2.1479393590464012`*^14-1.161216`*^6 (-21.886285616587216`+ϵ) (-3.7727364315103915`+ϵ) (-0.23637646216748398`+ϵ) ϵ (3.307201230918598`+ϵ) (16.775597041897242`+(-7.236743258517178`+ϵ) ϵ) (20.52239485893051`+(-6.867319003013444`+ϵ) ϵ) (13.572876528225972`+(-4.591698686889706`+ϵ) ϵ) (12.197792278673223`+ϵ (2.568361097256957`+ϵ)) (2.9210725268154074`+ϵ (2.652105067017786`+ϵ)))))))))))/((-1.4176813637503598`+ϵ) (15.019838529989403`+(-6.35395560793376`+ϵ) ϵ) (3456.` (ϵ^3)+1452.` (ϵ^4)-2478.` (ϵ^5)-804.` (ϵ^6)+588.` (ϵ^7)-108.` (ϵ^8)) (1.5028169479131313`+ϵ (2.327192527239676`+ϵ)))]/.{ϵ->0}]

Hessian is positive. Then O is a saddle-node.

Fixed point (SubMinus[I])

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

FullSimplify[Factor[(5.04077372765966`*^-38 (958581.`-1.407955497567615`*^21 ϵ-7.41687984230189`*^35 (ϵ^2)-6.422570016319116`*^49 (ϵ^3)+3.669401671092041`*^50 (ϵ^4)-8.345898676839894`*^50 (ϵ^5)+1.0470778356706771`*^51 (ϵ^6)-7.779442834223173`*^50 (ϵ^7)+3.6346563059804175`*^50 (ϵ^8)-1.0389253797126667`*^50 (ϵ^9)+1.099642382844937`*^49 (ϵ^10)+2.9019988908888147`*^48 (ϵ^11)-9.373275748102487`*^47 (ϵ^12)+6.526577903829689`*^46 (ϵ^13)+2.8470927754562774`*^44 (ϵ^14)+3.8032816255757364`*^41 (ϵ^15)+3.332821687237274`*^39 (ϵ^16)))/((-2.9306683067442765`+ϵ) (-0.8116351884645938`+ϵ) (-0.00004033724301898326`+ϵ) (138.9678509656974`+ϵ) (1.627042006913238`*^-9+0.000040337242948004715` ϵ+(ϵ^2)))]/.{ϵ->0}]

Hessian is negative. Then (SubMinus[I]) is a focus.

Fixed point (SubPlus[I])

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

FullSimplify[Factor[(-3.8999101264562414`*^13+ϵ (-1.6186166222262816`*^14-167.99999999999997` ϵ (6.598024378479714`+ϵ) (20.67208944081457`+ϵ) (41.08825248219877`+(-7.431565924907058`+ϵ) ϵ) (3.844348445783756`+(-0.24032362746951685`+ϵ) ϵ) (8.407523925116479`+ϵ (4.371349014561617`+ϵ)) (32.26469998355986`+ϵ (11.268728219984522`+ϵ)) (200.17226691197163`+ϵ (27.369787009409784`+ϵ)) (1425.1742043294985`+ϵ (56.29278402751161`+ϵ))))/((1.0000402728954167`+ϵ) (4.615490948418588`+ϵ) (16.08233080796594`+ϵ) (27.557012239554513`+ϵ) (0.9999597287264405`+ϵ (1.999959727104609`+ϵ)))]/.{ϵ->0}]

Hessian is negative. Then (SubPlus[I]) is a focus.