Forced Duffing Oscillator
Authors' Data: Housam Binous & Nasri Zakia
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reference :

www.wesleyan.edu/physics/courses/phys313/proj/NasimDuffing.pdf 

http://www.yorku.ca/marko/ComPhys/Duffing/Duffing.html

Time series and phase-space showing harmonic motion and 
closed phase space

γ=0.1;ϵ=0.25;Γ=1;Ω=1;ω=1;

sol=NDSolve[{D[x[t],{t,2}]+γ  D[x[t],{t,1}]-(ω^2) x[t]+ϵ ((x[t])^3)==Γ Cos[Ω t],x[0]== 1,x'[0]== 0},{x[t]},{t,0,4000},MaxSteps->Infinity]

Plot[x[t]/.sol,{t,100,200},PlotStyle-> RGBColor[1,0,0]]

f=First[D[x[t]/.sol,{t,1}]]

Plot[f,{t,100,200},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 10000]

ParametricPlot[{First[x[t]/.sol],f},{t,100,200},PlotStyle-> RGBColor[1,0,0]]

Poincare section  (This is done by taking one data point from phase space per period of the driving). As expected from the phase space diagram, we get a single point.

tbl=Table[{First[x[t]/.sol],f},{t,1000,4000,2 Pi/Ω}];

ListPlot[tbl,PlotStyle-> RGBColor[1,0,0],PlotRange->{{-3,3},{-3,3}}]

Time series and phase-space showing chotic behavior and strange attractor when forcing is increased

γ=0.1;ϵ=0.25;Γ=2.5;Ω=2;ω=1;

sol=NDSolve[{D[x[t],{t,2}]+γ  D[x[t],{t,1}]-(ω^2) x[t]+ϵ ((x[t])^3)==Γ Cos[Ω t],x[0]== 1,x'[0]== 0},{x[t]},{t,0,4000},MaxSteps->Infinity]

Plot[x[t]/.sol,{t,100,200},PlotStyle-> RGBColor[1,0,0]]

f=First[D[x[t]/.sol,{t,1}]]

Plot[f,{t,100,200},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 10000]

ParametricPlot[{First[x[t]/.sol],f},{t,100,200},PlotStyle-> RGBColor[1,0,0]]

Poincare section  (This is done by taking one data point from phase space per period of the driving). We get fractals.

tbl=Table[{First[x[t]/.sol],f},{t,100,4000,2 Pi/Ω}];

ListPlot[tbl,PlotStyle-> RGBColor[1,0,0],PlotRange-> All]

Interesting phase-space showing limit cycle when we drop the linear term.

γ=1;ϵ=0.25;Γ=2.5;Ω=0.25;

sol=NDSolve[{D[x[t],{t,2}]+γ  D[x[t],{t,1}]+ϵ ((x[t])^3)==Γ Cos[Ω t],x[0]== 1,x'[0]== 0},{x[t]},{t,0,200}]

f=First[D[x[t]/.sol,{t,1}]]

ParametricPlot[{First[x[t]/.sol],f},{t,100,200},PlotStyle-> RGBColor[1,0,0]]