Maple Text Commands for Chapter 8 -  

POINCARE MAPS AND NONAUTONOMOUS SYSTEMS IN THE PLANE 

>	with(DEtools):with(plots):

 

Warning, the name changecoords has been redefined 

Example 1: Poincare first returns for the ODE  

>	sys:=diff(r(t),t)+r(t)^2: returns:=dsolve({sys,r(0)=1},numeric,range=0..16*Pi): evalf(seq(r[i]=returns(2*i*Pi),i=1..8),5);

 

Example 5: Hamiltonian system with two-degrees of freedom  

>	omega1:=2:omega2:=2:H:=(omega1/2)*(p1^2+q1^2)+(omega2/2)*(p2^2+q2^2): hamilton_eqs(H);

 

Figures 8.5(a) and (b): Projections of the Poincare surface-of-section. 

>	poincare(H,t=-100..100,{[0,0.5,1.5,0.5,0]},stepsize=0.1,iterations=4, scene=[p1=-1.5..1.5,q1=-1.5..1.5,q2=-1.5..1.5],3);

 

>	poincare(H,t=0..30,{[0,0.5,1.5,0.5,0]},stepsize=0.005,iterations=3,scene=[p1,q1]);

 

Example 6: The Henon-Heiles Hamiltonian  

>	H:=(p1^2+q1^2+p2^2+q2^2)/2+q1^2*q2-q2^3/3:

 

Figures 8.6(a) and (b): Three-dimensional and two-dimensional surface-of-section.  

>	poincare(H,t=-100..100,{[0,0.06,0.1,-0.2,-0.2]},stepsize=0.1,iterations=4, scene=[p1=-0.3..0.3,q1=-0.3..0.3,q2=-0.3..0.3],3);

 

>	poincare(H,t=-100..100,{[0,0.06,0.1,-0.2,-0.2]},stepsize=0.1,iterations=4,3);

 

Figure 8.11(b): Phase portrait for the Duffing equation  

>	Gamma:=0.5:omega:=1.25:k:=0.3: DEplot([diff(x(t),t)=y(t),diff(y(t),t)=x(t)-k*y(t)-(x(t))^3+Gamma*cos(omega*t)],[x(t),y(t)],t=0..300,[[x(0)=1,y(0)=0.5]],x=-2..2,y=-2..2,stepsize=0.1,linecolor=blue,thickness=1);

 

Figure 8.11(b): Poincare first returns for the Duffing system. 

>

ff:=dsolve({diff(x(t),t)=y(t),diff(y(t),t)=x(t)-k*y(t)-(x(t))^3+Gamma* cos(omega*t),x(0)=1,y(0)=0.5},{x(t),y(t)}, type=numeric,method=classical,output=procedurelist): pt:=array(0..10000):x1:=array(0..10000):y1:=array(0..10000): imax:=4000: for i from 0 to imax do x1[i]:=eval(x(t),ff(i*2*Pi/omega)): y1[i]:=eval(y(t),ff(i*2*Pi/omega)): end do: pts:=[[x1[n],y1[n]]$n=10..imax]: # Plot the points on the Poincare section # pointplot(pts,style=point,symbol=circle,symbolsize=1,color=blue,axes=BOXED,scaling=CONSTRAINED,font=[TIMES,ROMAN,15]);

 

Figure 8.14: Bifurcation diagram for the Duffing system. 

>

G:=array(0..10000):Y:=array(0..10000): x1:=array(0..10000):y1:=array(0..10000):r:=array(0..10000): # Increase Gamma from 0 to 0.45 # jmax:=1000:k:=0.3:omega:=1.25:step:=0.00045:interval:=jmax*step: x1[0]:=1:y1[0]:=0:r[0]:=1: for j from 0 to jmax do G[j]:=step*j: ff := dsolve({diff(x(t),t)=y(t),diff(y(t),t)=-k*y(t)+x(t)-(x(t))^3+G[j]*cos( omega*t),x(0)=x1[j],y(0)=y1[j]},{x(t),y(t)}, type=numeric, output=procedurelist); x1[j+1]:=eval(x(t),ff(2*Pi/omega)): y1[j+1]:=eval(y(t),ff(2*Pi/omega)): r[j+1]:=sqrt((x1[j+1])^2+(y1[j+1])^2): end do: l :=[[G[n],r[n]] $n=0..jmax]: P1:=plot(l, x=0..interval,y=0..2, style=point,symbol=circle,symbolsize=1,color=blue): # Decrease Gamma from 0.45 to 0 # Gb:=array(0..10000): xb:=array(0..10000):yb:=array(0..10000):rb:=array(0..10000): xb[0]:=x1[jmax+1]:yb[0]:=y1[jmax+1]:rb[0]:=sqrt((xb[0])^2+(yb[0])^2): for j from 0 to jmax do Gb[j]:=interval-step*j: ff := dsolve({diff(x(t),t)=-y(t),diff(y(t),t)=-(-k*y(t)+x(t)-(x(t))^3+Gb[j]* cos(omega*t)),x(0)=xb[j],y(0)=yb[j]},{x(t),y(t)}, type=numeric, output=procedurelist); xb[j+1]:=eval(x(t),ff(-2*Pi/omega)): yb[j+1]:=eval(y(t),ff(-2*Pi/omega)): rb[j+1]:=sqrt((xb[j+1])^2+(yb[j+1])^2): end do: l := [[Gb[n], rb[n]] $n=0..jmax]: P2:=plot(l,x=0..interval,y=0..2, style=point,symbol=circle,symbolsize=1,color=blue): display({P1,P2},labels=['Gamma','r']);

 

End of Chapter 8 Commands