Last time we went over the orbits and everything nice that comes with it. Whilst implementing the code for those calculations I’ve observed something peculiar. The orbits the objects are following are more or less the same as the ones calculated in the initial state, but not quite.
Naturally, I went debugging, and everything seemed fine. Still, the bodies would slowly drift out of their orbits. Trying different things I implemented a quick time warp functionality and observed the cause of all the new issues, time.

As you are probably aware by now, games run in discrete intervals. Commonly known as framerate or FPS, the small time step between frames (dt) is the root of all problems when it comes to accurate physics. You might’ve seen code similar to the one below when applying velocity to position and acceleration to velocity.

position += velocity * dt;
velocity += acceleration *dt;

This is called “explicit Euler” integration. Named after the same Euler that named that “e” constant that I can never remember what is for. The problem with this is that whilst it is easy to understand and explain, it is wrong. Intuitively, it makes sense that if I have a body moving at ten kilomteres per hour, after a delta time (dt) of one hour, the body would have moved ten kilometers. And that is correct, as long as the velocity is constant. If the velocity varies say, due to an acceleration being applied, then the larger the delta time, the bigger the error.

Let’s consider the following example to hopefully clarify why this is. Assume we have a stationary body with a unitary mass. Let’s apply an acceleration of ten meters per second per second for ten seconds and see what we get.

    t=0:    position = 0      velocity = 0
    t=1:    position = 0      velocity = 10
    t=2:    position = 10     velocity = 20
    t=3:    position = 30     velocity = 30
    t=4:    position = 60     velocity = 40
    t=5:    position = 100    velocity = 50
    t=6:    position = 150    velocity = 60
    t=7:    position = 210    velocity = 70
    t=8:    position = 280    velocity = 80
    t=9:    position = 360    velocity = 90
    t=10:   position = 450    velocity = 100

Based on this calculation, it looks like the body would have travelled 450 meters after 10 seconds. If you were to apply the formula for movement under constant acceleration it would become obvious that that is wrong since the actual position should be 500 meters away after 10 seconds.
Essentially, what needs to be pointed out here is that the error has less to do with the length of the simulation than it does with the delta time at which it is run.

We can do better. To get around the error we can use the Velocity Verlet integration method. This method will yield the correct result, regardless of the delta time at which the simulation is run, as long as the acceleration is constant. You can find below the code that updates the position and velocity accordingly.

const Vector3 Acceleration = new Vector3(1,0,0);

Vector3 Velocity;
Vector3 OldVelocity;
Vector3 Position;

public void Update(float dt)
{
    OldVelocity = Velocity;
    Velocity += Acceleration * dt;
    Position += (OldVelocity + Velocity) * 0.5f * dt;
}

That is much better but we still have a problem. Our acceleration is not constant, not even with circular orbits. So far, we’ve covered a first order integrator, Euler, for when velocity is contant and a second order integrator, Velocity Verlet, for when acceleration is constant. What can we do if nothing is constant.

Every single solution I found has drawbacks. Higher order simplectic integrators lack accuracy whilst the RK4 (Runge-Kutta Fourth Order) integrator doesn’t conserve energy which means our orbits will decay in time.
In the end I decided to go with RK4 for now. See here a great article on RK4 for the implementation details. That article will have everything you need to know in order to be able to write your own RK4 implementation and there’s no point in me repeating it. The only tricky bit for me was to remember to use the Velocity Verlet integration for velocity into position after using RK4 to integrate the acceleration in order to get accurate values.

See below a short video illustrating correct and accurate orbit motion happening. The physics engine is running 15 times per second, with a time warp factor of 5.