Mofatt model 5D (Ch 17, page 248)

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

<<VectorFieldPlots`;

System

Parameters

α=1.5;ω=1;η=4;κ=1;

ξ=100;

{(Subscript[x, min]),(Subscript[x, max])}={0,10};
{(Subscript[y, min]),(Subscript[y, max])}={0,10};
{(Subscript[z, min]),(Subscript[z, max])}={3,13};

f[x_,y_,z_,u_,w_]:=α(-η x + ω y z)

g[x_,y_,z_,u_,w_]:=-η y + ω x z

h[x_,y_,z_,u_,w_]:=κ(u - z - x y)

i[x_,y_,z_,u_,w_]:= - u + ξ z - w z

j[x_,y_,z_,u_,w_]:=- w + u z

Vector field

V={f[x,y,z,u,w],g[x,y,z,u,w],h[x,y,z,u,w],i[x,y,z,u,w],j[x,y,z,u,w]}

Fixed Points

(Subscript[x, 1])=.
(Subscript[y, 1])=.
(Subscript[z, 1])=.
(Subscript[u, 1])=.
(Subscript[w, 1])=.

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

fp1=Sort[fp[[1]]];fp2=Sort[fp[[2]]];fp3=Sort[fp[[3]]];fp4=Sort[fp[[4]]];fp5=Sort[fp[[5]]];fp6=Sort[fp[[6]]];fp7=Sort[fp[[7]]];

{(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])}={0,0,0};
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])}={-4.419209404939518`,4.419209404939518`,-4};
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])}={4.419209404939518`,4.419209404939518`,4};
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])}}]];

{(Subscript[x, 4]),(Subscript[y, 4]),(Subscript[z, 4])}={-4.419209404939518`,-4.419209404939518`,4 };
fixedpoint4=Show[Graphics3D[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 4]),(Subscript[y, 4]),(Subscript[z, 4])}]},PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])},{(Subscript[z, min]),(Subscript[z, max])}}]];

{(Subscript[x, 5]),(Subscript[y, 5]),(Subscript[z, 5])}={4.419209404939518`,-4.419209404939518`,-4};
fixedpoint5=Show[Graphics3D[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 4]),(Subscript[y, 4]),(Subscript[z, 4])}]},PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])},{(Subscript[z, min]),(Subscript[z, max])}}]];

{(Subscript[x, 6]),(Subscript[y, 6]),(Subscript[z, 6])}={0.,0.,-9.9498743710662`};
fixedpoint6=Show[Graphics3D[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 4]),(Subscript[y, 4]),(Subscript[z, 4])}]},PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])},{(Subscript[z, min]),(Subscript[z, max])}}]];

{(Subscript[x, 7]),(Subscript[y, 7]),(Subscript[z, 7])}={0.,0.,9.9498743710662`};
fixedpoint7=Show[Graphics3D[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 4]),(Subscript[y, 4]),(Subscript[z, 4])}]},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,u,w],x]]

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

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

a14=Simplify[D[f[x,y,z,u,w],u]]

a15=Simplify[D[f[x,y,z,u,w],w]]

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

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

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

a24=Simplify[D[g[x,y,z,u,w],u]]

a25=Simplify[D[g[x,y,z,u,w],w]]

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

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

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

a34=Simplify[D[h[x,y,z,u,w],u]]

a35=Simplify[D[h[x,y,z,u,w],w]]

a41=Simplify[D[i[x,y,z,u,w],x]]

a42=Simplify[D[i[x,y,z,u,w],y]]

a43=Simplify[D[i[x,y,z,u,w],z]]

a44=Simplify[D[i[x,y,z,u,w],u]]

a45=Simplify[D[i[x,y,z,u,w],w]]

a51=Simplify[D[j[x,y,z,u,w],x]]

a52=Simplify[D[j[x,y,z,u,w],y]]

a53=Simplify[D[j[x,y,z,u,w],z]]

a54=Simplify[D[j[x,y,z,u,w],u]]

a55=Simplify[D[j[x,y,z,u,w],w]]

J={{a11,a12,a13,a14,a15},{a21,a22,a23,a24,a25},{a31,a32,a33,a34,a35},{a41,a42,a43,a44,a45},{a51,a52,a53,a54,a55}}

MatrixForm[J]

Fixed Points Stability - Eigenvalues

Eigenvalues[J]/.{x-> 0,y->0,z->0,u->0,w->0}//N

Eigenvalues[J]/.{w->94.11764705882354`,u->-23.529411764705884`,x->-4.419209404939518`,y->4.419209404939518`,z->-4.}//N

Eigenvalues[J]/.{w->94.11764705882354`,u->-23.529411764705884`,x->4.419209404939518`,y->-4.419209404939518`,z->-4.}//N

Eigenvalues[J]/.{w->94.11764705882354`,u->23.529411764705884`,x->-4.419209404939518`,y->-4.419209404939518`,z->4.}//Chop

Eigenvalues[J]/.{w->94.11764705882354`,u->23.529411764705884`,x->4.419209404939518`,y->4.419209404939518`,z->4.}//N

Eigenvalues[J]/.{w->99.`,u->-9.9498743710662`,x->0.,y->0.,z->-9.9498743710662`}//N

