restart; Poincare Map of Forced Differential Equation based on worksheet by Richard Enns and George McGuire
<Text-field style="Heading 1" layout="Heading 1">Introduction</Text-field> Topic: In this file, we present a method of generating Poincare sections with large numbers of points. This is particularly useful in looking at the shapes and internal structures of strange attractors, although it can also be applied to finding large period solutions. Warning: This file take about 1 hour to run on a 100 Mhz 16 Meg RAM computer. Reference: 1. Nonlinear Physics with Maple for Scientists and Engineers.
<Text-field style="Heading 1" layout="Heading 1">Calculation of kk separate graphs, each with numpts points</Text-field> with(plots): unprotect(gamma): gamma:=0.25;F:=0.40;alpha:=-1;beta:=1; In the line below kk represents the number of 50 point graphs in the total plot. Reduce the value of kk from 200 to some smaller number, if you wish to reduce the amount of time it takes the file to run. ic:=x(0)=0.09,y(0)=0,phi(0)=1; #the initial conditions vars:=[x(t),y(t),phi(t)]; sys:=diff(x(t),t)=y(t), diff(y(t),t)=-2*gamma*y(t)-alpha*x(t)-beta*x(t)^3+F*cos(phi(t)), diff(phi(t),t)=1; sol:=dsolve({sys,ic}, vars,numeric,output=listprocedure); assign(sol): x:=x(t): y:=y(t): phi:=phi(t): numpts:=50;# number of points in one plot kk:=200; # number of plots, numpts x kk = total number of points , 10,000 for k from 0 to kk-1 do pts:=seq([x(2*Pi*(k*numpts+i)),y(2*Pi*(k*numpts+i))],i=1..numpts); plt(k):=plot([pts],view=[-1.5..1.5,-.5..1],color=black): od: xx:=`xx`:yy:=`yy`:zz:=`zz`: xx:=x(2*Pi*numpts):yy:=y(2*Pi*numpts):zz:=phi(2*Pi*numpts): unassign(x,y,phi); ic:=x(0)=xx,y(0)=yy, phi(0)=zz;
<Text-field style="Heading 1" layout="Heading 1">Plot of the Poincare section</Text-field> display([seq(plt(j),j=0..kk-1)],symbol=POINT,style=point,view=[-1..1.6,-.4..0.8]); Modify the view values in the following command line if you wish to examine specific areas of the plot. Changing the view values allows you to zoom in on smaller and smaller regions to see the fractal-like structures. display([seq(plt(j),j=1..kk-1)],symbol=POINT,style=point,view=[-1..0,-.15..0.4]);