Orbital modelling, three bodies, leapfrog numerical integration

Question

Apologies for the length of this question.

I'm modelling the motion of a single planet around the Sun using the equations:

$$dv_{x_{i}}/dt=a(t)=-GM_{SUN}x_{i}/r^{3},$$

where $x_{1}=x$, $x_{2}=y$ and $x_{3}=z$. Also,

$$r=\sqrt{x^{2}+y^{2}+z^{2}}.$$

I'm using the leapfrog numerical method of solving these equations (from Feynman, here), where:

$$x\left(t+\epsilon\right)=x\left(t\right)+\epsilon v\left(t+\epsilon/2\right)$$

$$v\left(t+\epsilon/2\right)=v\left(t-\epsilon/2\right)+\epsilon a\left(t\right)$$

$$v\left(\epsilon/2\right)=v\left(0\right)+\left(\epsilon/2\right)a\left(0\right).$$

Unable to code, I've put together a spreadsheet to model the orbit of Mars and (separately) Jupiter. I'm using units of au and Julian years. I'm using a time interval of $\epsilon = 6 \text{ hours} = 0.000684 \text{ yr}$. I've used (first row) initial state vector data $\left(x,y,z,v_{x},v_{y},v_{z}\right)$ from JPL Horizons, with an arbitrary start date of 2000-Jan-01 12:00:00.0000 TDB. $GM_{SUN} = 39.478418\,\text{au}^{3}\text{yr}^{-2}$. The results are pretty accurate for both Mars and Jupiter.

So far so good. The problem starts when I try to factor in Jupiter's gravitational effect on Mars' orbit by using the equations:

$$\frac{dv_{ix}}{dt}=\sum_{j=1}^{N}-\frac{GM_{j}\left(x_{i}-x_{j}\right)}{r_{ij}^{3}}$$

and so on for $y$ and $z$. Here, the index $i$ refers to Mars, and the index $j$ refers to the Sun and Jupiter. Also,

$$r=\sqrt{\left(x_{i}-x_{j}\right)^{2}+\left(y_{i}-y_{j}\right)^{2}+\left(z_{i}-z_{j}\right)}.$$

To keep things simple, I'm ignoring the gravitational effect of Jupiter on the Sun. I find the $x$ component of Mars' acceleration toward Jupiter at $t\left(0\right)$, which I'll call $MJa_{x}\left(0\right)$ using

$$MJa_{x}\left(0\right)=-GM_{JUPITER}\frac{\left(x_{MARS}\left(0\right)-x_{JUPITER}\left(0\right)\right)}{\left(r_{MARS-JUPITER}\left(0\right)\right)^{3}},$$

where $x_{MARS}\left(0\right)$ and $x_{JUPITER}\left(0\right)$ values are from JPL Horizons and $GM_{JUPITER}$ equals $0.037685\,au^{3}yr^{-2}$. And then, summing the accelerations, find

$$v_{x}\left(\epsilon/2\right)=v_{x}\left(0\right)+\left(\epsilon/2\right)\left(a_{x}\left(0\right)+MJa_{x}\right(0)),$$

where $v_{x}\left(0\right)$ is the $x$-component of Mars' velocity vector wrt the Sun (obtained from JPL), and $a_{x}\left(0\right)$ is the $x$-component of Mars' acceleration vector solely due to the Sun's gravitational attraction. The corrected position of Mars at $t\left(1\right)$ is then given by

$$x_{MARS}\left(1\right)=x_{MARS}\left(0\right)+\epsilon v_{x}\left(\epsilon/2\right).$$

I then find the next value of $x_{MARS}\left(1\right)-x_{JUPITER}\left(1\right)$ using my previously calculated $x_{JUPITER}\left(1\right)$ value (from the Jupiter orbit spreadsheet). Next, I repeat the above to find $MJa_{x}\left(1\right)$, the $x$ component of Mars' acceleration toward Jupiter at $t\left(1\right)$. And so on.

I was rather pleased with all this except for the fact that it doesn't work, i.e. the results are consistently less accurate than the original Sun-Mars model.

Any idea where I'm going wrong? Thanks.

EDIT - I've now amended my model by using $x_{MARS}-x_{JUPITER}$ values calculated directly from my existing Mars and Jupiter orbit spreadsheets (ie simply subtract $x_{JUPITER}$ from $x_{MARS}$. With those positions known, I then don't need to worry about calculating velocities as I can then use $x_{MARS}-x_{JUPITER}$ to directly find the accelerations $MJa_{x}$, etc. I've also taken into account the gravitational effect of Jupiter on the Sun as implied in the above $\frac{dv_{ix}}{dt}$ equation. The results are ambiguous at best. Am I doing this numerical integration correctly?

Answer 1

I don't have enough information to determine what is wrong with your spreadsheet, but I wrote the same program in Fortran, trying to follow your description. My results seem to be quite accurate, considering that Leapfrog is only second order and that we are neglecting Saturn, the Earth and the other planets. You can try and compare my results with yours, to understand what is the problem.

I chose to use Fortran because it is very readable, almost like writing math formulae directly, therefore you shouldn't have problems comparing my formulae with yours. Lines starting with an exclamation mark (!) are comments.

The program is the following.

program ThreeBody
implicit none
! UNITS used:
! Length: AU
! Time: Julian year

