What is wrong with my gravity simulation?

Question

As per advice given to me in this answer, I have implemented a Runge-Kutta integrator in my gravity simulator.

However, after I simulate one year of the solar system, the positions are still off by cca 110 000 kilometers, which isn't acceptable.

My initial data was provided by NASA's HORIZONS system. Through it, I obtained position and velocity vectors of the planets, Pluto, the Moon, Deimos and Phobos at a specific point in time.

These vectors were 3D, however, some people told me that I could ignore the third dimension as the planets aligned themselves in a plate around the Sun, and so I did. I merely copied the x-y coordinates into my files.

This is the code of my improved update method:

"""
Measurement units:

[time] = s
[distance] = m
[mass] = kg
[velocity] = ms^-1
[acceleration] = ms^-2
"""

class Uni:
    def Fg(self, b1, b2):
        """Returns the gravitational force acting between two bodies as a Vector2."""

        a = abs(b1.position.x - b2.position.x) #Distance on the x axis
        b = abs(b1.position.y - b2.position.y) #Distance on the y axis

        r = math.sqrt(a*a + b*b)

        fg = (self.G * b1.m * b2.m) / pow(r, 2)

        return Vector2(a/r * fg, b/r * fg)

    #After this is ran, all bodies have the correct accelerations:
    def updateAccel(self):
        #For every combination of two bodies (b1 and b2) out of all bodies:
        for b1, b2 in combinations(self.bodies.values(), 2):
            fg = self.Fg(b1, b2) #Calculate the gravitational force between them

            #Add this force to the current force vector of the body:
            if b1.position.x > b2.position.x:
                b1.force.x -= fg.x
                b2.force.x += fg.x
            else:
                b1.force.x += fg.x
                b2.force.x -= fg.x


            if b1.position.y > b2.position.y:
                b1.force.y -= fg.y
                b2.force.y += fg.y
            else:
                b1.force.y += fg.y
                b2.force.y -= fg.y

        #For body (b) in all bodies (self.bodies.itervalues()):
        for b in self.bodies.itervalues():
            b.acceleration.x = b.force.x/b.m
            b.acceleration.y = b.force.y/b.m
            b.force.null() #Reset the force as it's not needed anymore.

    def RK4(self, dt, stage):
        #For body (b) in all bodies (self.bodies.itervalues()):
        for b in self.bodies.itervalues():
            rd = b.rk4data #rk4data is an object where the integrator stores its intermediate data

            if stage == 1:
                rd.px[0] = b.position.x
                rd.py[0] = b.position.y
                rd.vx[0] = b.velocity.x
                rd.vy[0] = b.velocity.y
                rd.ax[0] = b.acceleration.x
                rd.ay[0] = b.acceleration.y

            if stage == 2:
                rd.px[1] = rd.px[0] + 0.5*rd.vx[0]*dt
                rd.py[1] = rd.py[0] + 0.5*rd.vy[0]*dt
                rd.vx[1] = rd.vx[0] + 0.5*rd.ax[0]*dt
                rd.vy[1] = rd.vy[0] + 0.5*rd.ay[0]*dt
                rd.ax[1] = b.acceleration.x
                rd.ay[1] = b.acceleration.y

            if stage == 3:
                rd.px[2] = rd.px[0] + 0.5*rd.vx[1]*dt
                rd.py[2] = rd.py[0] + 0.5*rd.vy[1]*dt
                rd.vx[2] = rd.vx[0] + 0.5*rd.ax[1]*dt
                rd.vy[2] = rd.vy[0] + 0.5*rd.ay[1]*dt
                rd.ax[2] = b.acceleration.x
                rd.ay[2] = b.acceleration.y

            if stage == 4:
                rd.px[3] = rd.px[0] + rd.vx[2]*dt
                rd.py[3] = rd.py[0] + rd.vy[2]*dt
                rd.vx[3] = rd.vx[0] + rd.ax[2]*dt
                rd.vy[3] = rd.vy[0] + rd.ay[2]*dt
                rd.ax[3] = b.acceleration.x
                rd.ay[3] = b.acceleration.y

            b.position.x = rd.px[stage-1]
            b.position.y = rd.py[stage-1]

    def update (self, dt):
        """Pushes the uni 'dt' seconds forward in time."""

        #Repeat four times:
        for i in range(1, 5, 1):
            self.updateAccel() #Calculate the current acceleration of all bodies
            self.RK4(dt, i) #ith Runge-Kutta step

        #Set the results of the Runge-Kutta algorithm to the bodies:
        for b in self.bodies.itervalues():
            rd = b.rk4data
            b.position.x = b.rk4data.px[0] + (dt/6.0)*(rd.vx[0] + 2*rd.vx[1] + 2*rd.vx[2] + rd.vx[3]) #original_x + delta_x
            b.position.y = b.rk4data.py[0] + (dt/6.0)*(rd.vy[0] + 2*rd.vy[1] + 2*rd.vy[2] + rd.vy[3])

            b.velocity.x = b.rk4data.vx[0] + (dt/6.0)*(rd.ax[0] + 2*rd.ax[1] + 2*rd.ax[2] + rd.ax[3])
            b.velocity.y = b.rk4data.vy[0] + (dt/6.0)*(rd.ay[0] + 2*rd.ay[1] + 2*rd.ay[2] + rd.ay[3])

        self.time += dt #Internal time variable

