Chua's cubic model 5D (Ch 17, page 245)

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

<<VectorFieldPlots`;

System

α=11.934;
β=14.47;
γ=406.5;
(Subscript[γ, 1])=-0.0152;
(Subscript[γ, 2])=41000;
γ=406.5;

a=-1.246;
b=-0.6724;

d=1.75;

(Subscript[c, 1])=-((35((((d^2)-1))^2)(a-b))/(16(d^7)));
(Subscript[c, 2])=(((45(d^4)-50(d^2)+21)(a-b))/(16(d^5)))+b+1;

Nonlinear (cubic) function

F[x_ ]=(Subscript[c, 1]) (x^3)+(Subscript[c, 2])x;
Plot[F[x],{x,-2,2}]

Parameters

{(Subscript[x, min]),(Subscript[x, max])}={-2.5,2.5};
{(Subscript[y, min]),(Subscript[y, max])}={-.5,.5};
{(Subscript[z, min]),(Subscript[z, max])}={-3.5,3.5};

f[x_,y_,z_,u_,w_]:=α(y -  ((Subscript[c, 1]) (x^3)+(Subscript[c, 2])x))

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

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

i[x_,y_,z_,u_,w_]:=-γ(z + w)

j[x_,y_,z_,u_,w_]:=(Subscript[γ, 2])(u + (Subscript[γ, 1])w)

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]]]*);

{(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])}={1.7376726950364059`,0.026017164070679045`,-1.7116555309657269`};
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])}={0,0,0};
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])}={-1.7376726950364059`,-0.026017164070679045`,1.7116555309657269`};
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,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-> (Subscript[x, 1]),y->(Subscript[y, 1]),z->(Subscript[z, 1]),u->(Subscript[u, 1]),w->(Subscript[w, 1])}//N

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

Eigenvalues[J]/.{x-> (Subscript[x, 3]),y->(Subscript[y, 3]),z->(Subscript[z, 3]),u->(Subscript[u, 3]),w->(Subscript[w, 3])}//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,99},Method->Automatic];phaseplot=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solution],{t,20,99},PlotPoints->2000,PlotRange->All,PlotStyle->{RGBColor[1,0,0]},ImageSize->{400,400},BoxRatios->{1,1,1},AxesLabel->{"X","Y","Z"}]

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

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

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

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

plotU=Plot[Evaluate[u[t]/.solution],{t,0,99},PlotRange->{-0.1,0.1},PlotStyle->{Cyan,Thickness[0.005]}]

plotW=Plot[Evaluate[w[t]/.solution],{t,0,99},PlotRange->{(Subscript[z, min]),(Subscript[z, max])},PlotStyle->{Pink,Thickness[0.005]}]

Show[plotX,plotY,plotZ,plotU,plotW]

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

Since the functional Jacobian matrix is locally stationary the nonstationary part of the flow curvature manifold may be removed 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,(Subscript[x, min]),-1},{y,(Subscript[y, min]),(Subscript[y, max])},{z,0,(Subscript[z, max])},PlotPoints->{5,5,5},ColorFunction->Function[{x,y,z,f},ColorData["Rainbow"][z]],Mesh->None,BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"},BoxRatios->{1,1,1}]

flowmanifoldR=ContourPlot3D[ϕ==0,{x,1,(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},{z,(Subscript[z, min]),0},PlotPoints->{5,5,5},ColorFunction->Function[{x,y,z,f},ColorData["Rainbow"][z]],Mesh->None,BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"},BoxRatios->{1,1,1}]

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