A. Dynamical System Study 4D

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

<<VectorFieldPlots`;

System

Parameters

{(Subscript[x, min]),(Subscript[x, max])}={,};
{(Subscript[y, min]),(Subscript[y, max])}={,};
{(Subscript[z, min]),(Subscript[z, max])}={,};

f[x_,y_,z_,u_]:=

g[x_,y_,z_,u_]:=

h[x_,y_,z_,u_]:=

i[x_,y_,z_,u_]:=

Vector field

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

Fixed Points & Nullclines

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

fp1=Sort[fp[[1]]];fp2=Sort[fp[[2]]];fp3=Sort[fp[[3]]];(*fp4=Sort[fp[[4]]]*);

{(Subscript[u, 1]),(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])}={Extract[fp1,{1}][[2]],Extract[fp1,{2}][[2]],Extract[fp1,{3}][[2]],Extract[fp1,{4}][[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[u, 2]),(Subscript[x, 2]),(Subscript[y, 2]),(Subscript[z, 2])}={Extract[fp2,{1}][[2]],Extract[fp2,{2}][[2]],Extract[fp2,{3}][[2]],Extract[fp2,{4}][[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[u, 3]),(Subscript[x, 3]),(Subscript[y, 3]),(Subscript[z, 3])}={Extract[fp3,{1}][[2]],Extract[fp3,{2}][[2]],Extract[fp3,{3}][[2]],Extract[fp3,{4}][[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])}}]];

{(Subscript[u, 4]),(Subscript[x, 4]),(Subscript[y, 4]),(Subscript[z, 4])}={Extract[fp4,{1}][[2]],Extract[fp4,{2}][[2]],Extract[fp4,{3}][[2]],Extract[fp4,{4}][[2]]};
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])}}]];

Jacobian

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

J={{a11,a12,a13,a14},{a21,a22,a23,a24},{a31,a32,a33,a34},{a41,a42,a43,a44}}

MatrixForm[J]

Fixed Points Stability - Eigenvalues

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

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

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

Phase Portrait

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

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

z'[t]==h[x[t],y[t],z[t],u[t]],z[0]==.1,

u'[t]==i[x[t],y[t],z[t],u[t]] ,u[0]==.1},{x,y,z,u},

{t,0,1000},Method->"ExplicitRungeKutta",MaxSteps->Infinity,StartingStepSize->0.0001];phaseplot=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solution],{t,10,1000},PlotPoints->5000,PlotRange->All,PlotStyle->{RGBColor[1,0,0]},ImageSize->{400,400},BoxRatios->{1,1,1},AxesLabel->{"X","Y","Z"}]

phaseplotXY=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,10,1000},PlotPoints->5000,PlotRange->All,PlotStyle->{RGBColor[1,0,0]},ImageSize->{400,400},AspectRatio->1,AxesLabel->{"X","Y"}]

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]}]

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

Show[plotX,plotY]

Show[phaseplot,fixedpoint1,fixedpoint2,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]};

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

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

flowmanifold=ContourPlot3D[ϕ==0,{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},{z,(Subscript[z, min]),(Subscript[z, max])},Mesh->None,ContourStyle->Directive[Magenta,Opacity[0.3],Specularity[White,30]],PlotPoints->{7,7,7},BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"},BoxRatios->{1,1,1}]

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

Flow Curvature Invariance

Normal Vector

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

Invariance

Ψ=FullSimplify[V . η]

invariantpart=ContourPlot3D[Ψ==0,{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},{z,(Subscript[z, min]),(Subscript[z, max])},Mesh->None,ContourStyle->Directive[Red,Opacity[0.3],Specularity[White,30]],PlotPoints->{7,7,7},BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"},BoxRatios->{1,1,1}]

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

Flow

tableplot={};
For[t=1800,t≤2000,t=t+ .5,
x=Evaluate[{x[t],y[t],z[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t],z[t]}/.solution][[1]][[2]];
z=Evaluate[{x[t],y[t],z[t]}/.solution][[1]][[3]];
vx=f[x,y,z];
vy=g[x,y,z];
vz=h[x,y,z];
tmptable={{x,y,z},{vx,vy,vz}};
tableplot=Append[tableplot,tmptable];
];
flow=ListPlotVectorField3D[tableplot ,ScaleFunction->(.025&),ScaleFactor->2,VectorHeads->True,ColorFunction->(Hue[0.66]&)]