# Three body problem in C++

I am in a begginers programming course and we got a little project. I chose to simulate the three body problem using the Euler method. Even though the system is chaotic there are some special cases which produce symmetric trajectories. I get most of them, but the most important one doesn't come out. It's called figure 8 and supposed to be stable.

I'm sure that the initial conditions are fine. And the program should be as well, as I get the other expected results. So the question is whether there is an error in my program or the Euler method doesn't work for this specific case. (I would then change to the Runge-Kutta method)

I appreciate any kind of help!

#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

double const  ms=1.0000000000000000;
double const  me=1.0000000000000000;
double const  mm=1.0000000000000000;
double const  G=1.000000000000000000;
double const  dt=0.0001000000000000;

/* 1=Sonne,2=Erde, 3=Mond */

double fsx(double x1, double y1, double x2, double y2, double x3, double y3){
double dv;
dv=G*(me/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((x2-x1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(mm/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((x3-x1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)));
/*Projektion der Gravitationskraft auf Körper 1 auf x-Richtung */
return dv;
}

double fsy(double x1, double y1, double x2, double y2, double x3, double y3){
double dv=G*(me/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((y2-y1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(mm/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((y3-y1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1))); /*Projektion der Gravitationskraft auf Körper 1 auf y-Richtung */
return dv;
}

double fex(double x1, double y1, double x2, double y2, double x3, double y3){
double dv;
dv=G*(ms/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((x2-x1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(mm/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((x3-x1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)));
/*Projektion der Gravitationskraft auf Körper 1 auf x-Richtung */
return dv;
}

double fey(double x1, double y1, double x2, double y2, double x3, double y3){
double dv=G*(ms/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((y2-y1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(mm/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((y3-y1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1))); /*Projektion der Gravitationskraft auf Körper 1 auf y-Richtung */
return dv;
}

double fmx(double x1, double y1, double x2, double y2, double x3, double y3){
double dv;
dv=G*(me/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((x2-x1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(ms/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((x3-x1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)));
/*Projektion der Gravitationskraft auf Körper 1 auf x-Richtung */
return dv;
}

double fmy(double x1, double y1, double x2, double y2, double x3, double y3){
double dv=G*(me/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((y2-y1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(ms/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((y3-y1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1))); /*Projektion der Gravitationskraft auf Körper 1 auf y-Richtung */
return dv;
}

double euler(double (f)(double,double,double,double,double,double), double x1, double y1, double x2, double y2, double x3, double y3, double vi){
double vf=vi + f(x1, y1, x2, y2, x3, y3)*(dt);
return vf;
}

int main() {
float xf3, yf3, tmax;
ofstream out("project1.dat");

ti=0.00000000000000000, xi1=-0.97000436, yi1=0.24308753, vxi1=0.466203685, vyi1=0.43236573, tmax=3.000000000000, xi2=-xi1, vxi2=vxi1, yi2=-yi1, vyi2=vyi1, xi3=0.0000000000000, vxi3=-2*vxi1, yi3=0.0000000000000000, vyi3=-2*vyi1;
for(tf=0.0;tf<=tmax;tf=tf+dt){
vxf1=euler(fsx,xi1,yi1,xi2,yi2,xi3,yi3,vxi1);
vyf1=euler(fsy,xi1,yi1,xi2,yi2,xi3,yi3,vyi1);
vxf2=euler(fex,xi2,yi2,xi1,yi1,xi3,yi3,vxi2); /* Gravitationskraft und damit auch fx und fy sind symmetrisch -> es reichen 2 Funktionen */
vyf2=euler(fey,xi2,yi2,xi1,yi1,xi3,yi3,vyi2);
vxf3=euler(fmx,xi3,yi3,xi2,yi2,xi1,yi1,vxi3);
vyf3=euler(fmy,xi3,yi3,xi2,yi2,xi1,yi1,vyi3);
xf1=xi1+(vxf1)*(dt);
yf1=yi1+(vyf1)*(dt);
xf2=xi2+(vxf2)*(dt);
yf2=yi2+(vyf2)*(dt);
xf3=xi3+(vxf3)*(dt);
yf3=yi3+(vyf3)*(dt);

out << tf << " "<< xf1 << " " << yf1 <<" "<< vxf1 << " " << vyf1<<  " "<< xf2 << " " << yf2 <<" "<< vxf2 << " " << vyf2<< " "<< xf3 << " " << yf3 <<" "<< vxf3 << " " << vyf3 << endl;
ti=tf;
vxi1=vxf1;
vyi1=vyf1;
vxi2=vxf2;
vyi2=vyf2;
vxi3=vxf3;
vyi3=vyf3;

xi1=xf1;
yi1=yf1;
xi2=xf2;
yi2=yf2;
xi3=xf3;
yi3=yf3;
}
out.close();
}


The grafic (it's supposed to be a perfect 8):

• You should expect Euler to drift for this type of problem due to phase and amplitude error problems with explicit Euler. Range-Kutta has better properties as it pertains to amplitude and phase errors, but you might also want to look into some integration schemes that are symplectic integrators. – spektr May 19 '16 at 1:25
• Thanks a lot for the fast answer @choward! Actually I found some program where they were using the Euler method (not in c++) and apparently they got the expected results. Thas why I was hoping it would work for me in c++ as well. The thing is that we should solve our problem by using the methods we have learned in class. But you haven´t seen any other error in the code, right? – Maxi May 19 '16 at 2:02
• Ah I see. Well maybe a bug in your code exists then. Will have to investigate! – spektr May 19 '16 at 2:07
• Can you write the equations you used for your iterations? – nicoguaro May 19 '16 at 3:23
• Without even looking over your code, I have an advice: lose the Euler method and use a symplectic integrator, like leapfrog or verlet. – Adrian Roman May 20 '16 at 11:04

So something must be up with your code. I know this because when I just wrote a C++ code to solve this problem with my simulation framework, I got the following result using Explicit Euler:

I would recommend breaking frequent computations into functions, especially functions that can do vector related math. You might consider representing the positions and velocities of each body with an array, like:

double p1[2], v1[2]; // These represent position and velocity of body 1


Then you could write functions to compute a vector magnitude and a unit vector like:

double hypot( double x, double y ){
double t; // modified based on code from wikipedia on "Hypot" page
x = abs(x);
y = abs(y);
if( x == 0.0 ){ return y; }
if( y == 0.0 ){ return x; }
t = min(x,y);
x = max(x,y);
t = t/x;
return x*sqrt(1.0+t*t);
}

double mag( double * vec ){
return hypot(vec[0],vec[1]);
}

void toUnit( double * vec ){
double mag_ = mag(vec);
vec[0] /= mag_; vec[1] /= mag_;
}


And use the functions by simply doing:

double r[2] = {p1[0]-p2[0],p1[1]-p2[1]};
double magr = mag(r);
toUnit(r);
double accel_mag = G*m1/(magr*magr);
double dv2dt[2] = {accel_mag*r[0],accel_mag*r[1]};


Getting a habit of building frequently used computations into independent functions will not only make it take less time to build new codes, but debugging your problems will also be easier since if you have a mistake in some of the math, you only have to fix maybe one function even though this function is used in many places.

• Instead of your "mag" function, it might be better to use "hypot" which is more accurate for cases where the components are very large or very small. – Lysistrata May 19 '16 at 13:01
• @Lysistrata yeah I forgot about the "hypot" function, so good thought. I will add that in. – spektr May 19 '16 at 13:45