A. Dynamical System Study 2D

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

ClearAll["`*"]
ClearAll[Evaluate[Context[]<>"*"]]

System

Parameters

w1=0.2;w3=0.1; w4=0.0079;
f[x_,y_]:=-3+y^2
g[x_,y_]:=x^2-2

Vector field

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

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

Fixed Points & Nullclines

fp=NSolve[{f[x1,y1]==0,g[x1,y1]==0},{x1,y1}]

fixedpoint={};
For[i=1,i≤ Length[fp],i=i+1,
 (Subscript[fp, i])=Sort[fp[[i]]];
{(Subscript[x, i]),(Subscript[y, i])}={Extract[(Subscript[fp, i]),{1}][[2]],Extract[(Subscript[fp, i]),{2}][[2]]};

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

fixedpoint=Append[fixedpoint,(Subscript[fixpoint, i])]
]

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,fixedpoint]

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

(Subscript[x, 1])

Fixed Points Stability - Eigenvalues

SistemEigen={};
For[i=1,i≤ Length[fp],i=i+1,
 (Subscript[fp, i])=Sort[fp[[i]]];
{(Subscript[x, i]),(Subscript[y, i])}={Extract[(Subscript[fp, i]),{1}][[2]],Extract[(Subscript[fp, i]),{2}][[2]]};
(Subscript[eigen, i])=Eigensystem[J]/.{x-> (Subscript[x, i]),y->(Subscript[y, i])};
SistemEigen=Append[SistemEigen,(Subscript[eigen, i])]
]

SistemEigen

Phase Portrait

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

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

{t,0,10},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,100},PlotPoints->300,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]

Show[phaseplot,fixedpoint,nullclines]

plotX=Plot[Evaluate[x[t]/.solution],{t,0,100},PlotRange->{(Subscript[x, min]),(Subscript[x, max])},PlotStyle->{RGBColor[0,0,1]}]

(*x en funcion de t*)

plotY=Plot[Evaluate[y[t]/.solution],{t,0,100},PlotRange->{(Subscript[y, min]),(Subscript[y, max])},PlotStyle->{RGBColor[0,1,0]}]

Show[plotX,plotY]

vectorplot=VectorPlot[{f[x,y],g[x,y]},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},VectorPoints->20,VectorScale->{.03,Automatic,None},ColorFunction->(Hue[0.66]&),Frame->True,Axes->True,AxesLabel->{"X","Y"},AspectRatio->1];

Show[phaseplot,fixedpoint,nullclines,vectorplot,ImageSize->{400,400}]

data=Table[{f[x,y],g[x,y]},{x,(Subscript[x, min]),(Subscript[x, max]),0.15},{y,(Subscript[y, min]),(Subscript[y, max]),0.15}];

K=ListStreamPlot[data]

Show[phaseplot,fixedpoint,nullclines,K,ImageSize->Automatic]

Flow Curvature Manifold

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

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

Flow Curvature Invariance

Normal Vector

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

Invariance

Ψ=FullSimplify[V . η]

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

Show[phaseplot,fixedpoint,flowmanifold,invariantpart,ImageSize->{600,600}]

Flow

tableplot={};
For[t=10,t≤15,t=t+ .05,
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=VectorPlot[tableplot ,VectorPoints->20,VectorScale->{.03,Automatic,None},ColorFunction->(Hue[0.66]&),Frame->True,Axes->True,AxesLabel->{"X","Y"},AspectRatio->1]

Show[phaseplot,fixedpoint,nullclines,flow,ImageSize->{600,600}]

B. Balthazar Van der Pol (1926)

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

System

Parameters

ε=0.05;

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

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

g[x_,y_]:=-x

Vector field

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

Fixed Points & Nullclines

fp=Solve[{f[x1,y1]==0,g[x1,y1]==0},{x1,y1}]

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

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

Show[nullclines,fixedpoint]

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]==1,

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

{t,0,100},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,80,100},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];

Show[phaseplot,fixedpoint1,nullclines]

plotX=Plot[Evaluate[x[t]/.solution],{t,80,100},PlotRange->{(Subscript[x, min]),(Subscript[x, max])},PlotStyle->{RGBColor[0,0,1]}]

plotY=Plot[Evaluate[y[t]/.solution],{t,80,100},PlotRange->{(Subscript[y, min]),(Subscript[y, max])},PlotStyle->{RGBColor[0,1,0]}]

Show[plotX,plotY]

Vector Field

vectorfield=VectorFieldPlot[{f[x,y],g[x,y]},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},PlotPoints->20 ,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

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

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

Flow Curvature Invariance

Normal Vector

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

Invariance

Ψ=Factor[V . η]

invariantpart=ContourPlot[{Ψ==0},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},PlotPoints->50,ContourStyle->{{Magenta,Thickness[0.005],Dashing[{.025}]}},Frame->True,AspectRatio->1]

Show[phaseplot,fixedpoint1,flowmanifold,invariantpart,Frame->True,ImageSize->{600,600} ]

Attractive Region

attractiveL=RegionPlot[Ψ>0, {x, -2.1, -0.6}, {y, -2.1, 2.1},PlotStyle->Directive[Opacity[0.7],Cyan],Mesh->None]

attractiveR=RegionPlot[Ψ>0, {x, 0.6, 2.1}, {y, -2.1, 2.1},PlotStyle->Directive[Opacity[0.7],Cyan],Mesh->None]

Show[phaseplot,fixedpoint1,flowmanifold,nullclines,attractiveL,attractiveR,invariantpart,Frame->True,ImageSize->{600,600}]

Flow

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

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

Show[phaseplot,fixedpoint1,nullclines,flow,flowmanifold,invariantpart, ImageSize->{600,600}]