restart;
with(DEtools): with(plots):
sigma:=0.25: omega:=1:
Duffing:=diff(u(t),t)=v(t),
diff(v(t),t)=u(t)-u(t)^3-sigma*v(t)+Gamma*cos(omega*t):# Evaluating the Poincare return section at Gamma=0.5.
# There is only one point in the Poincare cross section suggests the existence of a periodic one orbit in the phase space as an attractor.
# I start i from i=500 to remove the transient behavior.
solA:=dsolve({eval(Duffing,Gamma=0.5),u(0)=0.5,v(0)=0.5},
{u(t),v(t)},
type=numeric, method=rkf45, maxfun=10^7,
output=Array([seq(2*Pi*i/omega,i=500..1500)])):
pointplot([solA[2][1][..,2],solA[2][1][..,3]],view=[-2..2,-2..2]);JSFH# The corresponding phase portrait of the forced Duffing oscillator for Gamma=0.5.
# There is only a single periodic one orbit in the phase space, agreeing withe the former result.
DEplot({eval(Duffing, Gamma = .5)},
{u(t),v(t)},
50..250,
{[0,0.5,0.5]},
stepsize=0.01, u=-2..2,v=-2..2,linecolor=burgundy, thickness=1
);JSFH# Evaluating the Poincare return section at Gamma=0.3.
# There are inifinitly many points in the Poincare cross section suggests the existence of a "chaotic attractor" in the phase space.
solB:=dsolve({eval(Duffing,Gamma=0.3),u(0)=0.5,v(0)=0.5},
{u(t),v(t)},
type=numeric, method=rkf45, maxfun=10^7,
output=Array([seq(2*Pi*i/omega,i=500..1500)])):
pointplot([solB[2][1][..,2],solB[2][1][..,3]],view=[-2..2,-2..2]);JSFHJSFH# The corresponding phase portrait of the forced Duffing oscillator for Gamma=0.3.
# In this case the trajectoris in the phase space is chaotic.
DEplot({eval(Duffing, Gamma = 0.3)},
{u(t),v(t)},
50..250,
{[0,0.5,0.5]},
stepsize=0.01, u=-2..2,v=-2..2,linecolor=burgundy, thickness=1
);# If Gamma=0, the Duffing oscillator does not explicitly depend on time anymore, the t=2*Pi*i/omega gives no longer the Poincare return map.
# Indeed, when Gamma is non-zero, the Duffing system is three-dimensional in the perspective of an autonomous system as:
Duffing3d:=diff(u(t),t)=v(t),
diff(v(t),t)=u(t)-u(t)^3-sigma*v(t)+Gamma*cos(omega*theta(t)),
diff(theta(t),t)=1:
# Letting theta(t)=2*Pi*i/omega, we reduce the 3d system to a 2d map.
# By contrast, if Gamma=0, the system is 2d itself, to obtain the Poincare return section, we may set u=v, and study the behavor of trajectories crossing the line u=v transversely. # In the case Gamma=0, I would like to plot intersections between trajectories and the Poincare return section, which is a line: u=v.
solC:=dsolv({eval(Duffing, Gamma=0),u(0)=0.5,v(0)=0.5},
{u(t),v(t)},
type=numeric, method=rkf45, maxfun=10^7,
# Here, I have no idea how to define the seq in this situation, could you please help me.
output=Array([seq(???, i=500..1500)])):
pointplot([solC[2][1][..,2],solC[2][1][..,3]],view=[-2..2,-2..2]);