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 with the following:

enter image description here

enter image description here

enter image description here


enter image description here

enter image description here

enter image description here

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)
    call calc_acc(masses,x1,a)
    call calc_acc(masses,x2,a)
    !----------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.

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

What exactly do you think the formula


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

  • $\begingroup$ 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. $\endgroup$ Oct 19 '20 at 18:43
  • $\begingroup$ I think that abs(vector) is [|x|,|y|,|z|], while norm2(vector) is Euclidian Distance of the vector $\endgroup$
    – maxbear123
    Oct 19 '20 at 19:44
  • $\begingroup$ 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. $\endgroup$ Oct 19 '20 at 20:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – maxbear123
    Oct 20 '20 at 2:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.