Example 12.4: First integral of quadratic system (Ch 12, page 204)

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

<<VectorFieldPlots`;

System

Parameters

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

f[x_,y_,z_]:=1 + ((x^2)/2) - (y^2)

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

h[x_,y_,z_]:=-x z

Vector field

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

Fixed Points

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

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

FullSimplify[y ((x^2)+2 (((1+y))^2)) (z^2)-38/.{x->0,y->1}]

Solve[-38+8 (z^2)==0]

{(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])}={0,1,((Sqrt[19])/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

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

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]==2,

z'[t]==h[x[t],y[t],z[t]],z[0]==1},{x,y,z},

{t,0,400},Method->Automatic];phaseplot1=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solution],{t,380,400},PlotPoints->2000,PlotStyle->{Red,Thickness[0.007]},PlotRange->All,ImageSize->{400,400},AspectRatio->1,AxesLabel->{"X","Y","Z"}]

Show[phaseplot,fixedpoint1]

Flow Curvature Manifold

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

ϕ=(y^a)(z^b)((((x^2)+2 (((-1+y))^2)))^c)((((x^2)+2 (((1+y))^2)))^d)

Flow Curvature Invariance

Normal vector

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

Invariance

Ψ=Factor[(V . η)]

Let's choose a = 1, b = 2, c = 0 and d = 1

ϕ=((((x^2)+2 (((-1+y))^2)))^0) (y^1) ((((x^2)+2 (((1+y))^2)))^1) (z^2)

FullSimplify[y ((x^2)+2 (((1+y))^2)) (z^2)/.{x->1,y->2,z->1}]

ϕ=y ((x^2)+2 (((1+y))^2)) (z^2)

flowmanifold=ContourPlot3D[ϕ==38,{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},{z,(Subscript[z, min]),(Subscript[z, max])},PlotPoints->{15,15,15},Mesh->None,ContourStyle->Directive[Orange,Opacity[0.5],Specularity[White,30]],BoxRatios->{1,1,1},Axes->True,AxesLabel->{"x","y","z"},AspectRatio->1]

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