Example 11.1: Van der Pol model (Ch 11, page 190)

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

<<VectorFieldPlots`;

System

Parameters

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

g[x_,y_]:=-x

Vector field

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

Fixed Points

fp=Solve[{f[(Subscript[x, 1]),(Subscript[y, 1])]==0,g[(Subscript[x, 1]),(Subscript[y, 1])]==0},{(Subscript[x, 1]),(Subscript[y, 1])}]

fp1=Sort[fp[[1]]];

{(Subscript[x, 1]),(Subscript[y, 1])}={Extract[fp1,{1}][[2]],Extract[fp1,{2}][[2]]};
fixedpoint1=Show[Graphics[{PointSize[0.025],Hue[1],Point[{(Subscript[x, 1]),(Subscript[y, 1])}]},PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}}]];

Jacobian

a11=Simplify[D[f[x,y],x]]

a12=Simplify[D[f[x,y],y]]

a21=Simplify[D[g[x,y],x]]

a22=Simplify[D[g[x,y],y]]

J={{a11,a12},{a21,a22}}

MatrixForm[J]

Acceleration

γ=J . V

Fixed Points Stability - Eigenvalues

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

Flow Curvature Method → Geometric Singular Perturbation Method

ϕ=FullSimplify[Det[{V,γ}]]

ϕ=6 (x^4)-2 (x^6)+9 x y+3 (x^3) y+9 (y^2)+9 (x^2) ϵ

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

Order (ϵ^0)

Coefficient (Subscript[a, 10])

(Subscript[a, 10])=-((D[ϕ,x])/(D[ϕ,y]))

FullSimplify[(Subscript[a, 10])/.{y->((x^3)/3)-x,ϵ-> 0}]

FullSimplify[DSolve[{ (Subscript[Y, 0])'[x]==-1+(x^2),(Subscript[Y, 0])[0]==0},(Subscript[Y, 0])[x],x]]

D[ (Subscript[Y, 0])[x],x]

(Subscript[Y, 0])[x_]:=((x^3)/3)-x

Order (ϵ^1)

Coefficient (Subscript[a, 10])

(Subscript[a, 01])=-((D[ϕ,ϵ])/(D[ϕ,y]))

FullSimplify[FullSimplify[(Subscript[a, 01])/.{y->(Subscript[Y, 0])[x]+ϵ (Subscript[Y, 1])[x]+(ϵ^2) (Subscript[Y, 2])[x]+(ϵ^3) (Subscript[Y, 3])[x]}]/.{ϵ->0}]

(Subscript[Y, 1])[x_]:=(x/(1-(x^2)))

(Subscript[Y, 1])[x]

Order (ϵ^2)

Coefficient (Subscript[a, 10])

(Subscript[a, 01])=-((D[ϕ,ϵ])/(D[ϕ,y]))/.{y->(Subscript[Y, 0])[x]+ϵ (Subscript[Y, 1])[x]+(ϵ^2) (Subscript[Y, 2])[x]+(ϵ^3) (Subscript[Y, 3])[x] }

FullSimplify[D[(Subscript[a, 01]),ϵ]/.{ϵ->0}]

Solve[2(Subscript[Y, 2])[x]==-((2 x)/(((-1+(x^2)))^3)),(Subscript[Y, 2])[x]]

(Subscript[Y, 2])[x_]:=(x/(((1-(x^2)))^3))

(Subscript[Y, 2])[x]

As expected there is a difference between the terms of second order in epsilon. So, according to the remark page 192 let's compute the Lie derivative of the flow curvature manifold.

Lie derivative of the Flow Curvature Manifold → Geometric Singular Perturbation Method

Normal vector

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

Invariance

Ψ=FullSimplify[V . η]

ϕ=Numerator[Together[Ψ]]

Order (ϵ^2) improved

Coefficient (Subscript[a, 10])

(Subscript[a, 01])=-((D[ϕ,ϵ])/(D[ϕ,y]))/.{y->(Subscript[Y, 0])[x]+ϵ (Subscript[Y, 1])[x]+(ϵ^2) (Subscript[Y, 2])[x]+(ϵ^3) (Subscript[Y, 3])[x] }

FullSimplify[D[(Subscript[a, 01]),ϵ]/.{ϵ->0}]

(Subscript[Y, 2])[x_]=.

Solve[2(Subscript[Y, 2])[x]==((2 (x+(x^3)))/(((-1+(x^2)))^4)),(Subscript[Y, 2])[x]]

(Subscript[Y, 2])[x_]:=((x(1+(x^2)))/(((1-(x^2)))^4))

(Subscript[Y, 2])[x]

Thus, as expected we find again the right second order term.