Integration by Example - Euler vs Verlet vs Runge-Kutta

Aug. 28, 2010
comments

In this post I'm testing different integration methods for a gravity simulation. The results can be inspected interactively in the canvas tags that accompany each test. Hover with the mouse over the illustration to start its simulation or click the illustration to reset the simulation.

I'm putting emphasis on gravity here because that's what's this post is about. Different integration methods exist, and for different equations they can be very accurate. I am focusing on relatively cheap and fast integration schemes and how their accuracies in a newtonian gravity situation compare.

Contents

Helpers

The "acceleration" method computes the acceleration between two positions. It does not include mass in the calculation for simplicty.

var G = 1500.0;
var acceleration = function(pos1, pos2){
    var direction = pos1.sub(pos2);
    var length = direction.length();
    var normal = direction.normalized();
    return normal.mul(G/Math.pow(length, 2));
};

In order to make the example code a little simpler, an object is provided called "OneBody". It expects as a first parameter the #id of the canvas object to draw to, the body that's simulated and a stepping function.

new OneBody('euler', {
    body: {}, // the simulated body
    step: function(center, body){
        // stepping here
    }
});

You can download the implementation of the helpers for yourself

Euler Integration

So the Euler integration is fairly simple to do. At each step you compute the force, add it to the velocity and add the velocity to the position. There's of course an integration error, as is evidenced by the orbital path that's not running in the same track.

new OneBody('euler', {
    body: {
        position: new Vec(70, 60),
        velocity: new Vec(0, 1)
    },
    step: function(center, body){
        var acc = acceleration(center, body.position);
        body.velocity.iadd(acc);
        body.position.iadd(body.velocity);
    }
});

If you decrease the time step you can get a more accurate Euler integration. Here I'm choosing to loop 8 times for each step, and the time step consequently is 1/8th.

var steps = 8;
var delta = 1/steps;
for(var i=0; i<steps; i++){
    var acc = acceleration(center, body.position);
    body.velocity.iadd(acc.mul(delta));
    body.position.iadd(body.velocity.mul(delta));
}

Verlet Integration

This integration scheme does not store a velocity but computes it on the fly. It has advantages over Euler integration when dealing with multiple bodies and constraints. Gotoandplay has a nice writeup of the method.

For one body however, it is identical to Euler Integration.

new OneBody('verlet', {
    body: {
        position: new Vec(70, 60),
        previous: new Vec(70, 59)
    },
    step: function(center, body){
        var acc = acceleration(center, body.position);

        var position = body.position.mul(2);
        position.isub(body.previous);
        position.iadd(acc);

        body.previous = body.position;
        body.position = position;
    }
});

Runge-Kutta Method

This integration method was proposed by C. Runge and M.W. Kutta. A nice introduction is supplied by gafferongames. It's way more complex then Euler or Verlet integration. In its basic form it also seems to be very inaccurate, way more inaccurate then Euler integration.

var Derivative = function(position, velocity){
    this.position = position ? position : new Vec(0, 0);
    this.velocity = velocity ? velocity : new Vec(0, 0);
    
    this.iadd = function(other){
        this.position.iadd(other.position);
        this.velocity.iadd(other.velocity);
        return this;
    }
    this.add = function(other){
        return new Derivative(
            this.position.add(other.position),
            this.velocity.add(other.velocity)
        )
    }

    this.mul = function(scalar){
        return new Derivative(
            this.position.mul(scalar),
            this.velocity.mul(scalar)
        )
    }
    this.imul = function(scalar){
        this.position.imul(scalar);
        this.velocity.imul(scalar);
        return this;
    }
}
var compute = function(center, initial, delta, derivative){
    var state = derivative.mul(delta).add(initial);

    return new Derivative(
        state.velocity,
        acceleration(center, state.position)
    );
}
var d0 = new Derivative();

new OneBody('rk4', {
    body: {
        position: new Vec(70, 60),
        velocity: new Vec(0, 1)
    },
    step: function(center, body){
        delta = 1;
        var d1 = compute(center, body, delta*0.0, d0);
        var d2 = compute(center, body, delta*0.5, d1);
        var d3 = compute(center, body, delta*0.5, d2);
        var d4 = compute(center, body, delta*1.0, d3);
        
        d2.iadd(d3).imul(2);
        d4.iadd(d1).iadd(d2).imul(1/6);
        
        body.position.iadd(d4.position.mul(delta));
        body.velocity.iadd(d4.velocity.mul(delta));
    }
});

However if you loop RK4 four times it is more accurate then running Euler 8 times. It also aproximates the true elipse form a bit better.

steps = 4;
delta = 1/steps;
for(var i=0; i<steps; i++){
    var d1 = compute(center, body, delta*0.0, d0);
    var d2 = compute(center, body, delta*0.5, d1);
    var d3 = compute(center, body, delta*0.5, d2);
    var d4 = compute(center, body, delta*1.0, d3);
    
    d2.iadd(d3).imul(2);
    d4.iadd(d1).iadd(d2).imul(1/6);
    
    body.position.iadd(d4.position.mul(delta));
    body.velocity.iadd(d4.velocity.mul(delta));
}

Verdict

I'm less then impressed with Verlet integration and even RK4 seems a bit lacking at least in the first iteration compared to simply running Euler more. The performance of zero loops RK4 is especially puzzling since it gets better than Euler integration quickly with more iterations. Euler, Verlet and RK4 integration have different performance for other differential equations (then the Newtonian physics), so try them out!

I'd like to thank Geoff for the help in debugging the code, and if any errors still remain, bug me about them and I will adjust the code and post.