Example 2.10: Unforced Duffing oscillator (Ch 8, page 144)

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

<<VectorFieldPlots`;

System

Parameters

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

f[x_,y_]:=y

g[x_,y_]:=x-(x^3)

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]]];fp2=Sort[fp[[2]]];fp3=Sort[fp[[3]]];

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

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

{(Subscript[x, 3]),(Subscript[y, 3])}={Extract[fp3,{1}][[2]],Extract[fp2,{2}][[2]]};
fixedpoint3=Show[Graphics[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 3]),(Subscript[y, 3])}]},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]}},DisplayFunction->Identity];

Show[nullclines,fixedpoint1,fixedpoint2,fixedpoint3]

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

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

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

Phase Portrait

Portrait de phase

x=.
y=.
t=.

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];phaseplotD=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,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];

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

x=.
y=.
t=.

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];phaseplotG=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,0,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];

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

x=.
y=.
t=.

solution=NDSolve[{x'[t]==f[x[t],y[t]],x[0]==(Sqrt[1.5]),

y'[t]==g[x[t],y[t]],y[0]==(Sqrt[1.5])},{x,y},

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

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

Show[phaseplot,flowG,flowD,phaseplot1,flow1,fixedpoint1,fixedpoint2,fixedpoint3,Frame->True,ImageSize->{600,600}]

Flow Curvature Manifold

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

Osculating Plane

X={x,y}

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

FullSimplify[(ϕ/Π)/.{x->ϵ,y->ϵ}]

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

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 . η]