我编写了以下程序来解决三体问题（太阳、地球、木星），初始化系统使得总角动量为 0。模拟的结果非常糟糕（太阳趋向于无穷大）。我看了一个多小时的代码试图找出问题所在。有人看到我的程序有什么问题吗？

非常感谢你。

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

编辑：

第一个答案的更正很好，但现在我改变了角动量的计算（它们是错误的）并且程序不再正常工作。代码如下：

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