Example 9.6: Volterra-Lotka predator-prey model 3D (Ch 9, page 156)

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

<<VectorFieldPlots`;

System

Parameters

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

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

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

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

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

Γ=FullSimplify[(J . γ +Dt[J] . V)/.{Dt[x]->f[x,y,z],Dt[y]->g[x,y,z],Dt[z]->h[x,y,z]}]

Fixed Points Stability - Eigenvalues

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

Flow Curvature Manifold

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

Thus, according to Converse Proposition 9.1, this dynamical system admits at least a linear invariant manifold.

Osculating Plane

X={x,y,z}

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

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

Flow Curvature Invariance

Normal vector

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

Invariance

Ψ=FullSimplify[(V . η)]

FullSimplify[(Ψ/ϕ)]

Collect[e x-a y+c (x+y)+(a+d) z+b (x+z),{x,y,z}]

Solve[{(b+c+e) ==0,(-a+c) ==0,(a+b+d) ==0},{a,b,c,d,e}]

ϕ=(x^a)  (y^b) (z^c)(((x +y))^d)(((y - z))^e)/.{a->c,b->-c-e,d->e}

Normal vector

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

Invariance

Ψ=FullSimplify[(V . η)]

Thus, a first integral for this dynamical system may be built according to Darboux Theorem 5 as:

ϕ = (x^c) (y^(-c-e)) (((x+y))^e) (((y-z))^e) (z^c)

Moreover, it may also be checked that

ϕ=x + y - z

is a first integral for this dynamical system.

Normal vector

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

Invariance

Ψ=FullSimplify[(V . η)]