5
$\begingroup$

I wrote the following program to solve the three body problem (Sun, Earth, Jupiter), initializing the system such that the total angular momentum is 0. The result of the simulation is really bad (and the Sun tend to infinity). I looked at the code for more than an hour trying to find the problem. Does anybody see what is wrong with my program?

Thank you very much.

#include <iostream>
using std::cout;
using std::endl;

#include <fstream>
using std::ofstream;

#include <cmath>
using std::sqrt;

int main()
{
    double time(0), dt(0.05);

    double G(1.9838e-29);

    double m1(6e24), m2(1.9e27), ms(1.99e30);

    double x1(1), y1(0), x2(5.2), y2(0), xs(0), ys(0);

    double T1(1), T2(11.8618);

    double vx1(0),vy1(2*M_PI*x1/T1), vx2(0), vy2(2*M_PI*x2/T2);

    double Rx( (m1*x1+m2*x2+ms*xs)/(m1+m2+ms) ), Ry( (m1*y1+m2*y2+ms*ys)/(m1+m2+ms) );

    x1 = x1 - Rx;
    x2 = x2 - Rx;
    xs = xs - Rx;
    y1 = y1 - Ry;
    y2 = y2 - Ry;
    ys = ys - Ry;

    double Lx1(x1*m1*vx1), Lx2(x2*m2*vx2);

    double vxs(0), vys( -(Lx1+Lx2)/(ms*xs) );

    double r12(0), r1s(0), r2s(0);

    ofstream output("TB_jupiter-earth.dat");

    output << time << ' ' << x1 << ' ' << y1 << ' ' << x2 << ' ' << y2 << ' ' << xs << ' ' << ys <<endl;

    while(time < 12)
    {
        r12 = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
        r1s = sqrt( (x1-xs)*(x1-xs) + (y1-ys)*(y1-ys) );
        r2s = sqrt( (x2-xs)*(x2-xs) + (y2-ys)*(y2-ys) );

        vx1 = vx1 - G*ms*dt*(x1-xs)/(r1s*r1s*r1s) - G*m2*(x1-x2)*dt/(r12*r12*r12);
        vy1 = vy1 - G*ms*dt*(y1-ys)/(r1s*r1s*r1s) - G*m2*(y1-y2)*dt/(r12*r12*r12);

        vx2 = vx2 - G*ms*dt*x2/(r2s*r2s*r2s) - G*m1*(x2-x1)*dt/(r12*r12*r12);
        vy2 = vy2 - G*ms*dt*y2/(r2s*r2s*r2s) - G*m1*(y2-y1)*dt/(r12*r12*r12);

        vxs = vxs - G*m1*(xs-x1)*dt*(r1s*r1s*r1s) - G*m2*(xs-x2)*dt*(r2s*r2s*r2s);
        vys = vys - G*m1*(ys-y1)*dt*(r1s*r1s*r1s) - G*m2*(ys-y2)*dt*(r2s*r2s*r2s);


        x1 = x1 + vx1*dt;
        y1 = y1 + vy1*dt;

        x2 = x2 + vx2*dt;
        y2 = y2 + vy2*dt;

        xs = xs + vxs*dt;
        ys = ys + vys*dt;

        time += dt;

        output << time << ' ' << x1 << ' ' << y1 << ' ' << x2 << ' ' << y2 << ' ' << xs << ' ' << ys <<endl;
    }

    output.close();

    return 0;
}

EDIT:

The correction of the first answer were good but now I changed the calculations of the angula momentum (they were wrong) and the program doesn't work properly anymore. The code is the following:

#include <iostream>
using std::cout;
using std::endl;

#include <fstream>
using std::ofstream;

#include <cmath>
using std::sqrt;

