Euler integration of the three-body problem

Question

Me and one of my colleague are trying to simulate the three-body problem with a C++ program in order to compare different integration schemes. At the moment we are trying to use the Euler's method, but we're stuck on a divergence problem. Using the attached code we obtain the result in the figure below, in which the three body go away from each other. We are thinking about this problem from several days, so we ask if we are doing something wrong with the physics, or the integration.

We also don't know how to handle collisions between the bodies

#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 = 50000;

class Planet{

public:

    double m;
    std::array <double, 3> x;
    std::array <double, 3> v;
    std::array <double, 3> a;
    double energy;
    
    Planet () {m = 0; x = {0, 0, 0}; v = {0, 0, 0}; a = {0, 0, 0};};
        void setPlanet(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 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)));
    }
}

int main(){
    
   
    double h = 0.01;

    Planet A;
    Planet B;
    Planet C;

    A.setPlanet(10, -20, 20, 0, 10, 10, 0.1);   
    B.setPlanet(10, 0, 0, 0, 0, 0, 1);
    C.setPlanet(10, 20, -20, 0.3, -10, -10, 0);

    
    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();



    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);
            
            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.v[j] += A.a[j] * h;
            B.v[j] += B.a[j] * h;
            C.v[j] += C.a[j] * h;     
        
    }
}

    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;
    }   

   return 0;
      
}
```

Answer 1

What you see is a result of Euler integration.

One thing you can do is change to a symplectic integration which does not add energy to the system.

You can quickly try this out, but moving the update for the velocities before the update for the positions.

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];

}

Try this and see if it improves the situation.

---

See Figure 2.3, page 12 of "Geometric Numerical Integration Structure-Preserving Algorithms for Ordinary Differential Equations, Second Edition", by Ernst Hairer Christian Lubich Gerhard Wanner, for the numerical integration of the Kepler problem

Our next experiment (Fig. 2.3) studies the conservation of invariants and the global error. The main observation is that the error in the energy grows linearly for the explicit Euler method, and it remains bounded and small (no secular terms) for the symplectic Euler method. The global error, measured in the Euclidean norm, shows a quadratic growth for the explicit Euler compared to a linear growth for the symplectic Euler. As indicated in Table 2.1 the implicit midpoint rule and the Stormer–Verlet scheme behave similarly to the symplectic Euler, but have a smaller error due to their higher order. We remark that the angular momentum L(p,q) is exactly conserved by the symplectic Euler, the Stormer–Verlet, and the implicit mid- ¨ point rule.

Answer 2

I won't go into the script, but I'd suggest to rely on mathematical libraries, to use/test many other integration schemes (as an example Runge-Kutta schemes), only by changing the call to the integration function of the library.

I'm not using C++ for quite a long time, and maybe things have changed, but you could try to have a look at boost library.

I leave here some references:

• boost odeint page: https://www.boost.org/doc/libs/1_76_0/libs/numeric/odeint/doc/html/index.html
• some easy working examples here: http://boccelliengineering.altervista.org/junk/boost_integration/boost_odeint.html

Answer 3

Some napkin calculations to see if the "physical reality" is sufficiently different from the results shown:

The total mass $M$ is $30$, with $G=10$ one gets $\mu=GM=300$. For the only weakly related situation of a small planet in a circular orbit of radius about $R=30$ around a sun with $GM=300$, we get an angular speed $ω$ and period $T$ of $$ ω^2=\frac{GM}{R^3}=\frac{300}{27000}=\frac1{90},\\ T=\frac{2\pi}{ω}\approx60. $$ with tangential speed of $Rω=\sqrt{10}\approx 3.2$. The given initial speeds of $\sqrt{200}\approx 14$ are much larger than that and will correspond to very long elliptical or more likely hyperbolic orbits.

Now add the already mentioned short-comings of the Euler method and the presented plot is the expected outcome.

---

As an easy variation, modify the gravitation constant $G$. Integration with numerical methods of higher accuracy gives a longer term bounded state for $G=400$.

This variable-step method required step sizes down to $5·10^{-3}$, so with the Euler method one probably will need even smaller step sizes to replicate the character of these orbits.

Answer 4

Euler Integration won’t give you an exact solution for a differential integration, but only an approximation. You can improve the precision by taking smaller steps, or by using more precise formulas, like Simpson, Runge-Kutta, Fehlberg. But they all give you approximations.

Euler integration gives you an approximation that tends to systematically increase the energy in the system. It assumes that the bodies move on a straight line for some time instead of being slowed down by gravity. So use smaller steps or a better formula to improve your approximation.