Example 2.9: First Limit cycle (Ch 2, page 43)

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

<<VectorFieldPlots`;

System

Parameters

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

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

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

Vector field

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

Fixed Points & Nullclines

(Subscript[x, 1])=.
(Subscript[y, 1])=.

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,Axes->True,DisplayFunction->$DisplayFunction]

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

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

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

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

tableplot={};
For[t=4,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];];
flow=ListVectorFieldPlot[tableplot ,ScaleFunction->(.05&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

x=.
y=.
t=.

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

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

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

tableplot={};
For[t=1,t≤2,t=t+ 0.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];];
flow1a=ListVectorFieldPlot[tableplot ,ScaleFunction->(.05&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

tableplot={};
For[t=1.9,t≤2,t=t+ 0.075,
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];];
flow1b=ListVectorFieldPlot[tableplot ,ScaleFunction->(.05&),ColorFunction->(Hue[0.66]&),Frame->True,AspectRatio->1]

Show[phaseplot,phaseplot1,fixedpoint1,flow,flow1a,flow1b,nullclines,DisplayFunction->$DisplayFunction]

vectorfield=VectorFieldPlot[{f[x,y],g[x,y]},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},PlotPoints->15,ScaleFunction->(.25&),ColorFunction->(Hue[0.66]&),Frame->True,Axes->True,AxesLabel->{"X","Y"},AspectRatio->1];

Show[phaseplot,fixedpoint1,nullclines,vectorfield,ImageSize->{600,600}]

Flow Curvature Manifold I

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

The flow curvature manifold corresponds to a circle. Let's first choose the radius R = 0.5.

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

flowmanifold=ContourPlot[{ϕ==0},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},ContourStyle->{{Green,Thickness[0.005],Dashing[{.025}]}},AspectRatio->1,Frame->True]

Flow Curvature Manifold II

Then, let's first choose the radius R = 1.5.

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

flowmanifold1=ContourPlot[{ϕ==0},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},ContourStyle->{{Green,Thickness[0.005],Dashing[{.025}]}},AspectRatio->1,Frame->True]

Poincaré Ring

PoincaréRing=RegionPlot[x^2+y^2>0.5&&x^2+y^2<1.5,{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},Mesh->None]

Show[PoincaréRing,phaseplot,phaseplot1,flowmanifold, flowmanifold1,flow,flow1a,flow1b,fixedpoint1,vectorfield,Frame->True,ImageSize->{600,600} ]

Osculating Plane

Let' s notice that the osculating plane passing through the unique fixed point located at the origin directly provides the limit cycle algebraic equation, i.e., the invariant manifold equation.

X={x,y};

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

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

limitcycle=ContourPlot[{ϕ==0},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},ContourStyle->{{Green,Thickness[0.005],Dashing[{.025}]}},AspectRatio->1,Frame->True]

Show[PoincaréRing,phaseplot,phaseplot1,limitcycle,flow,flow1a,flow1b,fixedpoint1,vectorfield,Frame->True,ImageSize->{600,600} ]