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],solA[..,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],solB[..,3]],view=[-2..2,-2..2]);JSFH JSFH # 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],solC[..,3]],view=[-2..2,-2..2]);