Chua's piecewise linear model 4D (Ch 17, page 240)

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

<<VectorFieldPlots`;

System

Piecewise linear function

F[x_ ,a_,b_]=Piecewise[{{b x +b- a,x<-1},{a x,-1≤x≤1},{b x + a-b,x>1}}]

Plot[F[x,-0.42,1.2],{x,-2,2}]

Parameters

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

f[x_,y_,z_,u_]:=(Subscript[α, 1])(z-F[x,a,b])

g[x_,y_,z_,u_]:=-(Subscript[α, 2])y-z-u

h[x_,y_,z_,u_]:=(Subscript[β, 1])( y-x- z)

i[x_,y_,z_,u_ ]:=(Subscript[β, 2])y

Vector field

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

Fixed Points (left side)

f[x_,y_,z_,u_]:=(z-(-a+b+b x)) (Subscript[α, 1])

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

(Subscript[x, 1])=.

(Subscript[y, 1])=.

(Subscript[z, 1])=.

(Subscript[u, 1])=.

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

{(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])}}]];

Jacobian (left side)

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 (left side)

Eigenvalues[J]/.{x-> (Subscript[x, 1]),y->(Subscript[y, 1]),z->(Subscript[z, 1]),(Subscript[α, 1])->2.1429,(Subscript[α, 2])->-0.18,(Subscript[β, 1])->0.0774,(Subscript[β, 2])->0.003,a->-0.42,b->1.2}//N

Phase Portrait

Parameters

(Subscript[α, 1])=2.1429;
(Subscript[α, 2])=-0.18;
(Subscript[β, 1])=0.0774;
(Subscript[β, 2])=0.003;

f[x_,y_,z_,u_]:=(Subscript[α, 1])(z-F[x,-0.42,1.2])

g[x_,y_,z_,u_]:=-(Subscript[α, 2])y-z-u

h[x_,y_,z_,u_]:=(Subscript[β, 1])( y-x- z)

i[x_,y_,z_,u_ ]:=(Subscript[β, 2])y

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,600,1000},PlotPoints->5000,PlotRange->All,PlotStyle->{RGBColor[1,0,0]},ImageSize->{400,400},BoxRatios->{1,1,1},AxesLabel->{"X","Y","Z"}]

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

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

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

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

plotU=Plot[Evaluate[u[t]/.solution],{t,100,600},PlotRange->{(Subscript[z, min]),(Subscript[z, max])},PlotStyle->{Cyan,Thickness[0.005]}]

Show[plotX,plotY,plotZ,plotU]

Flow Curvature Manifold (left side)

(Subscript[α, 1])=.
(Subscript[α, 2])=.
(Subscript[β, 1])=.
(Subscript[β, 2])=.

f[x_,y_,z_,u_]:=(z-(-a+b+b x)) (Subscript[α, 1])

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

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

V2=Dt[V1]/.{Dt[x]->f[x,y,z,u],Dt[y]->g[x,y,z,u], Dt[z]->h[x,y,z,u],Dt[u]->i[x,y,z,u],Dt[(Subscript[α, 1])]->0,Dt[(Subscript[α, 2])]->0,Dt[(Subscript[β, 1])]->0,Dt[(Subscript[β, 2])]->0,Dt[a]->0,Dt[b]->0}

V3=Dt[V2]/.{Dt[x]->f[x,y,z,u],Dt[y]->g[x,y,z,u], Dt[z]->h[x,y,z,u],Dt[u]->i[x,y,z,u],Dt[(Subscript[α, 1])]->0,Dt[(Subscript[α, 2])]->0,Dt[(Subscript[β, 1])]->0,Dt[(Subscript[β, 2])]->0,Dt[a]->0,Dt[b]->0}

V4=Dt[V3]/.{Dt[x]->f[x,y,z,u],Dt[y]->g[x,y,z,u], Dt[z]->h[x,y,z,u],Dt[u]->i[x,y,z,u],Dt[(Subscript[α, 1])]->0,Dt[(Subscript[α, 2])]->0,Dt[(Subscript[β, 1])]->0,Dt[(Subscript[β, 2])]->0,Dt[a]->0,Dt[b]->0}

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

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

Osculating plane (left side)

X={x+((-a+b)/(1+b)),y,z+((a-b)/(1+b)),u+((-a+b)/(1+b))}

Π=Factor[Det[{X,V1,V2,V3}]]/.{u->0}

FullSimplify[(ϕ/Π)]

FullSimplify[Det[J]]

This a consequence of Theorem 8.2 and Corollary 8.2 (Eq. 8.8) page 141.

ϕL=Factor[Det[K]/.{u->0,(Subscript[α, 1])->2.1429,(Subscript[α, 2])->-0.18,(Subscript[β, 1])->0.0774,(Subscript[β, 2])->0.003,a->-0.42,b->1.2}]

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

Flow Curvature Invariance (left side)

ϕ=Factor[Det[K]];

Normal vector

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

Invariance

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

Ψ=Factor[V . η];

FullSimplify[(Ψ/ϕ)]

FullSimplify[Tr[J]]

This a consequence of Theorem 8.2 and Corollary 8.2 (Eq. 8.7) page 141. Thus, the invariance of the osculating plane is stated ipso facto.

ΨL=FullSimplify[V . η/.{u->0,(Subscript[α, 1])->2.1429,(Subscript[α, 2])->-0.18,(Subscript[β, 1])->0.0774,(Subscript[β, 2])->0.003,a->-0.42,b->1.2}]

invariantpartL=ContourPlot3D[ΨL==0,{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),4},{z,-1.5,0.7},PlotPoints->{7,7,7},Mesh->None,ContourStyle->Directive[Orange,Opacity[0.5],Specularity[White,30]],BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"},BoxRatios->{1,1,1}]

Vector field

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

Fixed Points (right side)

f[x_,y_,z_,u_]:=(z-(a-b+b x)) (Subscript[α, 1])

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

(Subscript[x, 1])=.

(Subscript[y, 1])=.

(Subscript[z, 1])=.

(Subscript[u, 1])=.

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

{(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]]};
fixedpoint2=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])}}]];

Jacobian (right side)

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 (right side)

Eigenvalues[J]/.{x-> (Subscript[x, 1]),y->(Subscript[y, 1]),z->(Subscript[z, 1]),(Subscript[α, 1])->2.1429,(Subscript[α, 2])->-0.18,(Subscript[β, 1])->0.0774,(Subscript[β, 2])->0.003,a->-0.42,b->1.2}//N

Flow Curvature Manifold (right side)

(Subscript[α, 1])=.
(Subscript[α, 2])=.
(Subscript[β, 1])=.
(Subscript[β, 2])=.

f[x_,y_,z_,u_]:=(z-(a-b+b x)) (Subscript[α, 1])

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

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

V2=Dt[V1]/.{Dt[x]->f[x,y,z,u],Dt[y]->g[x,y,z,u], Dt[z]->h[x,y,z,u],Dt[u]->i[x,y,z,u],Dt[(Subscript[α, 1])]->0,Dt[(Subscript[α, 2])]->0,Dt[(Subscript[β, 1])]->0,Dt[(Subscript[β, 2])]->0,Dt[a]->0,Dt[b]->0}

V3=Dt[V2]/.{Dt[x]->f[x,y,z,u],Dt[y]->g[x,y,z,u], Dt[z]->h[x,y,z,u],Dt[u]->i[x,y,z,u],Dt[(Subscript[α, 1])]->0,Dt[(Subscript[α, 2])]->0,Dt[(Subscript[β, 1])]->0,Dt[(Subscript[β, 2])]->0,Dt[a]->0,Dt[b]->0}

V4=Dt[V3]/.{Dt[x]->f[x,y,z,u],Dt[y]->g[x,y,z,u], Dt[z]->h[x,y,z,u],Dt[u]->i[x,y,z,u],Dt[(Subscript[α, 1])]->0,Dt[(Subscript[α, 2])]->0,Dt[(Subscript[β, 1])]->0,Dt[(Subscript[β, 2])]->0,Dt[a]->0,Dt[b]->0}

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

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

Osculating plane (right side)

X={x+((a-b)/(1+b)),y,z+((-a+b)/(1+b)),u+((a-b)/(1+b))}

Π=Factor[Det[{X,V1,V2,V3}]]/.{u->0}

FullSimplify[(ϕ/Π)]

FullSimplify[Det[J]]

This a consequence of Theorem 8.2 and Corollary 8.2 (Eq. 8.8) page 141.

ϕR=Factor[Det[K]/.{u->0,(Subscript[α, 1])->2.1429,(Subscript[α, 2])->-0.18,(Subscript[β, 1])->0.0774,(Subscript[β, 2])->0.003,a->-0.42,b->1.2}]

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

Flow Curvature Invariance (left side)

ϕ=Factor[Det[K]];

Normal vector

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

Invariance

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

Ψ=Factor[V . η];

FullSimplify[(Ψ/ϕ)]

FullSimplify[Tr[J]]

This a consequence of Theorem 8.2 and Corollary 8.2 (Eq. 8.7) page 141. Thus, the invariance of the osculating plane is stated ipso facto.

ΨR=FullSimplify[V . η/.{u->0,(Subscript[α, 1])->2.1429,(Subscript[α, 2])->-0.18,(Subscript[β, 1])->0.0774,(Subscript[β, 2])->0.003,a->-0.42,b->1.2}]

invariantpartR=ContourPlot3D[ΨR==0,{x,(Subscript[x, min]),(Subscript[x, max])},{y,-4,(Subscript[y, max])},{z,-0.7,1.5},PlotPoints->{7,7,7},Mesh->None,ContourStyle->Directive[Orange,Opacity[0.5],Specularity[White,30]],BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"},BoxRatios->{1,1,1}]

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