2
$\begingroup$

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

(KE=1/2mv^2)

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)
    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.

Improve this question
$\endgroup$
2
  • 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, 2020 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$
    – Mauro Vanzetto
    Oct 19, 2020 at 10:25

1 Answer 1

Reset to default
0
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\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$
    – Lutz Lehmann
    Oct 19, 2020 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, 2020 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$
    – Lutz Lehmann
    Oct 19, 2020 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, 2020 at 2:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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