int main()
{
    double time(0), dt(0.005);

    double G(1.9838e-29);

    double m1(6e24), m2(500*1.9e27), ms(1.99e30);

    double x1(1), y1(0), x2(5.2), y2(0), xs(0), ys(0);

    double T1(1), T2(11.8618);

    double vx1(0),vy1(2*M_PI*x1/T1), vx2(0), vy2(2*M_PI*x2/T2);

    double Rx( (m1*x1+m2*x2+ms*xs)/(m1+m2+ms) ), Ry( (m1*y1+m2*y2+ms*ys)/(m1+m2+ms) );

    x1 = x1 - Rx;
    x2 = x2 - Rx;
    xs = xs - Rx;
    y1 = y1 - Ry;
    y2 = y2 - Ry;
    ys = ys - Ry;

    double Ly1(x1*m1*vy1), Ly2(x2*m2*vy2);

    double vxs(0), vys( -(Ly1+Ly2)/(ms*xs) );

    double r12(0), r1s(0), r2s(0);

    ofstream output("TB_jupiter-earth.dat");

    output << time << ' ' << x1 << ' ' << y1 << ' ' << x2 << ' ' << y2 << ' ' << xs << ' ' << ys <<endl;

    cout << x1 << ' ' << y1 << ' '<< vx1 << ' ' << vy1 << endl;
    cout << x2 << ' ' << y2 << ' '<< vx2 << ' ' << vy2 << endl;
    cout << xs << ' ' << ys << ' '<< vxs << ' ' << vys << endl;

    while(time < 4)
    {
        r12 = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
        r1s = sqrt( (x1-xs)*(x1-xs) + (y1-ys)*(y1-ys) );
        r2s = sqrt( (x2-xs)*(x2-xs) + (y2-ys)*(y2-ys) );

        vx1 = vx1 - G*ms*dt*(x1-xs)/(r1s*r1s*r1s) - G*m2*(x1-x2)*dt/(r12*r12*r12);
        vy1 = vy1 - G*ms*dt*(y1-ys)/(r1s*r1s*r1s) - G*m2*(y1-y2)*dt/(r12*r12*r12);

        vx2 = vx2 - G*ms*dt*(x2-xs)/(r2s*r2s*r2s) - G*m1*(x2-x1)*dt/(r12*r12*r12);
        vy2 = vy2 - G*ms*dt*(y2-ys)/(r2s*r2s*r2s) - G*m1*(y2-y1)*dt/(r12*r12*r12);

        vxs = vxs - G*m1*(xs-x1)*dt/(r1s*r1s*r1s) - G*m2*(xs-x2)*dt/(r2s*r2s*r2s);
        vys = vys - G*m1*(ys-y1)*dt/(r1s*r1s*r1s) - G*m2*(ys-y2)*dt/(r2s*r2s*r2s);


        x1 = x1 + vx1*dt;
        y1 = y1 + vy1*dt;

        x2 = x2 + vx2*dt;
        y2 = y2 + vy2*dt;

        xs = xs + vxs*dt;
        ys = ys + vys*dt;

        time += dt;

        output << time << ' ' << x1 << ' ' << y1 << ' ' << x2 << ' ' << y2 << ' ' << xs << ' ' << ys <<endl;
    }

    output.close();

    return 0;
}
Improve this question
$\endgroup$
  • $\begingroup$ "initializing the system such that the total angular momentum is 0" You want the total linear momentum to be zero (the keep the system in the same part of the coordinate space). There is no reason a priori to expect the total angular momentum to vanish. Hopefully you have just mistyped. $\endgroup$ – dmckee --- ex-moderator kitten Aug 26 '13 at 2:26
5
$\begingroup$

Your problems is that variants Euler is the incorrect method for gravitational systems. You need to look into Symplectic Integrators. This will solve the stability issue you are having since this type of integrators hold certain values constant (in this case energy).

Also, as a comment on your code, you should use temporary variables for holding derivates before you update state variables. As you have it now, updating variables one by one, mixes information from two time steps. Right now you are updating the velocities explicitly and the positions implicitly. Unless you have a good reason to do so, it is better to stick to a single method for a problem so you have known stability characteristics.

Edit: As the comments stated, you were using a symplectic integrator, so that is not the issue. You simply have a couple small typos in your code.

First, in the second term of your update for v2, you used x2 and y2 rather than (x2-xs) and (y2-ys). This shouldn't make much difference though since your sun should stay near 0,0.

Your main errors was that in the update for the suns velocity, you multiplied by the product of radii rather than dividing by it. The four lines are fixed here:

    vx2 = vx2 - G*ms*dt*(x2-xs)/(r2s*r2s*r2s) - G*m1*(x2-x1)*dt/(r12*r12*r12);
    vy2 = vy2 - G*ms*dt*(y2-ys)/(r2s*r2s*r2s) - G*m1*(y2-y1)*dt/(r12*r12*r12);

    vxs = vxs - G*m1*(xs-x1)*dt/(r1s*r1s*r1s) - G*m2*(xs-x2)*dt/(r2s*r2s*r2s);
    vys = vys - G*m1*(ys-y1)*dt/(r1s*r1s*r1s) - G*m2*(ys-y2)*dt/(r2s*r2s*r2s);
Improve this answer
$\endgroup$
  • $\begingroup$ Hi! The good reason too use the new velocities to update the old position is that I'm using the Euler-Cromer method, that is symplectic! ; ) $\endgroup$ – Pincopallino Aug 5 '13 at 16:29
  • $\begingroup$ Sorry I blanked on that one. I've converted your code to matlab and will take a closer look into what is going on. $\endgroup$ – Godric Seer Aug 5 '13 at 16:40
  • $\begingroup$ No problem! Thank you very much for your help! $\endgroup$ – Pincopallino Aug 5 '13 at 16:44
  • $\begingroup$ I think ive found your error, please look at the updated answer. $\endgroup$ – Godric Seer Aug 5 '13 at 16:48
  • 1
    $\begingroup$ OMG! So stupid! I was looking at the code for too much time (I should have done a break)... Thank you very much! ; ). $\endgroup$ – Pincopallino Aug 5 '13 at 17:10
-2
$\begingroup$

double G(1.9838e-29); ??

maybe use: G = 6.673e-11

and try something lake this: q = q*q*q; // one

f = 1/q; g = 2/q; h = 3/q; // two, three, ...

Improve this answer
$\endgroup$
  • $\begingroup$ The value of Newton's gravitational constant depends on the choice of units and the OP is clearly not working is SI base units, but in a system where the length unit is AU. $\endgroup$ – dmckee --- ex-moderator kitten Dec 3 '13 at 22:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.