# Discontinuities in a Poincare map for a double pendulum

I'm generating poincare sections of a double pendulum, and they mostly look okay, but some of them have weird discontinuities that seem wrong. The condition for these sections is the standard $$\theta_1 = 0$$ and $$\dot{\theta}_1 > 0$$. Looking at one of the maps, we see that most of the sections look fine, but there are some weird intersecting sections in the bottom:

If we look at just one of these sections, we see that the top part is being flipped upside down:

Indeed, if we just multiply the bottom bit by −1, it looks just fine:

This is weird, right? None of the analogous plots I've seen in the literature look like this; they all have symmetric limits along the vertical axis.

The fact that the other sections seem fine makes me think that I transcribed the equations correctly, but this problem has been extremely frustrating to diagnose. I've been able to figure out that all of the affected sections are those whose initial condition required $$\dot{\theta}_1<0$$ in order to get a condition with the correct total energy, but I have no idea why that would make them look like this.

Any ideas about what I'm missing?

If you want some minimal code which reproduces the error, this here's a quick C RK4 implementation with the relevant equations of motion,

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* hardwired parameters */

#define PI 3.14159265
#define N 4 /* number of equations to be solved */
#define G 9.8 /* acc'n due to gravity, in m/s^2 */
#define L1 1.0 /* length of pendulum 1 in m */
#define L2 1.0 /* length of pendulum 2 in m */
#define M1 1.0 /* mass of pendulum 1 in kg */
#define M2 1.0 /* mass of pendulum 2 in kg */

void runge_kutta(float xin, float yin[], float yout[], float h);
void derivs(float xin, float yin[], float dydx[]);

int main(int argc, char *argv[])
{

int i = 0, NSTEP;
float h, TMIN, TMAX, TH10, W10, TH20, W20;
float yin[N], yout[N];
float *t, *th1, *th2, *w1, *w2;

/* obtain command line values */

TMIN = atof(argv[1]);
TMAX = atof(argv[2]);
TH10 = atof(argv[3]);
W10 = atof(argv[4]);
TH20 = atof(argv[5]);
W20 = atof(argv[6]);
NSTEP = atoi(argv[7]);

/* allocate memory for arrays of values of time, angles 1 and 2,
and angular velocities 1 and 2 respectively */

t = (float *) malloc(NSTEP*sizeof(float));
th1 = (float *) malloc(NSTEP*sizeof(float));
w1 = (float *) malloc(NSTEP*sizeof(float));
th2 = (float *) malloc(NSTEP*sizeof(float));
w2 = (float *) malloc(NSTEP*sizeof(float));

/* stepsize for integration */

h = (TMAX - TMIN)/(NSTEP - 1.0);

/* Define array of t values */

for (i = 0; i < NSTEP; i++)
t[i] = TMIN + h*i;

/* initial values - convert all angles to radians */

th1[0] = TH10*PI/180.0;
w1[0] = W10*PI/180.0;
th2[0] = TH20*PI/180.0;
w2[0] = W20*PI/180.0;

/* perform the integration */

printf("%f %f %f %f %f\n", t[0], th1[0], w1[0], th2[0], w2[0]);
for (i = 0; i < NSTEP - 1; i++)
{
yin[0] = th1[i];
yin[1] = w1[i];
yin[2] = th2[i];
yin[3] = w2[i];
runge_kutta(t[i], yin, yout, h);
th1[i+1] = yout[0];
w1[i+1] = yout[1];
th2[i+1] = yout[2];
w2[i+1] = yout[3];

printf("%f %f %f %f %f\n", t[i+1], th1[i+1], w1[i+1], th2[i+1],
w2[i+1]);
}

return 0;

}