! Declare all the variables
integer, parameter:: dp=kind(0.d0)
! 3D position vector of Mars, the Sun and Jupiter
real(dp) :: r_m(3), r_s(3), r_j(3)
! 3D velocity vector of Mars, the Sun and Jupiter
real(dp) :: v_m(3), v_s(3), v_j(3)
! 3D acceleration vector of Mars, the Sun and Jupiter
real(dp) :: a_m(3), a_s(3), a_j(3)
! Mass of Jupiter, mass of the Sun
real(dp) :: GM_j = 0.037685, GM_s = 39.478418  ! au^3 yr^-2 
! Integration time step
real(dp) :: eps = 0.000684 ! 6 hours in Julian years
! index of current step, unit of the output file
integer :: i, outputFile

! Initial conditions from JPL Horizons 2000-Jan-01 12:00:00.0000 TDB
! Mars
r_m = [1.383579466402647E+00, 1.621204117825275E-02, -3.426152579691657E-02] ! au
v_m = [6.768779495222499E-04, 1.517984118259482E-02, 3.015574279972501E-04]  ! au/day
v_m = v_m*365.25 ! au/yr

! Jupiter
r_j = [3.994040256427818E+00, 2.935779784974931E+00, -1.015789849456240E-01]
v_j = [-4.562948121087195E-03, 6.435847855854325E-03, 7.548691257292337E-05]
v_j = v_j*365.25

! The Sun
r_s = [-7.137179161607904E-03, -2.795997495567612E-03, 2.062985061910257E-04]
v_s = [5.378460339618783E-06, -7.406916207973516E-06, -9.434293432137470E-08]
v_s = v_s*365.25

! open output file
open(newunit=outputFile, file="output.csv", status="replace")
! Write header to file
write(outputFile, '(A)', advance="no") "time,"
write(outputFile, '(A)', advance="no") "rx_m,ry_m,rz_m,rx_j,ry_j,rz_j,rx_s,ry_s,rz_s,"
write(outputFile, '(A)', advance="no") "vx_m,vy_m,vz_m,vx_j,vy_j,vz_j,vx_s,vy_s,vz_s,"
write(outputFile, '(A)')               "ax_m,ay_m,az_m,ax_j,ay_j,az_j,ax_s,ay_s,az_s"

! Main loop. Here we repeat the Leapfrog step 10000 times
do i = 1, 10000
    ! Calculate the acceleration
    ! Ignoring the effect of Mars on the Sun and on Jupiter
    a_m = -GM_j * (r_m - r_j)/norm2(r_m - r_j)**3 - GM_s * (r_m - r_s)/norm2(r_m - r_s)**3
    a_j = -GM_s * (r_j - r_s)/norm2(r_j - r_s)**3
    a_s = -GM_j * (r_s - r_j)/norm2(r_s - r_j)**3

    ! Calculate the velocity at the 1/2 step
    if (i == 1) then
        ! If we are at the first step, use eps/2 instead of eps
        v_m = v_m + eps/2 * a_m
        v_j = v_j + eps/2 * a_j
        v_s = v_s + eps/2 * a_s
    else
        v_m = v_m + eps * a_m
        v_j = v_j + eps * a_j
        v_s = v_s + eps * a_s
    end if
    
    ! Calculate the position
    r_m = r_m + eps * v_m
    r_j = r_j + eps * v_j
    r_s = r_s + eps * v_s

    ! Write the output to a csv file
    write(outputFile, '(*(G0.6,:,","))') i*eps, r_m, r_j, r_s, v_m, v_j, v_s, a_m, a_j, a_s
    
end do

close(outputFile)

end program

To run it, you can copy and paste it into a text file called ThreeBody.f. If you are on Linux, you just have to open a terminal in the directory where you have saved the file and run

gfortran ThreeBody.f && ./a.out

If you are on Windows or Mac, you will probably first need to install gfortran or any other compiler of your choice.

Alternatively, you can just Try it online!

The output of the program is a Comma Separated Value (CSV) file called output.csv that you can read as a spreadsheet with Excel or LibreOffice. It contains 28 columns:

• 1 time [yr]
• 2-4 3D position of Mars [au]
• 5-7 3D position of Jupiter [au]
• 8-10 3D position of the Sun [au]
• 11-19 same thing, but for the velocity [au/yr]
• 20-28 same thing, but for the acceleration [au/yr^2]

Just for reference, here is the 3D position of Mars for the first 10 steps of integration

And here it is after about 1 Julian year

Remarks on the accuracy of the result:

I don't know what accuracy you were expecting. In the comments I read that you considered $GM = 4\pi^2$ to be incorrect for your required level of accuracy. It best measured value can be found in the Wikipedia page for the Standard gravitational parameter, and, when expressed in units of Julian years and AU, it is equal to $4\pi^2$ to at least one part per 10000.

If your required level of accuracy is smaller than this, then you have to be careful that 6 hours are not $0.000684$ Julian years, but (with some more digits) $0.00068446269$ years. A better value for the gravitational parameter of the Sun would be $GM_\odot = 39.476926414 au^3 yr^{-2}$, and also the one of Jupiter could gain some more digits $GM_J = 0.037684447 au^3 yr^{-2}$. But, in any case, if you really need to know the position of Mars so precisely, then you probably need to include at least Saturn, the Earth and probably Venus in your calculation. They could have a far greater influence on the result, than the 6th decimal digit of the gravitational parameter.