5

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?

1
  • And yes, setting both bodies to same mass still generates velocity -.-
    – user2039981
    Commented May 11, 2014 at 11:59

1 Answer 1

24

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/

3
  • Hmm, in that do_physics, it is do_physics(1 / 8), shouldn't that be do_physics((1 / 60) * (1 / 8))
    – user2039981
    Commented Dec 17, 2014 at 15:45
  • Also where does that 100.0 come from?
    – user2039981
    Commented Dec 17, 2014 at 16:27
  • Gravity constant and time unit (and their relation to the length unit) are somewhat arbitrary in this simulation. It just works. 100 as the square of 10 gives a relatively small perturbation of the gravity potential. Its addition avoids divisions by zero or extremely large forces in exotic situations (while initializing a new disk). Commented Jan 27, 2015 at 22:23

Your Answer

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