void derivs(float xin, float yin[], float dydx[])
{

/* function to fill array of derivatives dydx at xin */

float den1, den2, del;

dydx[0] = yin[1];

del = yin[2]-yin[0];

den1 = (M1+M2)*L1 - M2*L1*cos(del)*cos(del);

dydx[1] = (M2*L1*yin[1]*yin[1]*sin(del)*cos(del)
+ M2*G*sin(yin[2])*cos(del) + M2*L2*yin[3]*yin[3]*sin(del)
- (M1+M2)*G*sin(yin[0]))/den1;

dydx[2] = yin[3];

den2 = (L2/L1)*den1;

dydx[3] = (-M2*L2*yin[3]*yin[3]*sin(del)*cos(del)
+ (M1+M2)*G*sin(yin[0])*cos(del)
- (M1+M2)*L1*yin[1]*yin[1]*sin(del)
- (M1+M2)*G*sin(yin[2]))/den2;

return;

}

void runge_kutta(float xin, float yin[], float yout[], float h)
{
/* fourth order Runge-Kutta - see e.g. Numerical Recipes */

int i;
float hh, xh, dydx[N], dydxt[N], yt[N], k1[N], k2[N], k3[N], k4[N];

hh = 0.5*h;
xh = xin + hh;

derivs(xin, yin, dydx); /* first step */
for (i = 0; i < N; i++)
{
k1[i] = h*dydx[i];
yt[i] = yin[i] + 0.5*k1[i];
}

derivs(xh, yt, dydxt); /* second step */
for (i = 0; i < N; i++)
{
k2[i] = h*dydxt[i];
yt[i] = yin[i] + 0.5*k2[i];
}

derivs(xh, yt, dydxt); /* third step */
for (i = 0; i < N; i++)
{
k3[i] = h*dydxt[i];
yt[i] = yin[i] + k3[i];
}

derivs(xin + h, yt, dydxt); /* fourth step */
for (i = 0; i < N; i++)
{
k4[i] = h*dydxt[i];
yout[i] = yin[i] + k1[i]/6. + k2[i]/3. + k3[i]/3. + k4[i]/6.;
}

return;

}


run this with a