Eigenvalues[J]/.{w->99.`,u->9.9498743710662`,x->0.,y->0.,z->9.9498743710662`}//N

Phase Portrait

solution = NDSolve[{x'[t] == f[x[t], y[t], z[t], u[t], w[t]], x[0] == 1,
   
   y'[t] == g[x[t], y[t], z[t], u[t], w[t]], y[0] == 1,
   
   z'[t] == h[x[t], y[t], z[t], u[t], w[t]], z[0] == 1,
   
   u'[t] == i[x[t], y[t], z[t], u[t], w[t]], u[0] == 1,
   
   w'[t] == j[x[t], y[t], z[t], u[t], w[t]], w[0] == 1}, {x, y, z, u, w},
  
  {t, 0, 300}, Method -> "ExplicitRungeKutta", MaxSteps -> Infinity, StartingStepSize -> 0.0001]; phaseplot = ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. solution], {t, 250, 300}, PlotPoints -> 6000, PlotRange -> All, PlotStyle -> {RGBColor[1, 0, 0]}, ImageSize -> {400, 400}, BoxRatios -> {1, 1, 1}, AxesLabel -> {"X", "Y", "Z"}]

phaseplotXZ=ParametricPlot[Evaluate[{x[t],z[t]}/.solution],{t,250,300},PlotPoints->5000,PlotRange->All,PlotStyle->{RGBColor[1,0,0]},ImageSize->{400,400},AspectRatio->1,AxesLabel->{"X","Z"}]

plotX=Plot[Evaluate[x[t]/.solution],{t,250,300},PlotRange->{-10,10},PlotStyle->{RGBColor[0,0,1]}]

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

plotZ=Plot[Evaluate[z[t]/.solution],{t,250,300},PlotRange->{-15,15},PlotStyle->{RGBColor[1,0,0]}]

plotU=Plot[Evaluate[u[t]/.solution],{t,250,300},PlotRange->{-40,40},PlotStyle->{Cyan,Thickness[0.005]}]

plotW=Plot[Evaluate[w[t]/.solution],{t,250,300},PlotRange->{40,140},PlotStyle->{Pink,Thickness[0.005]}]

Show[plotX,plotY,plotZ]

Show[phaseplot,fixedpoint1,fixedpoint2,fixedpoint3,ImageSize->{600,600}]

Flow Curvature Manifold

V1={f[x,y,z,u,w],g[x,y,z,u,w],h[x,y,z,u,w],i[x,y,z,u,w],j[x,y,z,u,w]}

V2=Dt[V1]/.{Dt[x]->f[x,y,z,u,w],Dt[y]->g[x,y,z,u,w], Dt[z]->h[x,y,z,u,w],Dt[u]->i[x,y,z,u,w],Dt[w]->j[x,y,z,u,w]};

V3=Dt[V2]/.{Dt[x]->f[x,y,z,u,w],Dt[y]->g[x,y,z,u,w], Dt[z]->h[x,y,z,u,w],Dt[u]->i[x,y,z,u,w],Dt[w]->j[x,y,z,u,w]};

V4=Dt[V3]/.{Dt[x]->f[x,y,z,u,w],Dt[y]->g[x,y,z,u,w], Dt[z]->h[x,y,z,u,w],Dt[u]->i[x,y,z,u,w],Dt[w]->j[x,y,z,u,w]};

V5=Dt[V4]/.{Dt[x]->f[x,y,z,u,w],Dt[y]->g[x,y,z,u,w], Dt[z]->h[x,y,z,u,w],Dt[u]->i[x,y,z,u,w],Dt[w]->j[x,y,z,u,w]};

K={V1,V2,V3,V4,V5};

ϕ=Factor[Det[K]/.{u->0,w->0}];

In high dimensional (hyper)space phase the hypersurface has a very "complex" structure which may hide the attractor. So, let's plot remove the nonstationary part of the flow curvature manifold in order to improve the figure.

V1={f[x,y,z,u,w],g[x,y,z,u,w],h[x,y,z,u,w],i[x,y,z,u,w],j[x,y,z,u,w]}

V2=J . V1;

V3=J . V2;

V4=J . V3;

V5=J . V4;

K={V1,V2,V3,V4,V5};

ϕ=Factor[Det[K]/.{u->0,w->0}]//N

flowmanifoldL=ContourPlot3D[ϕ==0,{x,1,10},{y,1,10},{z,1,11},PlotPoints->{10,10,10},BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"}]

flowmanifoldR=ContourPlot3D[ϕ==0,{x,-10,-1},{y,1,10},{z,-11,-1},PlotPoints->{10,10,10},BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"}]

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