The periodically forced damped Duffing oscillator

2-D formulation

Clear[diffX,diffY,x,y,f]
imagesize=345;maxT=30;
ω=1;
δ=0.02;

Manipulate[
diffX=y;diffY=x-(x^3)-δ y;
difX=y[t];difY=x[t]-((x[t])^3)-δ y[t]+f Cos[ω t];
splot=StreamPlot[{diffX,diffY},{x,-3,3},{y,-3,3},FrameLabel->{x[t],y[t]},ImageSize->imagesize,PlotLabel->"Phase portrait, NO forcing present",LabelStyle->15];
fixedpoints=Quiet[Solve[{diffX==0,diffY==0},{x,y}]];
coordinates={x,y}/.fixedpoints;
fixedp=Graphics[{PointSize[0.03],Red,Point[coordinates]}];
integration=NDSolve[{x'[t]==difX,y'[t]==difY,
x[0]==point[[1]],y[0]==point[[2]]},{x,y},{t,0,T}];
trajectory=ParametricPlot[Evaluate[{x[t],y[t]}/.integration],
{t,0,T},PlotStyle->Black];
solx=Plot[Evaluate[x[t]/.integration],
{t,0,T},PlotRange-> All,AxesLabel-> {t,x[t]},PlotLabel->"Displacement",ImageSize->imagesize,LabelStyle->15,PlotStyle->Black];
soly=Plot[Evaluate[y[t]/.integration],
{t,0,T},PlotRange-> All,AxesLabel-> {t,y[t]},PlotLabel->"Velocity",ImageSize->imagesize,LabelStyle->15,PlotStyle->Black];
Grid[{{Show[splot,trajectory,fixedp],soly},{solx}},Frame->All],
{{T,0.5 maxT,"Integration time, t"},1,maxT,Appearance->"Labeled"},{{f,0.84,"External forcing, F"},0,1.4,Appearance->"Labeled"},{{point,{0.1,0.1}},Locator},SaveDefinitions->True]

3-D formulation

tmax=40 π;
x0=0.1;
y0=-1;
ω=0.4;
δ=0.02;
f=0.81;

soln=NDSolve[{
x'[t]==y[t],
y'[t]==x[t]-((x[t])^3)-δ y[t]+f Cos[z[t]],
z'[t]==ω,
x[0]==x0,y[0]==y0,z[0]==ω},
{x,y,z},{t,0,tmax},MaxSteps->Infinity];
diffx=y;diffy=x-(x^3)-δ y;diffz=z;
fixedpoints=Quiet[Solve[{diffx==0,diffy==0,diffz==0},{x,y,z}]];
Print[Style["Coordinates of fixed points",18,Blue]]
coordinates={x,y,z}/.fixedpoints
Gfp=Graphics3D[{PointSize[0.03],Red,Point[coordinates]}];
pplot=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.soln],{t,0,tmax},PlotPoints->1000,AxesLabel->{x[t],y[t],z[t]},BoxRatios->{1, 1, 1.5},LabelStyle->20];Show[pplot,Gfp]
pplott=ParametricPlot3D[Evaluate[{x[t],y[t],t}/.soln],{t,0,tmax},PlotPoints->1000,AxesLabel->{x[t],y[t],t},BoxRatios->{1, 1, 1.5},LabelStyle->20];Show[pplott,Gfp]
px=Plot[Evaluate[x[t]/.soln],{t,0,tmax},AxesLabel-> {t,x[t]},PlotRange-> All,ImageSize->260,LabelStyle->20];
py=Plot[Evaluate[y[t]/.soln],{t,0,tmax},AxesLabel-> {t,y[t]},PlotRange-> All,ImageSize->260,LabelStyle->20];
pz=Plot[Evaluate[z[t]/.soln],{t,0,tmax},AxesLabel-> {t,z[t]},PlotRange-> All,ImageSize->260,LabelStyle->20];
Row[{px,py,pz}]

The Poincaré section

ω=1.12;
δ=0.15;
f=0.3;

x0=0;y0=0;maxt=50000;
data=Block[{δ0=δ,f0=f,ω0=ω},Reap[NDSolve[{
x'[t]== y[t],
y'[t]== x[t]-((x[t])^3)-δ y[t]+f Cos[ω t],
x[0]==x0,y[0]==y0,
WhenEvent[Mod[t,2 π/ω]==0,Sow[{x[t],y[t]}]]},{},{t,0,maxt},MaxSteps->∞]]][[-1,1]];
diffX=y;diffY=x-(x^3)-δ y;
fixedpoints=Quiet[Solve[{diffX==0,diffY==0},{x,y}]];
Print[Style["Coordinates of fixed points, non-forced",18,Blue]]
coordinates={x,y}/.fixedpoints
Gfp=Graphics[{PointSize[0.02],Red,Point[coordinates]}];
lp=ListPlot[data,ImageSize->Large,PlotRange->{{-1.5,1.5},All},PlotStyle->AbsolutePointSize[1],FrameLabel->{x[t],y[t]},PlotLabel->Style["Poincaré section",Large],Frame->True,LabelStyle->20];
Show[lp,Gfp]

Creation of this notebook was supported by IT Academy program of Information Technology Foundation for Education (HITSA).