$./pendulum.o TMIN TMAX TH10 W10 TH20 W20 NSTEP > pendulum.txt  The initial condition to create the single section above is $$TH_{10} = 0.0, W_{10} = -0.34590943, TH_{20} = 0.02453712, W_{20} = 3.76258371.$$ And here's some Python which will extract the Poincare map: def get_section(data, tolerance): th1 = data[:,1] th1d = data[:,2] th1_good = np.where(np.abs(th1) < tolerance) th1d_good = np.where(th1d > 0) poininds = np.intersect1d(th1_good, th1d_good) poinpoints = np.zeros((len(poininds), 2)) th2 = data[:,3] th2d = data[:,4] poinpoints[:,0] = th2[poininds] poinpoints[:,1] = th2d[poininds] return(poinpoints) data = np.genfromtxt('pendulum.txt') section = get_section(data, 0.0001) plt.scatter(section[:,0], section[:,1])  In response to the comment from @alex-trounev: Weirdly, this problem is there regardless if I use $$\theta_1$$ and $$\dot{\theta}_1$$ or $$\theta_2$$ and $$\dot{\theta}_2$$ for the surface of section in the Poincare map. In the following figure, the $$\theta_1$$ condition is in blue, while the $$\theta_2$$ condition is in orange: • looks like a modularity problem with$\pi $being an upper limit fpr$P_{\theta} $– AoZora Jul 17 at 7:09 • Check the initial data. In your code, the initial data is multiplied by$\pi / 180$. The oscillation amplitude for this data is about 0.01, and in your figure it is about 0.6. – Alex Trounev Jul 17 at 7:56 • Why are using float instead of double? – Kyle Kanos Jul 17 at 11:15 • I think$\pi$being an upper limit is just a coincidence for the section I happened to pull out. Other sections have the same problem higher and lower than what looks like$\pi$. – David Jul 17 at 11:51 • I'm voting to close this question as off-topic because it's about debugging code rather than physics – Kyle Kanos Jul 22 at 11:56 ## 1 Answer I calculated one section with data $$TH10 = 0; W10 = -.34590943; TH20 = .02453712; W20 = 3.76258371;$$ using Mathematica 12. We set for the upper pendulum the deflection angle $$\theta _1$$, and for the lower pendulum $$\theta _2$$. There are two possibilities to choose an event. If we put $$\theta _1=0, \theta _1'>0$$, then the pairs $$\theta _2, \theta _2'$$ form a figure on the plane as in Figure 1 on the left. This is what the author received. But if we put $$\theta _2=0, \theta _2'>0$$, then the pairs $$\theta _1, \theta _1'$$ form a figure on the plane as in Figure 1 on the right (I used to build a figure of only 4000 events. Therefore, separate points are visible.). The code for this problem in Wolfram is the next one (here are the Lagrangian Lgand equations e2,e3):  << VariationalMethods Needs["DifferentialEquationsNDSolveProblems"]; Needs["DifferentialEquationsNDSolveUtilities"]; tm = 10^4; L1 = 1; L2 = 1; m1 = 1; m2 = 1; g = 9.81; TH10 = 0; W10 = -.34590943; TH20 = .02453712; W20 = 3.76258371; x1 = L1*Sin[\[Theta]1[t]]; x2 = L1*Sin[\[Theta]1[t]] + L2*Sin[\[Theta]2[t]]; y1 = -L1*Cos[\[Theta]1[t]]; y2 = -L1*Cos[\[Theta]1[t]] - L2*Cos[\[Theta]2[t]]; v1 = m1*g*y1; v2 = m2*g*y2; V = v1 + v2; t1 = TrigReduce[0.5*m1*((D[x1, {t, 1}])^2 + (D[y1, {t, 1}])^2)]; t2 = TrigReduce[0.5*m2*((D[x2, {t, 1}])^2 + (D[y2, {t, 1}])^2)]; T = t1 + t2; Lg = T - V; e1 = EulerEquations[Lg, {\[Theta]1[t], \[Theta]2[t]}, {t}]; e2 = FullSimplify[e1[[1]]]; e3 = FullSimplify[e1[[2]]];  We build the section in two cases. First case $$\theta _1=0, \theta _1'>0$$ {sol, points} = Reap@NDSolve[{e2, e3, \[Theta]1[0] == TH10, \[Theta]2[0] == TH20, \[Theta]1'[0] == W10, \[Theta]2'[0] == W20, WhenEvent[\[Theta]1[t] == 0 && \[Theta]1'[t] > 0, Sow[{\[Theta]2[t], \[Theta]2'[t]}]]}, {\[Theta]1, \[Theta]2}, {t, 0, tm}] ListPlot[points, PlotRange -> All, Axes -> False, Frame -> True, FrameLabel -> {\[Theta], \[Theta]'}, PlotStyle -> Blue]  Second case $$\theta _2=0, \theta _2'>0$$ {sol, points} = Reap@NDSolve[{e2, e3, \[Theta]1[0] == TH10, \[Theta]2[0] == TH20, \[Theta]1'[0] == W10, \[Theta]2'[0] == W20, WhenEvent[\[Theta]2[t] == 0 && \[Theta]2'[t] > 0, Sow[{\[Theta]1[t], \[Theta]1'[t]}]]}, {\[Theta]1, \[Theta]2}, {t, 0, tm}] ListPlot[points, PlotRange -> All, Axes -> False, Frame -> True, FrameLabel -> {\[Theta], \[Theta]'}, PlotStyle -> Blue] ` • Wow, thanks for doing that! Is there any chance I could see the code you used to produce these plots? I still get a plot like the left even when I use the$\theta_2\$ condition. (see the edit to the original question) – David Jul 17 at 17:24
• @David See update to my answer. – Alex Trounev Jul 17 at 17:50