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

boolean addGravityEffect(double deltaT) {
    rVector = item2.pos - item1.pos
    double gF = G.mq.m2/r2
    fVector = gF in rVector direction
    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

  • 1
    $\begingroup$ The offset of Mercury would be about right, due to relativistic effects. I'd blame the inaccuracy of the moons positions on numerical accuracy errors, especially that they interact with each other a lot. $\endgroup$
    – SF.
    Sep 5, 2017 at 9:13
  • 1
    $\begingroup$ Also, if you calculate their position relative to the Sun yo For Io, you have semi-major axis of 0.003 AU orbital radius, at distance of 5.2 AU from the Sun. This works poorly with floating point values as the fixed, large offset from center of the system of coordinates prevents them to increase precision by changing the mantisse and changes occur at the far tail of the base with a lot of precision lost behind its tail. $\endgroup$
    – SF.
    Sep 5, 2017 at 9:24
  • 1
    $\begingroup$ Calculation duration in 6 months (01 Jan 2017 to 01 July 2017). $\endgroup$
    – vishu
    Sep 5, 2017 at 14:59
  • 1
    $\begingroup$ I am coding in Java with Vector3d ('double' variables). Suspecting precision problem, tried with Big Decimal(32 significant digits) for positions of Europa of Jupiter (only one object with BigDecimal since it takes too long to calculate). Still the error for Europa was very large similar to earlier calculation with 'double'. I am not using any integration methods - just numerical additions, updating positions and velocity every 2 seconds. $\endgroup$
    – vishu
    Sep 5, 2017 at 15:58
  • 1
    $\begingroup$ @uhoh Thanks for moving my additional data to the main question. $\endgroup$
    – vishu
    Sep 6, 2017 at 11:36

2 Answers 2


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.

  • $\begingroup$ Good point! I didn't even notice that. $\endgroup$
    – uhoh
    Sep 5, 2017 at 23:03
  • $\begingroup$ I am taking data from horizons with Sun-centre at 20170101 and let the objects react and move. Initially position and velocity of Sun are both 000. During calculation The sun does move. After completing the calculations till 20170701, all the resultant data are modified back with respect to the new position and velocity of the Sun as base. This is then compared with data from Horizons for 20170701 again with Sun-centre. $\endgroup$
    – vishu
    Sep 6, 2017 at 11:35
  • $\begingroup$ What I am recommending is that you get all of your data relative to the Solar System Barycenter (@0). Not the Sun (@10). If you do it correctly, the Sun will not have an initial position and velocity of zero. Then you don't need to rebase to the new Sun position and velocity at the end. $\endgroup$
    – Mark Adler
    Sep 6, 2017 at 12:01
  • $\begingroup$ But is it not that Barycenter itself is dependant on the relative positions of the planet. I can understand that Earth-Moon Barycenter is fixed relative to Earth and Moon. But with multiple planets dodn't it change with planets' postions. $\endgroup$
    – vishu
    Sep 6, 2017 at 12:27
  • 3
    $\begingroup$ No. That's the whole point of finding and using the barycenter. It is the center of mass of all of the solar system bodies, and so by conservation of momentum, it does not move in the inertial frame of a closed solar system, regardless of how the bodies in the solar system move. $\endgroup$
    – Mark Adler
    Sep 6, 2017 at 13:21

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.

enter image description here

enter image description here

enter image description here

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 = [], []


    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)


        t += h


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


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

if 1 == 1:
    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')

if 1 == 1:
    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.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])
  • $\begingroup$ @vishu I will keep adding to this, if something isn't clear please leave a mesasge. It is good practice for me to write this all out - sometimes explaining clearly can be a challenge. Anybody Else please comments and suggested improvements on clarity are welcome! $\endgroup$
    – uhoh
    Sep 10, 2017 at 12:00
  • $\begingroup$ @vishu is this looking at all helpful? Adding General Relativity to $f$ will correct the problem with Mercury and it's just a few more lines of code actually - I'll try to add it somehow in the next day or so. But to get everything working you need to either implement something like this or use a standard ODE solver in Java. Please feel free to add comments, or ask me to rewrite this or add more if it helps. $\endgroup$
    – uhoh
    Sep 13, 2017 at 8:13
  • $\begingroup$ Now I have fairly good understanding what you are saying (learning Python as a side effect :) ). I will try out with RK4 after 25th September, since I will not abel tuse my laptop for the next 10 days. $\endgroup$
    – vishu
    Sep 13, 2017 at 10:40
  • $\begingroup$ @vishu OK that's really great! Running python the way I show here is very very slow, you would not normally do it this way, but it's an easy way to "prototype". If you use an integration library like scipy.odeint or 'scipy.ode` it's better. $\endgroup$
    – uhoh
    Sep 13, 2017 at 10:46
  • 1
    $\begingroup$ I was taking average acceleration with the previous step, instead with the next future step would have been similar to Modified Euler. With the small increment time (2s and 250 ms) should have partly offset the inaccuracy, but looks like not enough for the orbits with higher curvature of the moons, which I understood from your explanations. Thanks for that. I will try out with RK4 after 25th September $\endgroup$
    – vishu
    Sep 13, 2017 at 10:52

Your Answer

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

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