The algorithm is as follows:

1. Update the accelerations of all bodies in the system
2. RK4(first step)
3. goto 1
4. RK4(second)
5. goto 1
6. RK4(third)
7. goto 1
8. RK4(fourth)

Did I mess something up with my RK4 implementation? Or did I just start with corrupted data (too few important bodies and ignoring the 3rd dimension)?

How can this be fixed?

---

Explanation of my data etc...

All of my coordinates are relative to the Sun (i.e. the Sun is at (0, 0)).

./my_simulator 1yr

Earth position: (-1.47589927462e+11, 18668756050.4)

HORIZONS (NASA):

Earth position: (-1.474760457316177E+11,  1.900200786726017E+10)

I got the 110 000 km error by subtracting the Earth's x coordinate given by NASA from the one predicted by my simulator.

relative error = (my_x_coordinate - nasa_x_coordinate) / nasa_x_coordinate * 100
               = (-1.47589927462e+11 + 1.474760457316177E+11) / -1.474760457316177E+11 * 100
               = 0.077%

The relative error seems miniscule, but that's simply because Earth is really far away from the Sun both in my simulation and in NASA's. The distance is still huge and renders my simulator useless.

Answer 1

110 000 km absolute error means what relative error?

I got the 110 000 km value by subtracting my predicted Earth's x coordinate with NASA's Earth x coordinate.

I'm not sure what you're calculating here or what you mean by "NASA's Earth x coordinate". That's a distance from what origin, in what coordinate system, at what time? (As far as I know, the earth moves in orbit around the sun, so its x-coordinate w.r.t. a coordinate system centered at the sun is changing all the time.)

In any case, you calculated an absolute error of 110,000 km by subtracting your calculated value from "NASA's Earth x coordinate". You seem to think this is a bad answer. What's your expectation? To hit it spot on? To be within a meter? One km? What's acceptable to you and why?

You get a relative error by dividing your error difference by "NASA's Earth x coordinate". Think of it as a percentage. What value do you get? If it's 1% or less, congratulate yourself. That would be quite good.

You should know that floating point numbers aren't exact on computers. (You can't represent 0.1 exactly in binary any more than you can represent 1/3 exactly in decimal.) There are going to be errors. Your job as a simulator is to understand the errors and minimize them as best you can.

You could have a stepsize problem. Try reducing your time step size by half and see if you do better. If you do, it says that your results have not converged. Reduce by half again until you achieve acceptable error.

Your equations might be poorly conditioned. Small initial errors will be amplified over time if that's true.

I'd suggest that you non-dimensionalize your equations and calculate the stability limit step size. Your intuition about what a "small enough" step size should be might surprise you.

I'd also read a bit more about the many body problem. It's subtle.

You might also try a numerical integration library instead of your integration scheme. You'll program your equations and give them to an industrial strength integrator. It could give some insight into whether or not it's your implementation or the physics that causes the problem.

Personally, I don't like your implementation. It'd be a better solution if you'd done it with mathematical vectors in mind. The "if" test for the relative positions leaves me cold. Vector mechanics would make the signs come out naturally.

UPDATE:

OK, your relative errors are pretty small.

Of course the absolute error does matter - depending on your requirements. If you're landing a vehicle on a planet you don't want to be off by that much.

So you need to stop making assumptions about what constitutes too small a step size and do what you must to drive the errors to an acceptable level.

Are all the quantities in your calculation 64-bit IEEE floating point numbers? If not, you'll never get there.

A 64 bit floating point number has about 16 digits of accuracy. If you need more than that, you'll have to use an infinite precision object like Java's BigDecimal or - wait for it - rescale your equations to use a length unit other than kilometers. If you scale all your distances by something meaningful for your problem (e.g., the diameter of the earth or the length of the major/minor axis of the earth's orbit) you might do better.