Example 10.2: Centre Manifold Approximation 3D (Ch 10, page 166)

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

<<VectorFieldPlots`;

System

Parameters

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

f[x_,y_,z_]:=- y + x z

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

h[x_,y_,z_]:=- z - ((x^2) +  (y^2)) +  (z^2)

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

fp1=Sort[fp[[1]]];fp2=Sort[fp[[2]]];

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

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

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

Show[phaseplot,fixedpoint1,fixedpoint2,BoxRatios->{1,1,1}]

plotX=Plot[Evaluate[x[t]/.solution],{t,550,600},PlotRange->All,PlotStyle->{RGBColor[0,0,1]}]

plotY=Plot[Evaluate[y[t]/.solution],{t,550,600},PlotRange->All,PlotStyle->{RGBColor[0,1,0]}]

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

Flow Curvature Manifold

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

Centre Manifold Approximation by Flow Curvature Method (Proposition 10.1)

Coefficient (Subscript[a, 10]) (K10)

K10=Factor[-((D[ϕ,x])/(D[ϕ,z]))/.{z->(Subscript[a, 20]) (x^2) + (Subscript[a, 11]) x y + (Subscript[a, 02]) (y^2)}]

Coefficient (Subscript[a, 01]) (K01)

K01=Factor[-((D[ϕ,y])/(D[ϕ,z]))/.{z->(Subscript[a, 20]) (x^2) + (Subscript[a, 11]) x y + (Subscript[a, 02]) (y^2)}]

Coefficient (Subscript[a, 20])

FullSimplify[FullSimplify[(1/2)D[K10,x]/.{y->0}]/.{x->0}]

Solve[(Subscript[a, 20])==-2-(Subscript[a, 20]),(Subscript[a, 20])]

Coefficient (Subscript[a, 11])

FullSimplify[FullSimplify[D[K10,y]/.{y->0}]/.{x->0}]

Solve[(Subscript[a, 11])==-2 (Subscript[a, 11]),(Subscript[a, 11])]

FullSimplify[FullSimplify[D[K01,x]/.{x->0}]/.{y->0}]

Solve[(Subscript[a, 11])==-2 (Subscript[a, 11]),(Subscript[a, 11])]

Coefficient (Subscript[a, 02])

FullSimplify[(1/2)FullSimplify[D[K01,y]/.{x->0}]/.{y->0}]

Solve[(Subscript[a, 02])==-2-(Subscript[a, 2]),(Subscript[a, 02])]

Centre Manifold Approximation by Flow Curvature Method (Identification, cf. remark page 166)

Let' s plug inside the flow curvature manifold the centre manifold approximation with its unknown coefficients. An identification order by order leads to :

ϕ=Factor[ ((x^2)+(y^2)) ((x^6)+(y^6)+2 (-1+z) z-(y^4) (-1+z) z+(x^4) (3 (y^2)+z-(z^2))-2 (y^2) (1+z (-1+z (2+z)))+(x^2) (3 (y^4)-2 (y^2) (-1+z) z-2 (1+z (-1+z (2+z)))))/.{z->(Subscript[a, 20]) (x^2) + (Subscript[a, 11]) x y + (Subscript[a, 02]) (y^2)}]

Collect[ϕ,{x,y,z}]

Coefficient[ϕ,x^4]/.{y-> 0}

Solve[-2-2 (Subscript[a, 20])==0,(Subscript[a, 20])]

Coefficient[ϕ,y^4]/.{x-> 0}

Solve[-2-2 (Subscript[a, 2])==0,(Subscript[a, 2])]

Coefficient[ϕ,x y^3]

Solve[-2 (Subscript[a, 11])==0,(Subscript[a, 11])]