# Chaotic Hamiltonian system poincare surfaces depend on the integrator

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:

• Are you performing sanity checks, such as computing the energy along the trajectories to confirm it's constant (and the one you chose initially), or that the values that define the Poincaré section are being respected? – stafusa Feb 27 '20 at 0:19
• So what is your actual question? It is to be expected that you will get different results with different interrogators since the are all approximations taking finite steps in what should be a continuum. You mentioned step size. This is a critical parameter in any ODE solver and either needs to be fixed small enough to never cause problems or regulated after each step with some sort of convergence check and an update. Are you using canned solvers in MATLAB, Maple, or did you write your own? I'd recommend reading numerical recipes. – ggcg Feb 28 '20 at 15:25

## 1 Answer

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.

• Thank you for your comment. When I check with $g_1=0$ the curves do coincide, that is why I am so confused. I edited the post to include that plot. – mooD Feb 26 '20 at 19:27
• @mooD Could you also try plotting a surface of section defined by $p = 0$? I.e. could you plot $P,X$ in a section? The problem with the equations of motion could still be hidden in the $P,X$ sector since it is decoupled at this point. – Void Feb 27 '20 at 17:04
• You are totally right, the problem is in the $P,X$ sector, I realized that the trajectories of the harmonic oscillator depend on the type of symplectic integrator I use. My guess is that it is not totally wrong but you have to be aware that some integrator's trajectories change a little bit. I edited the post again to show what I mean by this. – mooD Feb 28 '20 at 2:27
• @mooD There are basically two possibilities. 1) your integration time steps are too large (the time step should be a small fraction of the period of the integrator and the differences should become negligible for decreasing time steps), or 2) there is a coding error. 2) could consist of a mistake on your side, or of a mistake in the package you are using. – Void Feb 28 '20 at 19:04