Example 10.1: Centre Manifold Approximation 2D (Ch 10, page 163)

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

<<VectorFieldPlots`;

System

Parameters

a=-0.25;

{(Subscript[x, min]),(Subscript[x, max])}={-4,4};
{(Subscript[y, min]),(Subscript[y, max])}={-4,4};

f[x_,y_]:=x y

g[x_,y_]:=-y + a (x^2)

Vector field

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

Fixed Points & Nullclines

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[0.33],Point[{(Subscript[x, 1]),(Subscript[y, 1])}]},PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}}]];

nullclines=ContourPlot[{f[x,y]==0,g[x,y]==0},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},ContourStyle->{{Green,Thickness[0.005]}}];

Show[nullclines,fixedpoint1]

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])}

Phase Portrait

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==0.01,

y'[t]==g[x[t],y[t]],y[0]==6},{x,y},

{t,0,500},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot1=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,500},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True]

tableplot={};
For[t=0,t≤5,t=t+1.25,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow1=ListVectorFieldPlot[tableplot ,ColorFunction->(Hue[0.66]&),ScaleFactor->0.1,MaxArrowLength->5,Frame->True,AspectRatio->1]

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==-0.01,

y'[t]==g[x[t],y[t]],y[0]==6},{x,y},

{t,0,1500},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot2=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,1500},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True]

tableplot={};
For[t=0,t≤5,t=t+1.25,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow2=ListVectorFieldPlot[tableplot ,ColorFunction->(Hue[0.66]&),ScaleFactor->0.1,MaxArrowLength->5,Frame->True,AspectRatio->1]

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==0,

y'[t]==g[x[t],y[t]],y[0]==10},{x,y},

{t,0,50},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot3=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,50},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True]

tableplot={};
For[t=2,t≤4,t=t+3,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow3=ListVectorFieldPlot[tableplot ,ColorFunction->(Hue[0.66]&),ScaleFactor->0.1,MaxArrowLength->5,Frame->True,AspectRatio->1]

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==0,

y'[t]==g[x[t],y[t]],y[0]==-10},{x,y},

{t,0,50},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot4=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,50},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True]

tableplot={};
For[t=2,t≤4,t=t+3,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow4=ListVectorFieldPlot[tableplot ,ColorFunction->(Hue[0.66]&),ScaleFactor->0.1,MaxArrowLength->5,Frame->True,AspectRatio->1]

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==4,

y'[t]==g[x[t],y[t]],y[0]==-2},{x,y},

{t,0,500},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot5=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,500},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True]

tableplot={};
For[t=0,t≤2,t=t+.5,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow5=ListVectorFieldPlot[tableplot ,ColorFunction->(Hue[0.66]&),ScaleFactor->0.1,MaxArrowLength->5,Frame->True,AspectRatio->1]

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==-4,

y'[t]==g[x[t],y[t]],y[0]==-2},{x,y},

{t,0,500},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot6=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,500},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True]

tableplot={};
For[t=0,t≤2,t=t+.5,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{vx,vy}};
tableplot=Append[tableplot,tmptable];];
flow6=ListVectorFieldPlot[tableplot ,ColorFunction->(Hue[0.66]&),ScaleFactor->0.1,MaxArrowLength->5,Frame->True,AspectRatio->1]

Flow Curvature Manifold

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

Centre Manifold Approximation

ϕ=y+0.25 (x^2)

centremanifold=ContourPlot[{ϕ==0},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},PlotPoints->25,ContourStyle->{{Green,Thickness[0.005]}},AspectRatio->1,Frame->True]

Center Manifold Invariance

Normal Vector

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

Invariance

Ψ=FullSimplify[V . η]

K=FullSimplify[(Ψ/ϕ)]

Neglecting order terms of order greater than 2 the centre manifold is invariant according to Darboux Invariance Theorem.

FullSimplify[(((x^2) (-0.25)- y)/(0.25 (x^2)+y))]

Center Manifold Plot

The following "artificial system" has been introduced in order to have an oriented centre manifold.

f[x_,y_]:=1

g[x_,y_]:=2a x

Vitesse

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

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==0,

y'[t]==g[x[t],y[t]],y[0]==0},{x,y},

{t,0,50},Method->Automatic];phaseplotD=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,10},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Green,Thickness[0.005]}},AspectRatio->1,Frame->True]

x=.
y=.
t=.

tableplot={};
For[t=0,t≤10,t=t+1,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{x,y},{-vx,-vy}};
tableplot=Append[tableplot,tmptable];];
flowD=ListVectorFieldPlot[tableplot ,ColorFunction->None,ScaleFactor->0.1,MaxArrowLength->5,Frame->True,AspectRatio->1]

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==0,

y'[t]==g[x[t],y[t]],y[0]==0},{x,y},

{t,0,10},Method->Automatic];phaseplotG=ParametricPlot[Evaluate[{-x[t],y[t]}/.solution],{t,0,10},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Green,Thickness[0.005]}},AspectRatio->1,Frame->True]

x=.
y=.
t=.

tableplot={};
For[t=0,t≤10,t=t+1,
x=Evaluate[{x[t],y[t]}/.solution][[1]][[1]];
y=Evaluate[{x[t],y[t]}/.solution][[1]][[2]];
vx=f[x,y];
vy=g[x,y];
tmptable={{-x,y},{vx,-vy}};
tableplot=Append[tableplot,tmptable];];
flowG=ListVectorFieldPlot[tableplot ,ColorFunction->None,ScaleFactor->0.1,MaxArrowLength->5,Frame->True,AspectRatio->1]

Show[phaseplot1,phaseplot2,phaseplot3,phaseplot4,phaseplot5,phaseplot6,centremanifold,fixedpoint1,flow1,flow2,flow3,flow4,flow5,flow6,flowG,flowD,ImageSize->{600,600}]

Centre Manifold Approximation by Flow Curvature Method (Proposition 10.1)

Parameters

f[x_,y_]:=x y

g[x_,y_]:=-y + α (x^2)

Vector field

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

Flow curvature manifold

ϕ=x (-(y^3)-(x^2) y α-(x^2) (y^2) α+(x^4) (α^2))

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

K10=FullSimplify[-((D[ϕ,x])/(D[ϕ,y]))/.{y-> (Subscript[a, 20])(x^2)+(Subscript[a, 30])(x^3)+(Subscript[a, 40])(x^4)+(Subscript[a, 50])(x^5)}]

Coefficient (Subscript[a, 20]) (K20)

K20=Factor[(1/2)D[K10,x]]

FullSimplify[K20/.{x->0}]

Solve[{(Subscript[a, 2])==(1/2) (5 α-3 (Subscript[a, 20]))},{(Subscript[a, 20])}]

Coefficient (Subscript[a, 30]) (K30)

K30=Factor[(1/6)D[K10,x,x]]

FullSimplify[K30/.{x->0}]

Solve[{(Subscript[a, 30])==-(Subscript[a, 30])},{(Subscript[a, 30])}]

Coefficient (Subscript[a, 40]) (K40)

K40=Factor[(1/24)D[K10,x,x,x]]

FullSimplify[K40/.{x->0}]

-((5 α (Subscript[a, 20]))/2)-3 (Subscript[a, 20]^2)+((2 (Subscript[a, 20]^3))/α)-((3 (Subscript[a, 40]))/4)/.{(Subscript[a, 20])->α}

Solve[{(Subscript[a, 40])==-((7 (α^2))/2)-((3 (Subscript[a, 40]))/4)},{(Subscript[a, 40])}]

Coefficient (Subscript[a, 50]) (K50)

K50=Factor[(1/120)D[K10,x,x,x,x]]

FullSimplify[K50/.{x->0}]

-((2 (5 (α^2)+12 (α-(Subscript[a, 20])) (Subscript[a, 20])) (Subscript[a, 30])+3 α (Subscript[a, 50]))/(5 α))/.{(Subscript[a, 30])->0}

Solve[{(Subscript[a, 5])==-((3 (Subscript[a, 50]))/5)},{(Subscript[a, 50])}]

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

Let' s plug inside the flow curvature manifold the centre manifold approximation with its unknown coefficients. An identification order by order leads to :

ϕ=Factor[x (-(y^3)-(x^2) y α-(x^2) (y^2) α+(x^4) (α^2))/.{y-> (Subscript[a, 20])(x^2)+(Subscript[a, 30])(x^3)+(Subscript[a, 40])(x^4)+(Subscript[a, 50])(x^5)}]

Collect[ϕ,x]

Coefficient[ϕ,x^5]

Solve[(α^2)-α (Subscript[a, 20])==0,(Subscript[a, 20])]

Coefficient[ϕ,x^6]

Solve[-α (Subscript[a, 30])==0,(Subscript[a, 30])]

Coefficient[ϕ,x^7]

-α (Subscript[a, 20]^2)-(Subscript[a, 20]^3)-α (Subscript[a, 40])/.{(Subscript[a, 20])->α}

Solve[-2 (α^3)-α (Subscript[a, 40])==0,(Subscript[a, 40])]

Coefficient[ϕ,x^8]

-2 α (Subscript[a, 20]) (Subscript[a, 30])-3 (Subscript[a, 20]^2) (Subscript[a, 30])-α (Subscript[a, 50])/.{(Subscript[a, 30])->0}

Solve[-α (Subscript[a, 50])==0,(Subscript[a, 50])]