Chaotic Hamiltonian system poincare surfaces depend on the integrator

Question

First question on StackOverflow so go easy on me. I have a Hamiltonian system that consists on the following Hamiltonian:

$H(p,x;\textbf{P,X})=\frac{p^2}{2m}-a\frac{x^2}{2}+b\frac{x^4}{4}+x\sum_{n=1}^N g_n X_n+\sum_{n=1}^N \Big(\frac{P_n^2}{2M_n}+M_n \omega_n^2\frac{X_n^2}{2} \Big)$

Even for $n=1$ this is a chaotic system. The problem is that I have this Hamiltonian system that has some chaos, and I need to study it for a long time scale. I am having some problems integrating the equations of motion because the resulting trajectories depend on time step and integrator order, I thought that if I made Poincaré surface sections and they coincided then for the long term case different trajectories did not matter because the overall behavior of the system was defined by the Poincaré surfaces. But the big problem is that the Poincaré surfaces I got also depend on the type and order of the symplectic integrator I am using, so I am a little confused on what to trust now because different integrators give me different results for the same time scale.

I am currently using OrdinaryDiffEq.jl on Julia for the different symplectic integrators. Here is the plot of the Poincaré surface section for the case $X_1=0$ with positive momentum for two different types of symplectic integrators (Ruth3 and McAte4) with different initial conditions (with initial energy $E=15.0$) and parameters $m=1.$, $M_1=0.1$, $\omega_1=0.7071$, $g_1=0.05$, $a=0.25$ and $b=0.01$:

If I change the coupling constant to $g_1=0.01$ the difference is more notorious: The orbits differ a lot for this case.

For the case $g_1=0.0$ the curves do coincide:

edit: So the problem is in the harmonic oscillator side as different integrator's trajectories differ a little bit even thought they are symplectic, here I show the poincaré sections for $p=0$ for $X$ and $P$ for two different symplectic integrators:

Now I realize that for the Julia package consisting of different symplectic integrators (OrdinaryDiffEq.jl), the trajectories of the Harmonic Oscillator Hamiltonian differ a little bit from one another, here I show the phase portrait that illustrates my point:

Answer 1

A symplectic integrator can be understood as approximating the original system with a discrete Hamiltonian system with a Hamiltonian that should be close to the original one everywhere in phase space. This means, amongst other things, that the hypersurfaces of constant energy and other (even approximate) integrals of motion should be almost exactly matched in a Poincaré surface of section. This is obviously not the case in your second plot. There you see invariant structures crossing in a manner that is definitely not just due to numerical noise. This is a sign that either your equations of motion are not right, your integrator is not right, or you are not running the two integrators with initial conditions with the same energy.

To troubleshoot this I recommend to take the case $g_1 = 0$ where a single curve should emerge on the surface for every initial condition. If it does not, you are evolving the wrong equations or your integrator is badly written. From there you should be able to pin down the problem.