Calculating the planets and moons based on Newtons's gravitational force

Question

I am trying to understand how to calculate the orbits of solar system bodies in an n-body framework, based on the pair-wise gravitational interaction between the objects. At present, I am considering 44 objects (sun, planets, major moons and major asteroids).

I am starting with the state vectors (position and velocity) of each of the objects with Sun as centre obtained from telnet ssd.jpl.nasa.gov 6775 (JPL Horizons) 01-Jan-2017 at 00:00 UTC and would like to let the system evolve for 4344h, 01-July-2017 at 00:00h.

I have written a program to do this in Java, and so far the results do not seem to be even reasonably close to what they should be, comparing with the state vectors obtained from Horizons. After every 2 second time step the net gravitational forces on each body from all of the others is calculated, and then in one shot all the velocities and positions are updated based on the accelerations from those net forces. Then I compare the final updated position vectors from the application with data obtained from Horizons after correcting for the Sun's updated position.

The comparison shows that the positions of Earth and the outer planets are have position error of less than 50km (In fact, the farther planets it is less then 10km). Where as for Mercury the error is 250km. And the moons of Jupiter and Saturn are off by 50,000 to 300,000 km!

In my application, I am not differentiating between the Sun, planets and moons, so I am not sure why there should be so much more error for the moons. I have tried decreasing the step size, from 2 seconds down to 0.25 seconds, but there is not significant improvement.

What might be the problems that I should investigate here? Are there things that clearly need improvement right away? Or perhaps there are tests I can to do help diagagnose the primary sources of error?

---

EDIT: Here is the gist of the calculation method as requested in comments:

while (nowT < endT) {
    doOneStep(step, nowT)
    nowT += stepT
}

allItemLinks is collections of ItemLink - links between objects. in this case Gravity link between all object pairs. For n objects, there will be n.(n+1)/2 links

doOneStep(double deltaT, double nowT){
    initForces fo all items to 0,0,0
    for each ItemLink **allItemLinks**)
        inf.evalForce(deltaT, false)
    updatePosAndVel(deltaT, nowT, true)
}

In ItemLink:

evalForce(double deltaT, boolean bFinal) {
    addGravityEffect(deltaT);
}

boolean addGravityEffect(double deltaT) {
    rVector = item2.pos - item1.pos
    double gF = G.mq.m2/r2
    fVector = gF in rVector direction
    item1.addForce(fVector)
    similarly for item2 to item1
}

allItems is a collection of Item objects (Sun, planets and moons)

void updatePosAndVel(double deltaT, double nowT) {
    for each Item of **allItems** updatePandV(deltaT, nowT);
}

In Item:

netForce, nowAcc, effectiveAcc, deltaV, newPos etc. are all Vector3d

updatePAndV(double deltaT, double nowT, boolean bFinal){
        nowAcc = netForce / mass
        effectiveAcc = mean of lastAcc and nowAcc
        deltaV = effectiveAcc * deltaT
        meanV ...
        newPos = oldPos + meanV * deltaT
}

Not working on gravitation fields but with direct forces due to inter-object gravity.

With the above code, I am able to get stable orbits. Even the orbit times of the moons are good. Get nice Saturn set with the Cycloidal motions of moons and the Uranus set with helical motion of the moons around Uranus. I do not know as to how to send files or images for this discussion

Answer 1

Aside from numerical issues, "With Sun as centre" may be part of your problem. Get all the data from Horizons relative to the Solar System Barycenter, not the Sun, which moves relative to the barycenter. That barycenter is an inertial frame of reference, whereas the center of the Sun is not. Also make sure that you are putting in the initial position and velocity of the Sun and letting it move, if you aren't already.

Answer 2

I'll post the numerical methods in this answer, and the full-blown solar system calculation (including relativity and possible effects from the Sun's oblate shape) as a second answer. There's just too much to put it all in one answer.

The method you describe using looks like the Euler Method or what can be called Euler Forward method. At each time step $dt$ you calculate all of the forces and their resulting net accelerations $\mathbf{a}$, and then just increment the velocities $\mathbf{v}$ by $\mathbf{a} dt$ and all positions $\mathbf{x}$ by $\mathbf{v} dt$. You need absurdly tiny time steps to get even close with this. You mentioned 2 seconds and 250 milliseconds, a bit short for solar system timescales.

In the script below, I've written the Euler Forward method and two more Runge-Kutta low order methods usually called RK2 and RK4. For each method a simplified (fake) Mercury orbit around a fixed Sun is calculated for 100 days, with the number of iterations varied from 10 to 10,000. Also, for each, I use the SciPy library solver odeint() with a relative precision tolerance of 1E-12 per step. odeint is a python wrapper for lsoda from the FORTRAN library odepack.

You can see that even RK4 agrees with odeint to the level of meters after 100 days if you use a time step of about 15 minutes, and the Euler Forward method (what you use) will need an absurd number of steps to even approach this.

Numerical techniques are not the only issue here, I will post a second answer in a few days with the rest of what you need. I may do it under a separate question rather than post two answers to the same question.

But this should get you started either coding RK4 in java, or looking for a java numerical library like the Apache one mentioned in comments.

The simplest standard way to solve an orbital problem using an ODE solver is to put all of the cartesian coordinates and velocities in one long vector, call it $y$, and write a single function for the time derivative:

$$\dot{y} = f(t, y)$$

