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

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 <double, 3> x;
std::array <double, 3> v;
std::array <double, 3> 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<double, N_STEPS> 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<N_STEPS-1; i++){
//      for(int j=0; j<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<<A.x[0]<<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<double,3> vA = A.v;
std::array<double,3> xA = A.x;
std::array<double,3> vB = B.v;
std::array<double,3> xB = B.x;
std::array<double,3> vC = C.v;
std::array<double,3> xC = C.x;
double t;

for(int i=0; i<N_STEPS-1; i++){
for(int j=0; j<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<N_STEPS-1; i++){
output_file_A << x_A[0][i] << ";" << x_A[1][i] << ";" << x_A[2][i]<< std::endl;
output_file_B << x_B[0][i] << ";" << x_B[1][i] << ";" << x_B[2][i]<< std::endl;
output_file_C << x_C[0][i] << ";" << x_C[1][i] << ";" << x_C[2][i]<< std::endl;
}
output_file_A.close();
output_file_B.close();
output_file_C.close();

return 0;

}



• Your question doesn't include a question mark. Usually having one of those makes it clearer what you are asking. Nov 6, 2022 at 19:35
• Have you tested RK4 for something simple, e.g. $u' = \lambda u$? Nov 6, 2022 at 19:38
• Test both methods on a problem for which you know the solution. If the two methods give similar results, test them on your problem, on a short time interval with small timesteps, to check if the results are similiar as well. Nov 6, 2022 at 20:43
• Your RK4 is wrong. You can not apply it per component/coordinate, you have to apply it to all components/coordinates at once. Meaning the m1,k1,... should have type array<double,3*3> to contain the data for all coordinates of all bodies. Or you need k1Ax,...,k4Cz,m1Ax,...,m4Cz. A better implementation would totally separate the physical model from the numerical method. The method does not need to know the interna of the physical model, and the model should be agnostic of the method used. The interface are the vector operations on the state vectors. Nov 7, 2022 at 5:29
• In the first segment the bodies B and C are in free fall towards each other, slightly perturbed by the distant body A that also moves away. This results in a near-collision shortly before time $t=8$. The outcome of this is very sensible to numerical errors, like insufficient sampling and step errors. Every difference after $t=8$ is more caused by the difference in exit directions than by further numerical errors. Nov 7, 2022 at 8:27

### 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"

• Incredible answers from you. I hope students somewhere are getting the benefit of such deep expertise. Nov 21, 2022 at 20:00