# Trouble Making 3rd-Order Sympletic Integrator for Planitary N-Body Problem (A Hamiltonian System)

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). I applied this to the N-body problem with the following:   (KE=1/2mv^2)   I have implemented this into Fortran 2008, where x, a, v, p, and m are all vectors of length 30, which hold the x,y,z position, x,y,z acceleration, x,y,z velocity, x,y,z momentum, and m,m,m respectively for 10 separate bodies in the solar system (Planets + Sun + Pluto).

Acceleration on each body is calculated as the sum of a=GM/(r^2) for x,y,z for each other body on each body.

Here is the integration part of the code:

!----------Looping Through Time-----------
do while(t<365.250000d0) ! Length of simulation in days
!----------Calculating Values-----------
call calc_acc(masses,x,a)
p1=p+(7.0d0/24.0d0)*h*m*a
x1=x+(2.0d0/3.0d0)*h*p1/m
call calc_acc(masses,x1,a)
p2=p1+(3.0d0/4.0d0)*h*m*a
x2=x1-(2.0d0/3.0d0)*h*p2/m
call calc_acc(masses,x2,a)
p=p2-(1.0d0/24.0d0)*h*m*a
x=x2+h*p/m
v=p/m
t=t+h
!----------Saving Values-----------
do bodnum=1,10,1
write((100+bodnum),*) t, x((1+3*(bodnum-1)):(3+3*(bodnum-1))), v((1+3*(bodnum-1)):(3+3*(bodnum-1)))
write((200+bodnum),*) x((1+3*(bodnum-1))), x((2+3*(bodnum-1))), x((3+3*(bodnum-1)))
end do
end do


The full program can be found here.

Please tell me what I am doing wrong.

• I didn't have enough reputation to hyperlink this before, but here is Ruth's paper. And here is a page on Energy Drift. – maxbear123 Oct 19 at 2:01
• Welcome to scicomp! A tip: You can use MathJax to typeset your mathematical formulas. This will make the question much easier to read. – Mauro Vanzetto Oct 19 at 10:25

to_add=(rj-ri)*big_g*masses(i)*masses(j)/((abs(rj-ri))**3.0d0)

• I'd think that abs applied to a vector produces the vector of absolute values. Replacing with sum((rj-ri)**2)**1.5 does also not change the general divergence. – Lutz Lehmann Oct 19 at 18:43
• Yes, that works too. Check again the computation of big_g, using the formula in the comment I arrive at 2.670d-54. Also check all the signs of the initial conditions, in the earth data there is a minus missing in the position coordinates. – Lutz Lehmann Oct 19 at 20:26