Example 11.2: Chua's system (Ch 11, page 195)

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

<<VectorFieldPlots`;

System

Parameters

f[x_,y_,z_]:=(1/ϵ)(z-α (x^3)-β (x^2)-μ x)

g[x_,y_,z_]:=-z

h[x_,y_,z_]:=-a x +y+b z

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

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

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[a]->0,Dt[b]->0,Dt[α]->0,Dt[β]->0,Dt[μ]->0,Dt[ϵ]->0}

Fixed Points Stability - Eigenvalues

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

Flow Curvature Method → Geometric Singular Perturbation Method

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

ϕ=(-6 a x (z^4) α+18 a (x^4) (z^3) (α^2)-18 a (x^7) (z^2) (α^3)+6 a (x^10) z (α^4)-2 a (z^4) β+24 a (x^3) (z^3) α β-42 a (x^6) (z^2) (α^2) β+20 a (x^9) z (α^3) β+6 a (x^2) (z^3) (β^2)-30 a (x^5) (z^2) α (β^2)+24 a (x^8) z (α^2) (β^2)-6 a (x^4) (z^2) (β^3)+12 a (x^7) z α (β^3)+2 a (x^6) z (β^4)-4 (a^2) (x^3) (z^2) α ϵ+9 a (x^2) y (z^2) α ϵ-6 x (y^2) (z^2) α ϵ+6 a b (x^2) (z^3) α ϵ-6 b x y (z^3) α ϵ-6 x (z^4) α ϵ+14 (a^2) (x^6) z (α^2) ϵ-33 a (x^5) y z (α^2) ϵ+21 (x^4) (y^2) z (α^2) ϵ-18 a b (x^5) (z^2) (α^2) ϵ+21 b (x^4) y (z^2) (α^2) ϵ+21 (x^4) (z^3) (α^2) ϵ-10 (a^2) (x^9) (α^3) ϵ+24 a (x^8) y (α^3) ϵ-15 (x^7) (y^2) (α^3) ϵ+12 a b (x^8) z (α^3) ϵ-15 b (x^7) y z (α^3) ϵ-15 (x^7) (z^2) (α^3) ϵ-(a^2) (x^2) (z^2) β ϵ+2 a x y (z^2) β ϵ-2 (y^2) (z^2) β ϵ+2 a b x (z^3) β ϵ-2 b y (z^3) β ϵ-2 (z^4) β ϵ+17 (a^2) (x^5) z α β ϵ-41 a (x^4) y z α β ϵ+28 (x^3) (y^2) z α β ϵ-23 a b (x^4) (z^2) α β ϵ+28 b (x^3) y (z^2) α β ϵ+28 (x^3) (z^3) α β ϵ-22 (a^2) (x^8) (α^2) β ϵ+54 a (x^7) y (α^2) β ϵ-35 (x^6) (y^2) (α^2) β ϵ+27 a b (x^7) z (α^2) β ϵ-35 b (x^6) y z (α^2) β ϵ-35 (x^6) (z^2) (α^2) β ϵ+4 (a^2) (x^4) z (β^2) ϵ-10 a (x^3) y z (β^2) ϵ+8 (x^2) (y^2) z (β^2) ϵ-6 a b (x^3) (z^2) (β^2) ϵ+8 b (x^2) y (z^2) (β^2) ϵ+8 (x^2) (z^3) (β^2) ϵ-15 (a^2) (x^7) α (β^2) ϵ+38 a (x^6) y α (β^2) ϵ-26 (x^5) (y^2) α (β^2) ϵ+19 a b (x^6) z α (β^2) ϵ-26 b (x^5) y z α (β^2) ϵ-26 (x^5) (z^2) α (β^2) ϵ-3 (a^2) (x^6) (β^3) ϵ+8 a (x^5) y (β^3) ϵ-6 (x^4) (y^2) (β^3) ϵ+4 a b (x^5) z (β^3) ϵ-6 b (x^4) y z (β^3) ϵ-6 (x^4) (z^2) (β^3) ϵ+2 (a^3) (x^5) α (ϵ^2)-7 (a^2) (x^4) y α (ϵ^2)+8 a (x^3) (y^2) α (ϵ^2)-3 (x^2) (y^3) α (ϵ^2)-2 (a^2) b (x^4) z α (ϵ^2)+5 a b (x^3) y z α (ϵ^2)-3 b (x^2) (y^2) z α (ϵ^2)+a (x^3) (z^2) α (ϵ^2)-3 (x^2) y (z^2) α (ϵ^2)-2 (a^2) b (x^7) (α^2) (ϵ^2)+5 a b (x^6) y (α^2) (ϵ^2)-3 b (x^5) (y^2) (α^2) (ϵ^2)+2 a (x^6) z (α^2) (ϵ^2)+2 a (b^2) (x^6) z (α^2) (ϵ^2)-3 (b^2) (x^5) y z (α^2) (ϵ^2)-3 b (x^5) (z^2) (α^2) (ϵ^2)+(a^3) (x^4) β (ϵ^2)-4 (a^2) (x^3) y β (ϵ^2)+5 a (x^2) (y^2) β (ϵ^2)-2 x (y^3) β (ϵ^2)-(a^2) b (x^3) z β (ϵ^2)+3 a b (x^2) y z β (ϵ^2)-2 b x (y^2) z β (ϵ^2)-2 x y (z^2) β (ϵ^2)-3 (a^2) b (x^6) α β (ϵ^2)+8 a b (x^5) y α β (ϵ^2)-5 b (x^4) (y^2) α β (ϵ^2)+4 a (x^5) z α β (ϵ^2)+3 a (b^2) (x^5) z α β (ϵ^2)-5 (b^2) (x^4) y z α β (ϵ^2)-5 b (x^4) (z^2) α β (ϵ^2)-(a^2) b (x^5) (β^2) (ϵ^2)+3 a b (x^4) y (β^2) (ϵ^2)-2 b (x^3) (y^2) (β^2) (ϵ^2)+2 a (x^4) z (β^2) (ϵ^2)+a (b^2) (x^4) z (β^2) (ϵ^2)-2 (b^2) (x^3) y z (β^2) (ϵ^2)-2 b (x^3) (z^2) (β^2) (ϵ^2)-(a^2) (x^5) α (ϵ^3)+2 a (x^4) y α (ϵ^3)-(x^3) (y^2) α (ϵ^3)+a b (x^4) z α (ϵ^3)-b (x^3) y z α (ϵ^3)-(x^3) (z^2) α (ϵ^3)-(a^2) (x^4) β (ϵ^3)+2 a (x^3) y β (ϵ^3)-(x^2) (y^2) β (ϵ^3)+a b (x^3) z β (ϵ^3)-b (x^2) y z β (ϵ^3)-(x^2) (z^2) β (ϵ^3)+18 a (x^2) (z^3) α μ-36 a (x^5) (z^2) (α^2) μ+18 a (x^8) z (α^3) μ+6 a x (z^3) β μ-48 a (x^4) (z^2) α β μ+42 a (x^7) z (α^2) β μ-12 a (x^3) (z^2) (β^2) μ+30 a (x^6) z α (β^2) μ+6 a (x^5) z (β^3) μ-a y (z^2) ϵ μ+10 (a^2) (x^4) z α ϵ μ-24 a (x^3) y z α ϵ μ+18 (x^2) (y^2) z α ϵ μ-14 a b (x^3) (z^2) α ϵ μ+18 b (x^2) y (z^2) α ϵ μ+18 (x^2) (z^3) α ϵ μ-16 (a^2) (x^7) (α^2) ϵ μ+40 a (x^6) y (α^2) ϵ μ-27 (x^5) (y^2) (α^2) ϵ μ+20 a b (x^6) z (α^2) ϵ μ-27 b (x^5) y z (α^2) ϵ μ-27 (x^5) (z^2) (α^2) ϵ μ+3 (a^2) (x^3) z β ϵ μ-7 a (x^2) y z β ϵ μ+8 x (y^2) z β ϵ μ-5 a b (x^2) (z^2) β ϵ μ+8 b x y (z^2) β ϵ μ+8 x (z^3) β ϵ μ-20 (a^2) (x^6) α β ϵ μ+52 a (x^5) y α β ϵ μ-38 (x^4) (y^2) α β ϵ μ+26 a b (x^5) z α β ϵ μ-38 b (x^4) y z α β ϵ μ-38 (x^4) (z^2) α β ϵ μ-5 (a^2) (x^5) (β^2) ϵ μ+14 a (x^4) y (β^2) ϵ μ-12 (x^3) (y^2) (β^2) ϵ μ+7 a b (x^4) z (β^2) ϵ μ-12 b (x^3) y z (β^2) ϵ μ-12 (x^3) (z^2) (β^2) ϵ μ-(a^2) (x^2) y (ϵ^2) μ+2 a x (y^2) (ϵ^2) μ-(y^3) (ϵ^2) μ+a b x y z (ϵ^2) μ-b (y^2) z (ϵ^2) μ-a x (z^2) (ϵ^2) μ-y (z^2) (ϵ^2) μ-2 (a^2) b (x^5) α (ϵ^2) μ+6 a b (x^4) y α (ϵ^2) μ-4 b (x^3) (y^2) α (ϵ^2) μ+4 a (x^4) z α (ϵ^2) μ+2 a (b^2) (x^4) z α (ϵ^2) μ-4 (b^2) (x^3) y z α (ϵ^2) μ-4 b (x^3) (z^2) α (ϵ^2) μ-(a^2) b (x^4) β (ϵ^2) μ+4 a b (x^3) y β (ϵ^2) μ-3 b (x^2) (y^2) β (ϵ^2) μ+4 a (x^3) z β (ϵ^2) μ+a (b^2) (x^3) z β (ϵ^2) μ-3 (b^2) (x^2) y z β (ϵ^2) μ-3 b (x^2) (z^2) β (ϵ^2) μ-(a^2) (x^3) (ϵ^3) μ+2 a (x^2) y (ϵ^3) μ-x (y^2) (ϵ^3) μ+a b (x^2) z (ϵ^3) μ-b x y z (ϵ^3) μ-x (z^2) (ϵ^3) μ-18 a (x^3) (z^2) α (μ^2)+18 a (x^6) z (α^2) (μ^2)-6 a (x^2) (z^2) β (μ^2)+24 a (x^5) z α β (μ^2)+6 a (x^4) z (β^2) (μ^2)+a x y z ϵ (μ^2)+(y^2) z ϵ (μ^2)+b y (z^2) ϵ (μ^2)+(z^3) ϵ (μ^2)-6 (a^2) (x^5) α ϵ (μ^2)+16 a (x^4) y α ϵ (μ^2)-13 (x^3) (y^2) α ϵ (μ^2)+8 a b (x^4) z α ϵ (μ^2)-13 b (x^3) y z α ϵ (μ^2)-13 (x^3) (z^2) α ϵ (μ^2)-2 (a^2) (x^4) β ϵ (μ^2)+6 a (x^3) y β ϵ (μ^2)-7 (x^2) (y^2) β ϵ (μ^2)+3 a b (x^3) z β ϵ (μ^2)-7 b (x^2) y z β ϵ (μ^2)-7 (x^2) (z^2) β ϵ (μ^2)+a b (x^2) y (ϵ^2) (μ^2)-b x (y^2) (ϵ^2) (μ^2)+2 a (x^2) z (ϵ^2) (μ^2)-(b^2) x y z (ϵ^2) (μ^2)-b x (z^2) (ϵ^2) (μ^2)+6 a (x^4) z α (μ^3)+2 a (x^3) z β (μ^3)-x (y^2) ϵ (μ^3)-b x y z ϵ (μ^3)-x (z^2) ϵ (μ^3))

Order (ε^0)

Coefficient (Subscript[a, 100]) (A100)

A100=Factor[-((D[ϕ,x])/(D[ϕ,z]))/.{z->(Subscript[Z, 0])+ϵ (Subscript[Z, 1])}]

Factor[Factor[A100/.{y->0,(Subscript[Z, 0])-> (x^3) α+(x^2) β+x μ}]/.ϵ->0]

The derivative of (Subscript[Z, 0]) with respect to x is thus found again.

Coefficient (Subscript[a, 010]) (A010)

A010=Factor[-((D[ϕ,y])/(D[ϕ,z]))/.{z->(Subscript[Z, 0])+ϵ (Subscript[Z, 1])}]

FullSimplify[A010/.{ϵ->0,y->0,(Subscript[Z, 0])-> (x^3) α+(x^2) β+x μ}]

The derivative of (Subscript[Z, 0]) with respect to y is thus found again.

Order (ε^1)

Coefficient (Subscript[a, 100]) (A100)

A100=Factor[-((D[ϕ,x])/(D[ϕ,z]))/.{z->(Subscript[Z, 0])+ϵ (Subscript[Z, 1])}]

Factor[Factor[D[A100,ϵ]/.{(Subscript[Z, 0])-> (x^3) α+(x^2) β+x μ,(Subscript[Z, 1])-> -((a x-y-b (x^3) α-b (x^2) β-b x μ)/(3 (x^2) α+2 x β+μ))}]/.{ϵ->0}]

Factor[D[-((a x-y-b (x^3) α-b (x^2) β-b x μ)/(3 (x^2) α+2 x β+μ)),x]]

The derivative of (Subscript[Z, 1]) with respect to x is thus found again.

Coefficient (Subscript[a, 010]) (A010)

A010=Factor[-((D[ϕ,y])/(D[ϕ,z]))/.{z->(Subscript[Z, 0])+ϵ (Subscript[Z, 1])}]

Factor[Factor[D[A010,ϵ]/.{(Subscript[Z, 0])-> (x^3) α+(x^2) β+x μ,(Subscript[Z, 1])->-((a x-y-b (x^3) α-b (x^2) β-b x μ)/(3 (x^2) α+2 x β+μ))}]/.{ϵ->0}]

The derivative of (Subscript[Z, 1]) with respect to y is thus found again.

Coefficient (Subscript[a, 001]) (A001)

A001=Factor[-((D[ϕ,ϵ])/(D[ϕ,z]))/.{z->(Subscript[Z, 0])+ϵ (Subscript[Z, 1])}]

Factor[Factor[A001/.{(Subscript[Z, 0])-> (x^3) α+(x^2) β+x μ}]/.{ϵ->0}]

Factor[Solve[-((2 a x-2 y-2 b (x^3) α-2 b (x^2) β-2 b x μ+3 (x^2) α (Subscript[Z, 1])+2 x β (Subscript[Z, 1])+μ (Subscript[Z, 1]))/(3 (x^2) α+2 x β+μ))==(Subscript[Z, 1]),(Subscript[Z, 1])]]

Thus, (Subscript[Z, 1]) is found again.

Order (ε^2)

Coefficient (Subscript[a, 001])

A001=Factor[-((D[ϕ,ϵ])/(D[ϕ,z]))/.{z->(Subscript[Z, 0])+ϵ (Subscript[Z, 1])+(ϵ^2) (Subscript[Z, 2])}]

K1=Factor[(1/2)D[A001,ϵ]/.{(Subscript[Z, 0])-> (x^3) α+(x^2) β+x μ,(Subscript[Z, 1])-> -((a x-y-b (x^3) α-b (x^2) β-b x μ)/(3 (x^2) α+2 x β+μ))}]

Factor[K1/.{ϵ->0}]

Factor[Solve[-(1/(2 (((3 (x^2) α+2 x β+μ))^4)))(-9 (a^2) (x^3) α+27 a (x^2) y α-18 x (y^2) α+27 a b (x^5) (α^2)-36 b (x^4) y (α^2)-27 (x^7) (α^3)-18 (b^2) (x^7) (α^3)+6 a x y β-6 (y^2) β+33 a b (x^4) α β-48 b (x^3) y α β-63 (x^6) (α^2) β-42 (b^2) (x^6) (α^2) β+6 a b (x^3) (β^2)-12 b (x^2) y (β^2)-48 (x^5) α (β^2)-30 (b^2) (x^5) α (β^2)-12 (x^4) (β^3)-6 (b^2) (x^4) (β^3)+3 (a^2) x μ-3 a y μ+24 a b (x^3) α μ-36 b (x^2) y α μ-45 (x^5) (α^2) μ-36 (b^2) (x^5) (α^2) μ+3 a b (x^2) β μ-12 b x y β μ-66 (x^4) α β μ-48 (b^2) (x^4) α β μ-24 (x^3) (β^2) μ-12 (b^2) (x^3) (β^2) μ-3 a b x (μ^2)-21 (x^3) α (μ^2)-18 (b^2) (x^3) α (μ^2)-15 (x^2) β (μ^2)-6 (b^2) (x^2) β (μ^2)-3 x (μ^3)+81 (x^8) (α^4) (Subscript[Z, 2])+216 (x^7) (α^3) β (Subscript[Z, 2])+216 (x^6) (α^2) (β^2) (Subscript[Z, 2])+96 (x^5) α (β^3) (Subscript[Z, 2])+16 (x^4) (β^4) (Subscript[Z, 2])+108 (x^6) (α^3) μ (Subscript[Z, 2])+216 (x^5) (α^2) β μ (Subscript[Z, 2])+144 (x^4) α (β^2) μ (Subscript[Z, 2])+32 (x^3) (β^3) μ (Subscript[Z, 2])+54 (x^4) (α^2) (μ^2) (Subscript[Z, 2])+72 (x^3) α β (μ^2) (Subscript[Z, 2])+24 (x^2) (β^2) (μ^2) (Subscript[Z, 2])+12 (x^2) α (μ^3) (Subscript[Z, 2])+8 x β (μ^3) (Subscript[Z, 2])+(μ^4) (Subscript[Z, 2]))==(Subscript[Z, 2]),(Subscript[Z, 2])]]

Thus, (Subscript[Z, 2]) is found again. Let's check.

FullSimplify[((1/(((3 (x^2) α+2 x β+μ))^4))(3 (a^2) (x^3) α-9 a (x^2) y α+6 x (y^2) α-9 a b (x^5) (α^2)+12 b (x^4) y (α^2)+9 (x^7) (α^3)+6 (b^2) (x^7) (α^3)-2 a x y β+2 (y^2) β-11 a b (x^4) α β+16 b (x^3) y α β+21 (x^6) (α^2) β+14 (b^2) (x^6) (α^2) β-2 a b (x^3) (β^2)+4 b (x^2) y (β^2)+16 (x^5) α (β^2)+10 (b^2) (x^5) α (β^2)+4 (x^4) (β^3)+2 (b^2) (x^4) (β^3)-(a^2) x μ+a y μ-8 a b (x^3) α μ+12 b (x^2) y α μ+15 (x^5) (α^2) μ+12 (b^2) (x^5) (α^2) μ-a b (x^2) β μ+4 b x y β μ+22 (x^4) α β μ+16 (b^2) (x^4) α β μ+8 (x^3) (β^2) μ+4 (b^2) (x^3) (β^2) μ+a b x (μ^2)+7 (x^3) α (μ^2)+6 (b^2) (x^3) α (μ^2)+5 (x^2) β (μ^2)+2 (b^2) (x^2) β (μ^2)+x (μ^3)))-((1/(((3 (x^2) α+2 x β+μ))^4))(3 (a^2) (x^3) α-9 a (x^2) y α+6 x (y^2) α-9 a b (x^5) (α^2)+12 b (x^4) y (α^2)+9 (x^7) (α^3)+6 (b^2) (x^7) (α^3)-2 a x y β+2 (y^2) β-11 a b (x^4) α β+16 b (x^3) y α β+21 (x^6) (α^2) β+14 (b^2) (x^6) (α^2) β-2 a b (x^3) (β^2)+4 b (x^2) y (β^2)+16 (x^5) α (β^2)+10 (b^2) (x^5) α (β^2)+4 (x^4) (β^3)+2 (b^2) (x^4) (β^3)-(a^2) x μ+a y μ-8 a b (x^3) α μ+12 b (x^2) y α μ+15 (x^5) (α^2) μ+12 (b^2) (x^5) (α^2) μ-a b (x^2) β μ+4 b x y β μ+22 (x^4) α β μ+16 (b^2) (x^4) α β μ+8 (x^3) (β^2) μ+4 (b^2) (x^3) (β^2) μ+a b x (μ^2)+7 (x^3) α (μ^2)+6 (b^2) (x^3) α (μ^2)+5 (x^2) β (μ^2)+2 (b^2) (x^2) β (μ^2)+x (μ^3)))]

Thus, (Subscript[Z, 2]) is found again provided that the whole expression of the (2^nd) curvature (torsion) is used. Sometimes for sake of simplicity in the computation only the stationary part of the torsion may have been considered.