Three body problem in C++

Question

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):

Answer 1

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.