vector - N-Body Gravity Simulation in JavaScript

Question

So, I am trying to create a N-Body Gravity simulation in JavaScript:

http://jsfiddle.net/4M94x/

var Circle = function(c, r, cor, cof) { // Fix CoR & CoF // Had to add code for JSFiddle link :P
    this.c = c
    this.r = r
    this.m = r * r * Math.PI
    this.v = new Vector()
    this.cor = cor
    this.cof = cof
}

The problem's that when you spawn (click) and put 2 balls (accidentally renamed "particles") next to each other they start to generate velocity and faster and faster push eachother.. How do I fix this, btw is my gravity implementation correct?

Answer 1

This is easy to explain: You are implementing Euler forward as solver for the ODE, and in mechanical systems the systematic error of Euler forward increases the energy. Euler backward decreases the energy, so a combination of alternating the explicit and implicit Euler methods would leave the energy a little more constant.

But then you can implement with the same or even less effort the second order symplectic methods which preserve energy and other physical invariants, either the (implicit) midpoint method or the Verlet-(Stoermer-Cromer-...-Newton-)method.

Or even higher order Runge-Kutta, which will also preserve the energy to a higher order, despite not being symplectic.

See Hairer on the Stoermer-Verlet-...-Newton method, postprint or preprint and the "Moving stars around" tutorial text using C++ or Ruby.

---

A note to the physics: All in all the implementation is very minimal and well readable. But the gravitational force is

g*m1*m2*(p2-p1)/norm(p2-p1)^3

as the negative gradient of

g*m1*m2/norm(p2-p1)

you are only using the square of the norm, where the force would be the negative gradient of the gravitational potential energy

 g*m1*m2*ln(norm(p2-p1))

which is right for flatland physics, but not for a 2D section of 3D space.

---

Working code

with velocity Verlet and energy preservation:

Add a new field a=Vector() to the circle object and replace the kitchen sink which is the update() function with the following collection of dedicated functions

function compute_forces() {
    for (var i = 0; i < particles.length; i++) {
        var p = particles[i];
        p.a.set(0);

        for (var j = 0; j < i; j++) {
            var p2 = particles[j];

            var d = p.c.sub(p2.c);
            var norm = Math.sqrt(100.0 + d.lengthSq());
            var mag = gravity / (norm * norm * norm);

            p.a.set(p.a.sub(d.mul(mag * p2.m)));
            p2.a.set(p2.a.add(d.mul(mag * p.m)));

        }
    }

}


function do_collisions() {
    for (var i = 0; i < particles.length; i++) {
        var p = particles[i];
        for (var j = 0; j < i; j++) {
            var p2 = particles[j];

            if (checkCCCol(p, p2)) {
                resCCCol(p, p2);
            }
        }
    }
}


function do_physics(dt) {
    // do velocity Verlet 
    // with consistent state at interval half points
    // x += 0.5*v*dt
    for (var i1 = 0; i1 < particles.length; i1++) {
        var p1 = particles[i1];
        p1.c.set(p1.c.add(p1.v.mul(0.5 * dt)));
    }
    // a = A(x)
    compute_forces();
    // v += a*dt
    for (var i2 = 0; i2 < particles.length; i2++) {
        var p2 = particles[i2];
        p2.v.set(p2.v.add(p2.a.mul(dt)));
    }
    // x += 0.5*v*dt
    for (var i3 = 0; i3 < particles.length; i3++) {
        var p3 = particles[i3];
        p3.c.set(p3.c.add(p3.v.mul(0.5 * dt)));
    }
    do_collisions();
}

function update() {

    for (var k = 0; k < 4; k++) {
        do_physics(1.0 / 4);
    }

    render();

    RAF(update);
}

See http://jsfiddle.net/4XVPH/

---

Altered example with particles coloured based on their mass(hopefully better displaying their interaction), one bug fixed, and some additional comments: http://jsfiddle.net/24mg6ctg/12/