# 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 '20 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 '20 at 10:25

What exactly do you think the formula

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


does, especially the denominator? For the correct physics it should be the third power of the Euclidean distance.

• 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 '20 at 18:43
• I think that abs(vector) is [|x|,|y|,|z|], while norm2(vector) is Euclidian Distance of the vector – maxbear123 Oct 19 '20 at 19:44
• 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 '20 at 20:26
• Vielen Dank. The problem was twofold: 1) I should have used norm2 instead of abs. 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. – maxbear123 Oct 20 '20 at 2:59