Runge-Kutta N-Body Implementation not working

Question

I've implemented in Java the Runge-Kutta 4 algorithm as described from this paper. http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf

However, the movements are erratic. Even with only two bodies, there are no periodic movements.

protected void integrateRK4(double time) {

    final double H = time;
    final double HO2 = H/2;
    final double HO6 = H/6;

    Vector2[] currentVelocities = new Vector2[objects.length];
    Vector2[] currentPositions = new Vector2[objects.length];
    Vector2[] vk1 = new Vector2[objects.length];
    Vector2[] vk2 = new Vector2[objects.length];
    Vector2[] vk3 = new Vector2[objects.length];
    Vector2[] vk4 = new Vector2[objects.length];
    Vector2[] rk1 = new Vector2[objects.length];
    Vector2[] rk2 = new Vector2[objects.length];
    Vector2[] rk3 = new Vector2[objects.length];
    Vector2[] rk4 = new Vector2[objects.length];


    for (int i=0; i<objects.length; i++) {
        currentVelocities[i] = objects[i].velocity().clone();
        currentPositions[i] = objects[i].position().clone();
    }

        vk1 = computeAccelerations(objects);
    for (int i=0; i<objects.length; i++) {
        rk1[i] = currentVelocities[i].clone();
    }

    for (int i=0; i<objects.length; i++) {
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(rk1[i], HO2)));
    }
        vk2 = computeAccelerations(objects);

    for (int i=0; i<objects.length; i++) {
        rk2[i] = Vectors2.prod(vk1[i], HO2);
    }


    for (int i=0; i<objects.length; i++) {
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(rk2[i], HO2)));
    }
        vk3 = computeAccelerations(objects);

    for (int i=0; i<objects.length; i++) {
        rk3[i] = Vectors2.prod(vk2[i], HO2);
    }


    for (int i=0; i<objects.length; i++) {
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(rk3[i], H)));
    }
        vk4 = computeAccelerations(objects);

    for (int i=0; i<objects.length; i++) {
        rk4[i] = Vectors2.prod(vk3[i], H);
    }


    for (int i=0; i<objects.length; i++) {
        objects[i].setVelocity(Vectors2.add(currentVelocities[i], Vectors2.prod(Vectors2.add(vk1[i], Vectors2.prod(vk2[i], 2), Vectors2.prod(vk3[i], 2), vk4[i]), HO6)));
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(Vectors2.add(rk1[i], Vectors2.prod(rk2[i], 2), Vectors2.prod(rk3[i], 2), rk4[i]), HO6)));
    }
}

Is my implementation incorrect somehow?

Note:

Vectors2 is my own implementation of Vectors, it is a first order Tensor of size 2.
All the static methods Vectors2.* return a copy of the vector.
All the non static methods called on instances of Vector2 modify the vector, same for objects[i].addVelocity() and objects[i].addPosition()

Vectors2.add(Vector2... vectors2) does element-wise addition.
Vectors2.prod(Vector2... vectors2) does element-wise multiplication.
Vectors2.prod(Vector2 vector2, double c) does element-wise multiplication.

computeAccelerations(PointBody[] objects) returns an array of acceleration vectors.
The variable objects is a property of the class NBodyIntegrator, which these functions are part of. It is an array of PointBody[].

To make sure that it wasn't some other bug, I've reduced the RK4 algorithm to a Explicit Euler algorithm by removing k2, k3, and k4. This one works as expected.

protected void integrateRK1(double time) {
    final double H = time;
    final double HO2 = H/2;

    Vector2[] currentVelocities = new Vector2[objects.length];
    Vector2[] currentPositions = new Vector2[objects.length];
    Vector2[] vk1;
    Vector2[] rk1 = new Vector2[objects.length];


    for (int i=0; i<objects.length; i++) {
        currentVelocities[i] = objects[i].velocity().clone();
        currentPositions[i] = objects[i].position().clone();
    }


        vk1 = computeAccelerations(objects);
    for (int i=0; i<objects.length; i++) {
        rk1[i] = currentVelocities[i].clone();
    }


    for (int i=0; i<objects.length; i++) {
        objects[i].setVelocity(Vectors2.add(currentVelocities[i], Vectors2.prod(Vectors2.add(vk1[i]), H)));
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(Vectors2.add(rk1[i]), H)));
    }
}

Consider a fourth-order [symplectic integrator](https://en.wikipedia.org/wiki/Symplectic_integrator), which will better preserve the total energy of the system. — David Eisenstat, Mar 17 '17 at 13:15
What are the implications of using runge-kutta vs symplectic integrators? I know that symplectic integrators are time-reversible but what are other advantages? — Bloc97, Mar 17 '17 at 13:35
What are the reasons people still use RK4 or even higher order RK algorithms? Are there some inherent advantages to RK? — Bloc97, Mar 17 '17 at 14:46

score 2 · Accepted Answer · answered Mar 17 '17 at 07:25

2

You are setting

rk1 = v0
pos2 = pos0 + rk1*h/2
vk2 = acc(pos2)

which is correct. But then you continue with

rk2 = vk1*h/2

where it should be

rk2 = v0 + vk1*h/2

and so on. Of course then the cumulative position update is also wrong.

answered Mar 17 '17 at 07:25

Lutz Lehmann

25,219
2
22
51

it should be rk3 = v0 + vk2 h/2 and rk4 = v0 + vk3 h/2 or rk3 = v1 + vk2 h/2 and rk4 = v2 + vk3 h/2 ? – Bloc97 Mar 17 '17 at 13:39
It is all offsets from the basis point, thus the first variant. The same as you did with the positions. – Lutz Lehmann Mar 17 '17 at 14:20
It is much better, but the orbits are still wrong, a binary system is stable in Explicit Euler, Midpoint and Leapfrog, but using RK4 it is not stable and the two objects go towards eachother. – Bloc97 Mar 17 '17 at 14:56

score -1 · Answer 2 · answered Mar 17 '17 at 09:11

-1

This is not a good implementation of RK integration.

You're using shared, mutable data. That's not thread safe.

Any proper implementation of numerical integration that I've ever seen abstracts the function to be integrated into a method that is passed into the integration scheme. You should be able to change integration schemes just by passing the function to a new routine.

Start with an interface for integration that takes a function and initial conditions as parameters.

answered Mar 17 '17 at 09:11

duffymo

305,152
44
369
561

I don't quite understand, all my properties are private, including Vector2, as it extends Tensor which has private properties. Vector2 has checks within everything to make sure it is a Tensor of order 1, size 2. – Bloc97 Mar 17 '17 at 13:43
objects isn't passed in; I assume it's a private, shared data member in the class. If it's mutable, it's not thread safe. – duffymo Mar 17 '17 at 14:38
I should transfer the integration to a static method in an abstract class? I kept objects[i] as a property of this integrator for performance reasons. – Bloc97 Mar 17 '17 at 14:48
Nope. I worry about your understanding of integration and objects. – duffymo Mar 17 '17 at 15:13
1

I take advice but I don't understand what I am doing wrong. Could you give an example or give a specific case where threading would cause problems? – Bloc97 Mar 17 '17 at 15:18

Runge-Kutta N-Body Implementation not working

2 Answers2