Double pendulum poincare section issues
I'm generating high-resolution Poincaré sections of a double pendulum, and I'm running into some issues in creating the initial conditions for a given section.
In general, I describe my pendulum with a state vector $(\theta_{1}, P_{\theta 1}, \theta_2, P_{\theta 2})$ for angles $\theta_1$, $\theta_2$ and momenta $P_{\theta 1}$, $P_{\theta 2}$. My prescription for generating initial conditions to get my Poincaré sections is:
1. Set the total energy $E$
2. Set $\theta_1 = P_{\theta 2} = 0$
3. Grid $\theta_2$ in a range from $-\pi$ to $\pi$
4. 394
Solve the Hamiltonian as a quadratic to determine $P_{\theta 1}$
This Hamiltonian is from this <PERSON> source, and was inverted in Mathematica to give
$$P_1 = \frac{1}{2} \left(\frac{P_2 \cos (\theta_1-\theta_2)}{\sin ^2(\theta_1-\theta_2)+1}-\sqrt{4 E+\frac{P_2^2 \cos ^2(\theta_1-\theta_2)}{\left(\sin ^2(\theta_1-\theta_2)+1\right)^2}-8 P_2^2+8 \cos (\theta_1)+4 \cos (\theta_2)}\right)$$
I think this is a fine procedure? There are a few arbitrary choices (like step 2; they could really be set to anything), but I'm not sure if this is doing anything necessarily bad.
However, there are a few different kinds of problems I'm encountering with this.
The first is for low-energy sections. The quadratic only has real solutions for certain ranges of $\theta_2$, and this precludes me from getting initial conditions which show the fixed points of the poincare section:
No matter what, the sections won't reduce to the fixed points because the Hamiltonian doesn't solve for those conditions. This figure is an attempt to recreate the figure from this SE article. It's getting the outermost sections, but the surfaces simply don't reduce to single points.
The second problem is with high-energy sections. I get these weird plots which look nothing like any poincare section I've seen in the literature:
Certainly, there's some structure there (but the convergence problem from before is still present), but there is a crazy amount of apparently chaotic behavior in the background. Shouldn't the chaotic behavior only be happening between the fixed points?
So, my questions are:
* Does my initial condition procedure seem okay?
* How can I generate the initial conditions for the fixed points?
* Does the second poincare section seem sane, or should I be looking elsewhere for a problem (like maybe my integrator is doing something weird or something)? 474