3rd Order Symplectic Integrator for N-Body Solar-System Problem

Question

I am doing a solar-system simulation. I am using Ruth's 3rd order sympletic integrator to avoid the problem of Energy Drift (which I had with RK4), but the the planets quickly leave orbit, and energy is by no means conserved (just like with RK4).

Here is Ruth's Integrator.

I applied this to the N-body problem here.

To get velocity, I just did: v=p/m

Have I correctly applied this algorithm to the gravitational N-body problem?

If you need more information on how I derived this, visit my question here.

Answer 1

The problem was twofold: 1) I should have used norm2 instead of abs when implementing in FORTRAN 2008.

i.e. This line: to_add=(rj-ri)big_gmasses(i)*masses(j)/((abs(rj-ri))**3.0d0) --> to_add=(rj-ri)big_gmasses(i)*masses(j)/((norm2(rj-ri))**3.0d0)

2. My data was messed up, so I tried switching to some one else's, but I forgot to change the units when I did.

So in summary, all the math and physics above should work.