Poincare Map for the Forced Pendulum

The pendulum with periodic forcing is a simple dynamical system with very complicated behavior.  It provides a good example of the chaotic dynamics resulting from transverse homoclinic points to a hyperbolic periodic orbit.

Here are some functions to integrate the differential equation and compute forward and backward Poincare maps.  The Poincare map is just the time 2 PI map of the flow.

ForcedPendulumOrbit[{q0_,v0_},tfinal_,{eps_,omega_}]:={q[t],v[t]}/.NDSolve[{q'[t]==v[t],v'[t]==-(1+eps*Cos[omega*t])Sin[q[t]],q[0]==q0,v[0]==v0},{q[t],v[t]},{t,0,tfinal}][[1]]

FPPMap[{q0_,v0_},{eps_,omega_}]:=With[{orb=ForcedPendulumOrbit[{q0,v0},2*Pi/omega,{eps,omega}]},
orb/.t->2*Pi/omega]

FPPInverseMap[{q0_,v0_},{eps_,omega_}]:=With[{orb=ForcedPendulumOrbit[{q0,v0},-2*Pi/omega,{eps,omega}]},
orb/.t->-2*Pi/omega]

FPPOrbit[n_,{q0_,v0_},{eps_,omega_}]:=If[n≥0,NestList[FPPMap[#,{eps,omega}]&,{q0,v0},n],NestList[FPPInverseMap[#,{eps,omega}]&,{q0,v0},-n]]

MapToCylinder[{q_,v_}]:= {-Pi+Mod[q+Pi,2*Pi],v}

FPPCOrbit[n_,{q0_,v0_},{eps_,omega_}]:=If[n≥0,Map[MapToCylinder,NestList[FPPMap[#,{eps,omega}]&,{q0,v0},n]],Map[MapToCylinder,NestList[FPPInverseMap[#,{eps,omega}]&,{q0,v0},-n]]]

Here is an interactive orbit plotter for the Poincare map.

MakeForcedPendulumPlot[{q0_,v0_},n_,{eps_,omega_}]:= With[{orb=FPPOrbit[n,{q0,v0},{eps,omega}]},
ListPlot[Map[MapToCylinder,orb],PlotRange->{{-Pi,Pi},{-3,3}}]]

Manipulate[MakeForcedPendulumPlot[p,n,{eps,omega}],{{p,{0,0}},Locator},{{n,50},0,500,1},{{eps,0.1},0,1},{{omega,1.618},0,5}]

Here is a phase portrait showing many orbits of the Poincare map for parameters eps = 0.1, omega = 1.618

orbittable = Table[FPPCOrbit[500,{0,0.1*i},{0.1,1.618}],{i,1,30}];
orbittable2 = Table[FPPCOrbit[500,{0,-0.1*i},{0.1,1.618}],{i,20,30}];
orbittable3 = Table[FPPCOrbit[500,{0.1*i,0},{0.1,1.618}],{i,14,21}];
poincareplot =ListPlot[orbittable~Join~orbittable2~Join~orbittable3,PlotRange->{{-Pi,Pi},{-3,3}},PlotStyle->{PointSize[0.002]}]

To find stable and unstable manifolds we will also need to integrate the variational equations to find the derivative matrix of the Poincare map.

ForcedPendulumVariationalOrbit[{q0_,v0_,dq0_,dv0_},tfinal_,{eps_,omega_}]:={q[t],v[t],dq[t],dv[t]}/.NDSolve[{q'[t]==v[t],v'[t]==-(1+eps*Cos[omega*t])Sin[q[t]],dq'[t]==dv[t],dv'[t]==-(1+eps*Cos[omega*t])(Cos[q[t]])*dq[t],q[0]==q0,v[0]==v0,dq[0]==dq0,dv[0]==dv0},{q[t],v[t],dq[t],dv[t]},{t,0,tfinal}][[1]]

FPPMapDerivative[{q0_,v0_},{eps_,omega_}]:=Module[{q1,v1,a,b,c,d},
{q1,v1,a,c}= ForcedPendulumVariationalOrbit[{q0,v0,1,0},2*Pi/omega,{eps,omega}]/.t->2*Pi/omega;
{q1,v1,b,d}= ForcedPendulumVariationalOrbit[{q0,v0,0,1},2*Pi/omega,{eps,omega}]/.t->2*Pi/omega;
{{q1,v1},{{a,b},{c,d}}}]

Here is the derivative matrix and the eigenvalues and eigenvectors of the fixed point (Pi,0)

{p,monodromymatrix} = FPPMapDerivative[{Pi,0},{0.1,1.618}]
evals =  Eigenvalues[monodromymatrix]
evects = Eigenvectors[monodromymatrix]

To draw the manifolds we will need to have function which find the image of a curve of initial points.

FPPCurveImage[iclist_,{eps_,omega_},n_]:=Module[{curvesegment,curve,i},
curvesegment = iclist;
curve = iclist;
For[i=1,i≤n,i++,
curvesegment = Map[FPPMap[#,{eps,omega}]&,curvesegment];
curve =  curve~Join~curvesegment;];
curve]

FPPInvCurveImage[iclist_,{eps_,omega_},n_]:=Module[{curvesegment,curve,i},
curvesegment = iclist;
curve = iclist;
For[i=1,i≤n,i++,
curvesegment = Map[FPPInverseMap[#,{eps,omega}]&,curvesegment];
curve =  curve~Join~curvesegment;];
curve]

To compute the unstable manifold, first find a point along the unstable eigenvector and compute its first image under the Poincare map.  Both these points are on approximately on the unstable manifold. Divide the segment between these points into lots of other points and compute their forward iterates to approximate the global unstable manifold.   Do the same for the stable manifold using the inverse map.

x0 = {Pi,0}+0.00002*evects[[1]];
x1 = FPPMap[x0,{0.1,1.618}];
n0 = 2000;
ictable1= Table[(1-i/n0)x0+i/n0*x1,{i,0,n0}];
iters = 4;
wu1=FPPCurveImage[ictable1,{0.1,1.618},iters];

x0 = {Pi,0}+0.00002*evects[[2]];
x1 = FPPInverseMap[x0,{0.1,1.618}];
n0 = 2000;
ictable1= Table[(1-i/n0)x0+i/n0*x1,{i,0,n0}];
iters = 4;
ws1=FPPInvCurveImage[ictable1,{0.1,1.618},iters];

We can use symmetry to find more branches of stable and unstable manifolds.  Then combine with the previous Poincare plot.  The stable and unstable manifolds

Trans[{q_,v_},tr_]:= {q+tr,v}
Reflect[{q_,v_}]:= {-q,-v}
wu2=Map[Trans[#,2*Pi]&,Map[Reflect,wu1]];
ws2= Map[Reflect,ws1];
wu3=Map[Trans[#,-2*Pi]&,wu1];

wsplot0=ListLinePlot[{wu3,ws1},PlotRange->All,
PlotStyle->{{Thick,Blue},{Thick,Red}},PlotLabel->"Branches of stable and unstable manifolds with transverse homoclinic points"]

wsuplot = ListLinePlot[{wu2,wu3,ws1,ws2},PlotRange->All,PlotStyle->{{Thick,Blue},{Thick,Blue},{Thick,Red},{Thick,Red}}]

Show[poincareplot,wsuplot,PlotLabel->"Poincare Map for the Forced Pendulum with Stable and Unstable Manifolds"]