Lorenz model (Ch 16, page 229)

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

<<VectorFieldPlots`;

System

Parameters

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

g[x_,y_,z_]:=-(1+σ)y+((x+σ y)/(1+σ))(z -ϵ)

h[x_,y_,z_]:=-β z + ((x+σ y)/(1+σ))(1+σ)(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])}={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])}}]];

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)

Γ=Factor[(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}]

Fixed Points Stability - Eigenvalues

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

Flow Curvature Manifold

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

Centre Manifold Approximation by Flow Curvature Method (Proposition 10.2)

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

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

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

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

Coefficient (Subscript[b, 10]) (B10)

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

Coefficient (Subscript[b, 01]) (B01)

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

Coefficient (Subscript[a, 20])

Factor[Factor[(1/2)D[A10,x]/.{ϵ->0}]/.{x->0}]

Solve[(Subscript[a, 20])==-(((Subscript[a, 20]) (-5+3 β (Subscript[b, 20])))/(2 (-1+β (Subscript[b, 20])))),(Subscript[a, 20])]

Coefficient (Subscript[a, 02])

Factor[Factor[(1/2)D[A01,ϵ]/.{x->0}]/.{ϵ->0}]

Solve[(Subscript[a, 02])==((σ (Subscript[a, 2]))/(2 (-σ+(Subscript[b, 2])+2 σ (Subscript[b, 2])+(σ^2) (Subscript[b, 2])))),(Subscript[a, 02])]/.{(Subscript[a, 20])->0}

Coefficient (Subscript[b, 20])

Factor[Factor[(1/2)D[B10,x]/.{ϵ->0}]/.{x->0}]

Solve[(Subscript[b, 20])==-((-5+3 β (Subscript[b, 20]))/(2 β)),(Subscript[b, 20])]

Coefficient (Subscript[b, 02])

Factor[Factor[(1/2)D[B01,ϵ]/.{x->0}]/.{ϵ->0}]

Solve[(Subscript[b, 02])==((σ (Subscript[b, 2]))/(-σ+2 (Subscript[b, 2])+4 σ (Subscript[b, 2])+2 (σ^2) (Subscript[b, 2]))),(Subscript[b, 02])]

Taking into account the preliminary results, the coefficients A10, A01, B10 and B01 may be refined in order to simplify the computation.

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

A10=Factor[-((D[ϕ,x])/(D[ϕ,y]))/.{y-> (Subscript[a, 11]) ϵ x ,z-> (Subscript[b, 11]) ϵ x + (1/β) (x^2)}];

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

A01=Factor[-((D[ϕ,ϵ])/(D[ϕ,y]))/.{y-> (Subscript[a, 11]) ϵ x ,z-> (Subscript[b, 11]) ϵ x + (1/β) (x^2)}];

Coefficient (Subscript[b, 10]) (B10)

B10=Factor[-((D[ϕ,x])/(D[ϕ,z]))/.{y-> (Subscript[a, 11]) ϵ x ,z-> (Subscript[b, 11]) ϵ x + (1/β) (x^2)}];

Coefficient (Subscript[b, 01]) (B01)

B01=Factor[-((D[ϕ,ϵ])/(D[ϕ,z]))/.{y-> (Subscript[a, 11]) ϵ x ,z-> (Subscript[b, 11]) ϵ x + (1/β) (x^2)}];

Coefficient (Subscript[a, 11])

Factor[Factor[D[A10,ϵ]/.{x->0}]/.{ϵ->0}]

Solve[(Subscript[a, 11])==-((2+(Subscript[a, 11])+2 σ (Subscript[a, 11])+(σ^2) (Subscript[a, 11]))/(((1+σ))^2)),(Subscript[a, 11])]

Coefficient (Subscript[b, 11])

Factor[Factor[D[B10,ϵ]/.{x->0}]/.{ϵ->0}]

Solve[(Subscript[b, 11])==-(((2+(Subscript[a, 11])+2 σ (Subscript[a, 11])+(σ^2) (Subscript[a, 11])) (Subscript[b, 11]))/(1+(Subscript[a, 11])+2 σ (Subscript[a, 11])+(σ^2) (Subscript[a, 11]))),(Subscript[b, 11])]

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

Taking into account the preliminary results, i.e., that the centre manifold is homogeneous polynomial of degree three, let' s plug inside the flow curvature manifold the centre manifold approximation with its unknown coefficients. An identification order by order leads to :

K=Expand[ϕ/.{y->((Subscript[a, 02]) (ϵ^2) + (Subscript[a, 11]) ϵ x + (Subscript[a, 20]) (x^2)),z->((Subscript[b, 02]) (ϵ^2) + (Subscript[b, 11]) ϵ x + (Subscript[b, 20]) (x^2))}];

FullSimplify[Coefficient[K,x^5ϵ^1]]

FullSimplify[Coefficient[K,x^4ϵ^2]]

FullSimplify[Coefficient[K,x^3ϵ^3]]

FullSimplify[Coefficient[K,x^2ϵ^4]]

FullSimplify[Coefficient[K,x^1ϵ^5]]

Solve[{(Subscript[a, 20]) (-1+β (Subscript[b, 20]))==0, (-1+β (((1+σ))^2) (Subscript[a, 20]) (Subscript[b, 11])+β (Subscript[b, 20])+(((1+σ))^2) (Subscript[a, 11]) (-1+β (Subscript[b, 20])))==0, (β ((((1+σ))^2) (Subscript[a, 20]) (Subscript[b, 2])+(1+(((1+σ))^2) (Subscript[a, 11])) (Subscript[b, 11]))+(((1+σ))^2) (Subscript[a, 2]) (-1+β (Subscript[b, 20])))==0,((1+(((1+σ))^2) (Subscript[a, 11])) (Subscript[b, 2])+(((1+σ))^2) (Subscript[a, 2]) (Subscript[b, 11]))==0},{(Subscript[a, 20]) ,(Subscript[a, 11]) ,(Subscript[a, 02]) ,(Subscript[b, 20]),(Subscript[b, 11]),(Subscript[b, 02])}]

There is no difference with the expected coefficients. This is due to the fact that the centre manifold is in this case of degree two. Another approach such as Euler's theorem may be also used.

Centre Manifold Approximation by Euler's theorem (cf. remark page 174)

Using Euler's theorem leads fo the following system:

x((Overscript[y, •])/(Overscript[x, •]))+μ((∂(Subscript[h, 1]))/(∂ϵ))=2(Subscript[h, 1])(x,ϵ)     Eq. (1)

x((Overscript[z, •])/(Overscript[x, •]))+μ((∂(Subscript[h, 2]))/(∂ϵ))=2(Subscript[h, 2])(x,ϵ)     Eq. (2)

Equation (1)

D[((Subscript[a, 02]) (ϵ^2) + (Subscript[a, 11]) ϵ x + (Subscript[a, 20]) (x^2)),ϵ]

K1=Expand[x g[x,y,z] +f[x,y,z]( ϵ (2 ϵ (Subscript[a, 2])+x (Subscript[a, 11]))-2((Subscript[a, 20]) (x^2) +(Subscript[a, 11]) x ϵ + (Subscript[a, 02]) (ϵ^2))) /.{y-> ((Subscript[a, 20]) (x^2) +(Subscript[a, 11]) x ϵ + (Subscript[a, 02]) (ϵ^2)),z-> ((Subscript[b, 20]) (x^2) +(Subscript[b, 11]) x ϵ + (Subscript[b, 02]) (ϵ^2))}]

Collect[K1,{x}]

Coefficient[K1,{x^3}]/.{ϵ->0}

Coefficient[K1,{x^2 ϵ^1}]

Coefficient[K1,{x^1 ϵ^2}]

Coefficient[K1,{ϵ^3}]/.{x->0}

Equation (2)

D[((Subscript[b, 20]) (x^2) +(Subscript[b, 11]) x ϵ + (Subscript[b, 02]) (ϵ^2)),ϵ]

K2=Expand[x h[x,y,z] +f[x,y,z]( ϵ (2 ϵ (Subscript[b, 2])+x (Subscript[b, 11]))-2((Subscript[b, 20]) (x^2) +(Subscript[b, 11]) x ϵ + (Subscript[b, 02]) (ϵ^2))) /.{y-> ((Subscript[a, 20]) (x^2) +(Subscript[a, 11]) x ϵ + (Subscript[a, 02]) (ϵ^2)),z-> ((Subscript[b, 20]) (x^2) +(Subscript[b, 11]) x ϵ + (Subscript[b, 02]) (ϵ^2))}]

Collect[K2,{x}]

Coefficient[K2,{x^3}]/.{ϵ->0}

Coefficient[K2,{x^2 ϵ^1}]

Coefficient[K2,{x^1 ϵ^2}]

Coefficient[K2,{ϵ^3}]/.{x->0}

Solve[{-(Subscript[a, 20])-σ (Subscript[a, 20])==0,-(1/(1+σ))-(Subscript[a, 11])-σ (Subscript[a, 11])==0,-(Subscript[a, 2])-σ (Subscript[a, 2])==0,1-β (Subscript[b, 20])==0,-β (Subscript[b, 11])==0,-β (Subscript[b, 2])==0},{(Subscript[a, 20]) ,(Subscript[a, 11]) ,(Subscript[a, 02]) ,(Subscript[b, 20]),(Subscript[b, 11]),(Subscript[b, 02])}]