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:
