Rikitake model (Ch 14, page 223)

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_]:=α  -  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]]];fp2=Sort[fp[[2]]];fp3=Sort[fp[[3]]];fp4=Sort[fp[[4]]];

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

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

{(Subscript[x, 4]),(Subscript[y, 4]),(Subscript[z, 4])}={Extract[fp4,{1}][[2]],Extract[fp4,{2}][[2]],Extract[fp4,{3}][[2]]};
fixedpoint4=Show[Graphics3D[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 4]),(Subscript[y, 4]),(Subscript[z, 4])}]},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],Dt[α]->0,Dt[β]->0, Dt[μ]->0}]

Flow Curvature Manifold

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

Case 1. Invariant Manifold: μ=α=0

ϕ=FullSimplify[Γ  . Cross[V,γ]/.{μ->0,α->0}]

ϕ=((y^2)+(((z-β))^2)) ((x^2)+(((z+β))^2))

Invariance

Normal vector

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

Invariance

Ψ=FullSimplify[(V . η)/.{μ->0,α->0}]

FullSimplify[(Ψ/ϕ)]

Thus, the nonlinear manifolds  ((y^2)+(((z-β))^2)) and ((x^2)+(((z+β))^2)) are first integral.

Case 2. Invariant Manifold:  μ=β=0

ϕ=FullSimplify[Γ  . Cross[V,γ]/.{μ->0,β->0}]

ϕ=((x^2)-(y^2))

Invariance

Normal vector

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

Invariance

Ψ=FullSimplify[(V . η)/.{μ->0,β->0}]

FullSimplify[(Ψ/ϕ)]

Thus, the nonlinear manifold  (x^2)-(y^2) is first integral.

Case 3. Invariant Manifold:  β=0

ϕ=FullSimplify[Γ  . Cross[V,γ]/.{β->0}]

ϕ=((x^2)-(y^2))

Invariance

Normal vector

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

Invariance

Ψ=FullSimplify[(V . η)/.{β->0}]

FullSimplify[(Ψ/ϕ)]

Thus, the nonlinear manifold  (x^2)-(y^2) is first integral.