So if you had two bodies and three dimensions, the vector $y$ would be:

$$ y=(x_1, y_1, z_1, x_2, y_2, z_2, v_{x1}, v_{y1}, v_{z1}, v_{x2}, v_{y2}, v_{z2})$$

having six elements per body. The derivative of $x_1$ is $v_{x1}$, and the derivative of $v_{x1}$ is the acceleration $a_{x1}$ due to all of the other bodies.

Say you have one body in a central force with a standard gravitational parameter $GM$. The rate of change of position is just the velocity,

$$\frac{d\mathbf{x}}{dt} = \mathbf{v}$$

and the rate of change of velocity is the acceleration, due to the force;

$$\frac{d\mathbf{v}}{dt} = \mathbf{a} = -GM \frac{\mathbf{r}}{r^3}$$

If you had multiple bodies, the $\mathbf{r}$ would be the distance between pairs of bodies and for each body you would sum over all of the others as you've described already.

So if you wrote out $f$, it would be:

$$ f=(v_{x}, v_{y}, v_{z}, -GMx/r^3, -GMy/r^3, -GMz/r^3).$$

The Euler Forward method that I believe you are using just iterates

$$y_{i+1} = h \ f(t, y_i) $$

with time steps of $h$. The improved RK2 method (shown below) would be written:

$$k_1 = f(t, y_i)$$ $$k_2 = f(t+\frac{h}{2}, y_i+\frac{h}{2}k_1)$$ $$y_{i+1} = y_n + h k_2 $$

and the ubiquitous RK4 method is written as

$$k_1 = f(t, y_i)$$ $$k_2 = f(t+\frac{h}{2}, y_i+\frac{h}{2}k_1)$$ $$k_3 = f(t+\frac{h}{2}, y_i+\frac{h}{2}k_2)$$ $$k_4 = f(t+h, y_i+k_3)$$

$$y_{i+1} = y_n + h (k_1 + 2k_2 + 2k_3 + k_4)/6 $$

Here is the section that shows the similarity and differences between the Euler Forward method (the simplest Runge-Kutta method) and two higher order RK methods. This is written for clarity, not speed obviously.

def deriv(t, X):

    x, v = X.reshape(2, -1)
    acc  = -GMs * x * ((x**2).sum())**-1.5

    return np.hstack((v, acc))

def derivv(X, t):
    return deriv(t, X) # historical reasons

def RK_all(F, t0, X0, n, h, method=None):

    hov2, hov6   = h/2.0, h/6.0
    t, X         = t0, X0
    answer, time = [], []

    answer.append(X)
    time.append(t)

    for i in range(n):
        k1 = F(t,        X            )
        k2 = F(t + hov2, X + hov2 * k1)
        k3 = F(t + hov2, X + hov2 * k2)
        k4 = F(t + h,    X + h    * k3)

        if   method == 'EulerFwd':
            X = X + h*k1        # X + h*F(t, X)

        elif method == 'RK2':
            X = X + h*k2

        elif method == 'RK4':
            X = X + hov6*(k1 + 2.*(k2+k3) + k4)

        else:
            pass

        t += h

        answer.append(X)
        time.append(t)

    return np.array(time), np.array(answer)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

GMs = 1.327E+20 # approx

X0 = np.array([-2.094E+10,  4.303E+10,  5.412E+09,
               -5.328E+04, -2.011E+04,  3.243E+03]) # approx

methodnames = ['RK4', 'RK2', 'EulerFwd']
niterations = [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000]

Time = 100*24*3600.  # total time

methdict = dict()

for methodname in methodnames:

    times, answers, ODEanswers, posdevs = [], [], [], []

    for n in niterations:

        h  = Time/float(n)
        t0 = 0.0

        time, answer    = RK_all(deriv, t0, X0, n, h, method=methodname)
        # recalculate using library ODE solver for same times, to compare
        ODEanswer, info = ODEint(derivv, X0, time,
                                 rtol=1E-12, full_output=True)
        posdev = np.sqrt((((answer - ODEanswer)[:,:3])**2).sum(axis=1))


        times.append(time)
        answers.append(answer)
        ODEanswers.append(ODEanswer)
        posdevs.append(posdev)

    methdict[methodname] = (times, answers, ODEanswers, posdevs)

if 1 == 1:
    plt.figure()
    for i, meth in enumerate(methodnames):
        plt.subplot(1, 3, i+1)
        for time, answer, ODEanswer, posdev in zip(*methdict[meth]):
            x, y, z = answer.T[:3]
            plt.plot(x, y)
        plt.ylim(-2.8E+11, 0.8E+11)
        plt.xlim(-1.2E+11, 0.8E+11)
        plt.title(meth, fontsize=16)
        plt.plot([0],[0], 'ok')
    plt.show()

if 1 == 1:
    plt.figure()
    for i, meth in enumerate(methodnames):
        plt.subplot(1, 3, i+1)
        for time, answer, ODEanswer, posdev in zip(*methdict[meth]):
            plt.plot(time/(24*3600.), posdev)
        plt.yscale('log')
        plt.ylim(1E-01, 1E+12)
        plt.title(meth+' vs odeint', fontsize=16)
    plt.suptitle('RKmethod - odeint (meters) vs time (days)', fontsize=18)
    plt.xticks([0, 20, 40, 60, 80, 100])
    plt.show()