RK4 integration of the three-bodies problem with C++

Question

first of all thank you for all the answers you gave me yesterday for the integration via Symplectic Euler's method of the three-body problem. We managed to implement both Euler's and Runge Kutta 4's method using C++, but the results are distinct. Since the three bodies move chaotically, we cannot understand which of the two methods is implemented correctly.

Below are the two pieces of code for Euler and RK4

#include <iostream>
#include <cmath>
#include <array>
#include <fstream>
#include "integrators.h"
static constexpr int DIM = 4;
static constexpr double G = 10;
static constexpr int N_BODIES = 3;
static constexpr int N_STEPS = 80000;
double distance(std::array<double, 3> r1, std::array<double, 3> r2){
    return sqrt(pow(r1[0]-r2[0],2)+pow(r1[1]-r2[1],2)+pow(r1[2]-r2[2],2));
}
class Planet{
public:
double m;
std::array &lt;double, 3&gt; x;
std::array &lt;double, 3&gt; v;
std::array &lt;double, 3&gt; a;
double energy;

Planet (double mass, double x_position, double y_position, double z_position, double x_velocity, double y_velocity, double z_velocity) {
    m = mass;
        x[0] = x_position;
        x[1] = y_position;
        x[2] = z_position;
        v[0] = x_velocity;
        v[1] = y_velocity;
        v[2] = z_velocity;
    };
    double getPositionX(void){
        return x[0];
    }
    double getPositionY(void){
        return x[1];
    }
    double getPositionZ(void){
        return x[2];
    }
    double getMass(void){
        return m;
    }
    double getVelocityX(void){
        return v[0];
    }
    double getVelocityY(void){
        return v[1];
    }
    double getVelocityZ(void){
        return v[2];
    }
    double getAccelX(void){
        return a[0];
    }
    double getAccelY(void){
        return a[1];
    }
    double getAccelZ(void){
        return a[2];
    }
    double getEnergy(void){
        return energy;
    }
    void computeEnergy(Planet planet1, Planet planet2){

        energy = 0.5 * m * (pow(v[0],2)+pow(v[1],2)+pow(v[2],2)) - G * m * planet1.m / distance(x, planet1.x) - G * m * planet2.m / distance(x, planet2.x);
    }

};
double acceleration(Planet A, Planet B, Planet C, int axe){
    //compute the acceleration along one axis of the body C
    double mass_A = A.getMass();
    double mass_B = A.getMass();
    double posx_A = A.getPositionX(); 
    double posx_B = B.getPositionX(); 
    double posx_C = C.getPositionX(); 
    double posy_A = A.getPositionY(); 
    double posy_B = B.getPositionY(); 
    double posy_C = C.getPositionY(); 
    double posz_A = A.getPositionZ(); 
    double posz_B = B.getPositionZ(); 
    double posz_C = C.getPositionZ();
    if (axe == 0)
        return (-1 * G * (mass_A * (posx_C-posx_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posx_C-posx_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
    else if (axe == 1){
        return (-1 * G * (mass_A * (posy_C-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posy_C-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
    }else if (axe == 2){
        return (-1 * G * (mass_A * (posz_C-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posz_C-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
    }
}
double F(double x, double v, double t, Planet A, Planet B, Planet C, int j ){
    double mass_A = A.getMass();
    double mass_B = B.getMass();
    double posx_A = A.getPositionX(); 
    double posx_B = B.getPositionX(); 
    double posx_C = C.getPositionX(); 
    double posy_A = A.getPositionY(); 
    double posy_B = B.getPositionY(); 
    double posy_C = C.getPositionY(); 
    double posz_A = A.getPositionZ(); 
    double posz_B = B.getPositionZ(); 
    double posz_C = C.getPositionZ();
if (j == 0){
    return (-1 * G * (mass_A * (x-posx_A) / pow(sqrt(pow(x-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posx_B) / pow(sqrt(pow(x-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 1) {
    return (-1 * G * (mass_A * (x-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(x-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(x-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 2) {
    return (-1 * G * (mass_A * (x-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(x-posz_A,2)), 3) + mass_B * (x-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(x-posz_B,2)), 3)));
}

}
int main(){
double h = 0.001;

Planet A(10, -10, 10, -11, -3, 0, 0);   // corpi allineati sull'asse delle x
Planet B(10, 0, 0, 0, 0, 0, 0);
Planet C(10, 10, 14, 12, 0, 0, 0);

std::array&lt;double, N_STEPS&gt; time;

double x_A[DIM][N_STEPS];
double x_B[DIM][N_STEPS];
double x_C[DIM][N_STEPS];
double vx_A;
double vy_A;
double vz_A;
double vx_B;
double vy_B;
double vz_B;
double vx_C;
double vy_C;
double vz_C;


double mass_A = A.getMass();
double mass_B = B.getMass();
double mass_C = C.getMass();


x_A[0][0] = A.getPositionX();
x_B[0][0] = B.getPositionX();
x_C[0][0] = C.getPositionX();

vx_A = A.getVelocityX();
vx_B = B.getVelocityX();
vx_C = C.getVelocityX();

x_A[1][0] = A.getPositionY();
x_B[1][0] = B.getPositionY();
x_C[1][0] = C.getPositionY();

vy_A = A.getVelocityY();
vy_B = B.getVelocityY();
vy_C = C.getVelocityY();

x_A[2][0] = A.getPositionZ();
x_B[2][0] = B.getPositionZ();
x_C[2][0] = C.getPositionZ();

vz_A = A.getVelocityZ();
vz_B = B.getVelocityZ();
vz_C = C.getVelocityZ();



//Function for the Euler method


// for (int i=0; i&lt;N_STEPS-1; i++){
//      for(int j=0; j&lt;DIM-1; j++){


//         A.a[j] = acceleration(B, C, A, j);
//         B.a[j] = acceleration(A, C, B, j);
//         C.a[j] = acceleration(B, A, C, j);

//         A.v[j] += A.a[j] * h;
//         B.v[j] += B.a[j] * h;
//         C.v[j] += C.a[j] * h;

//         x_A[j][i + 1] = x_A[j][i] + A.v[j] * h;
//         x_B[j][i + 1] = x_B[j][i] + B.v[j] * h;
//         x_C[j][i + 1] = x_C[j][i] + C.v[j] * h;

//         A.x[j] = x_A[j][i + 1];
//         B.x[j] = x_B[j][i + 1];
//         C.x[j] = x_C[j][i + 1];

//         A.computeEnergy(B, C);
//         B.computeEnergy(A, C);
//         C.computeEnergy(B, A);

//     } 

//     std::cout&lt;&lt;A.x[0]&lt;&lt;std::endl;
// }


// ---------------------------------------------------------------------------------------------
// RUNGE KUTTA 4

// Equation to Integrate
// x'' = -G(...)
// become
// x' = v
// v' = -G(...)

double m1;
double k1;
double m2;
double k2;
double m3;
double k3;
double m4;
double k4;

// initial conditions for the bodies
std::array&lt;double,3&gt; vA = A.v; 
std::array&lt;double,3&gt; xA = A.x;
std::array&lt;double,3&gt; vB = B.v; 
std::array&lt;double,3&gt; xB = B.x;
std::array&lt;double,3&gt; vC = C.v; 
std::array&lt;double,3&gt; xC = C.x;
double t;


for(int i=0; i&lt;N_STEPS-1; i++){
    for(int j=0; j&lt;DIM-1; j++){
        // body A
        m1 = h*vA[j];
        k1 = h*F(xA[j], vA[j], t, C, B, A, j);  

        m2 = h*(vA[j] + 0.5*k1);
        k2 = h*F(xA[j]+0.5*m1, vA[j]+0.5*k1, t+0.5*h, C, B, A, j);

        m3 = h*(vA[j] + 0.5*k2);
        k3 = h*F(xA[j]+0.5*m2, vA[j]+0.5*k2, t+0.5*h, C, B, A, j);

        m4 = h*(vA[j] + k3);
        k4 = h*F(xA[j]+m3, vA[j]+k3, t+h, C, B, A, j);

        xA[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
        vA[j] += (k1 + 2*k2 + 2*k3 + k4)/6;

        x_A[j][i+1] = xA[j];


        // body B
        m1 = h*vB[j];
        k1 = h*F(xB[j], vB[j], t, A, C, B, j);  //(x, v, t)

        m2 = h*(vB[j] + 0.5*k1);
        k2 = h*F(xB[j]+0.5*m1, vB[j]+0.5*k1, t+0.5*h, A, C, B, j);

        m3 = h*(vB[j] + 0.5*k2);
        k3 = h*F(xB[j]+0.5*m2, vB[j]+0.5*k2, t+0.5*h, A, C, B, j);

        m4 = h*(vB[j] + k3);
        k4 = h*F(xB[j]+m3, vB[j]+k3, t+h, A, C, B, j);

        xB[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
        vB[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
        x_B[j][i+1] = xB[j];


        // body C
        m1 = h*vC[j];
        k1 = h*F(xC[j], vC[j], t, A, B, C, j);  //(x, v, t)

        m2 = h*(vC[j] + 0.5*k1);
        k2 = h*F(xC[j]+0.5*m1, vC[j]+0.5*k1, t+0.5*h, A, B, C, j);

        m3 = h*(vC[j] + 0.5*k2);
        k3 = h*F(xC[j]+0.5*m2, vC[j]+0.5*k2, t+0.5*h, A, B, C, j);

        m4 = h*(vC[j] + k3);
        k4 = h*F(xC[j]+m3, vC[j]+k3, t+h, A, B, C, j);

        xC[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
        vC[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
        x_C[j][i+1] = xC[j];

        A.x[j] = xA[j];
        B.x[j] = xB[j];
        C.x[j] = xC[j];
    }

}



//----------------------------------------------------------------
// Print data on .csv
    std::ofstream output_file_A("positions_A.csv");
    std::ofstream output_file_B("positions_B.csv");
    std::ofstream output_file_C("positions_C.csv");
    output_file_A<<"x;y;z"<<std::endl;
    output_file_B<<"x;y;z"<<std::endl;
    output_file_C<<"x;y;z"<<std::endl;
for(int i = 0; i&lt;N_STEPS-1; i++){
    output_file_A &lt;&lt; x_A[0][i] &lt;&lt; &quot;;&quot; &lt;&lt; x_A[1][i] &lt;&lt; &quot;;&quot; &lt;&lt; x_A[2][i]&lt;&lt; std::endl;
    output_file_B &lt;&lt; x_B[0][i] &lt;&lt; &quot;;&quot; &lt;&lt; x_B[1][i] &lt;&lt; &quot;;&quot; &lt;&lt; x_B[2][i]&lt;&lt; std::endl;
    output_file_C &lt;&lt; x_C[0][i] &lt;&lt; &quot;;&quot; &lt;&lt; x_C[1][i] &lt;&lt; &quot;;&quot; &lt;&lt; x_C[2][i]&lt;&lt; std::endl;
}    
output_file_A.close();
output_file_B.close(); 
output_file_C.close();


return 0;
}

thanks for your help

Answer 1

General history, scalar ODE and systems

Up to order 4, the order conditions for explicit Runge-Kutta methods of scalar first-order equations are the same for the method for vector valued first-order systems, simply replacing the scalar state and derivative with the vector-valued state and derivative.

It is only in computing 5th order methods that the terms $f'(x)f''(x)[f(x),f(x)]$ and $f''(x)[f'(x)f(x),f(x)]$ and their associated order conditions, that are clearly different in the vector or system case, fall together in the scalar case. This detail led Kutta to design a 5-stage 5th order method that does not generalize to the vector case, where all 5th order methods require at least 6 stages.

On implementing coupled systems

To successfully implement the vector case of such a RK method while keeping the order conditions satisfied requires to treat the vectors and tangents as atomic objects, applying all vector operations to all coordinates at once. So if implementing \begin{align} \vec k_{1,x}&=\vec v,& \vec k_{1,v}&=\vec a(t,\vec x,\vec v)\\ \vec k_{2,x}&=\vec v+h/2·\vec k_{1,v},& \vec k_{2,v}&=\vec a(t+h/2,\vec x+h/2·\vec k_{1,x},\vec v+h/2·\vec k_{1,v}) \end{align} etc. you would first compute all forces/accelerations for all bodies, then compute an intermediary state of positions and velocities according to the second stage, then compute the forces/accelerations for this intermediary state and so on.

Variants on implementing composite systems

In my view, there are two design decisions on how to implement the state space. You can use a large flat vector and in each derivative computation read these data into the positions and velocities of the model instance, compute forces in the instance and then read the velocities and accelerations from the model instance into a flat derivatives vector.

The other variant is to treat the model as an abstract vector. The only operations that are needed are a copy constructor and the scaled addition of a tangent vector from a different instance, something like V.add_tangent(h,W) which would pass this operation from the system the bodies, and in the bodies would apply bv.x+=h*bw.v and bv.v+=h*bw.a.

For partitioned methods like the Verlet method one would have to provide separate update methods for position and velocities.

Reference solution in Python (as rapid prototyping)

The scipy solvers require the "one big flat array" variant, with a translation to the model data and back. For the present system, this can be implemented as

data

G=10
m=[10]*3
m=np.asarray(m)
x = [[-10,10,-11],[0,0,0],[10,14,12]]
v = [[-3,0,0],[0,0,0],[0,0,0]]
u0 = np.reshape([x,v],-1)
sys(0,u0)

ODE function

def sys(t,u):
    x,v = u.reshape([2,-1,3])
    dx = x[None,:,:]-x[:,None,:]
    r3 = np.sum(dx**2,axis=2)**1.5
    a = -G*np.sum(m[:,None,None]*dx/(1e-40+r3[:,:,None]), axis=0)
    return np.reshape([v,a],-1)

solver call

T=30; t = np.linspace(0,T,T*10+1)
res = solve_ivp(sys, t[[0,-1]], u0, t_eval=t, method='DOP853', atol=1e-8, rtol=1e-10)
print(res.message)

stereographic plot

vertical is the x-axis, horizontal y-axis, z-axis points "behind"