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,
```C
include
include
include
/ 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(argv1); TMAX = atof(argv2); TH10 = atof(argv3); W10 = atof(argv4); 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(NSTEPsizeof(float)); th1 = (float ) malloc(NSTEPsizeof(float)); w1 = (float ) malloc(NSTEPsizeof(float)); th2 = (float ) malloc(NSTEPsizeof(float)); w2 = (float ) malloc(NSTEPsizeof(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] = TH10PI/180.0; w1[0] = W10PI/180.0; th2[0] = TH20PI/180.0; w2[0] = W20PI/180. 935