速度Verlet在C中的实现困难(太阳漂移)

Question

我试图通过速度verlet来模拟地球-太阳系统，但不知怎么的，太阳不会绕着原点(质量减少的地方)运行，而是会漂移。我花了很长时间仔细研究我的算法，但找不到漏洞。

有人能知道这是怎么回事吗？

下面是一个模拟图：http://i.imgur.com/5l8GzZS.png

#include <stdio.h>
#include <math.h>

double xearth,yearth,vxearth,vyearth;
double xsun,ysun,vxsun,vysun;
double dt=0.5;
double fxearth,fyearth;
double fxsun,fysun;
double r;
double G;
double ms, ma;
double rx,ry;
double t;


main(){
    FILE * pFile;
    int n;
    xearth= -2.569651552438753*pow(10,-2);  /* in AU */
    yearth= -1.008909556982513;
    xsun= 2.563054664344734*pow(10,-4);
    ysun= 6.897319465467234*pow(10,-3);

    vxearth= 1.690809814669721*pow(10,-2); /* in AU per day */
    vyearth= -4.950293720310762*pow(10,-4);
    vxsun= -5.788119594348977*pow(10,-6);
    vysun= 3.335986886320253*pow(10,-6);

    G=1.488*pow(10,-34); /* G in AU, t in day */
    ms=1.9884*pow(10,30); /* kg */
    ma=5.9722*pow(10,24); /* kg */
    t=0;


    pFile = fopen ("/file.txt", "w");

    rx=xearth-xsun;
    ry=yearth-ysun;

    r=sqrt((rx*rx+ry*ry));

    fxearth= -G*ms*ma*(rx)/pow(r,3);
    fyearth= -G*ms*ma*(ry)/pow(r,3);

    fxsun= -G*ms*ma*(-rx)/pow(r,3);
    fysun= -G*ms*ma*(-ry)/pow(r,3);

    vxearth=vxearth+.5*dt/ma*fxearth;
    vyearth=vyearth+.5*dt/ma*fyearth;

    vxsun=vxsun+.5*dt/ms*fxsun;
    vysun=vysun+.5*dt/ms*fysun;


    for(n=1; n<60000; n++){

        xearth=xearth+dt*vxearth;
        yearth=yearth+dt*vyearth;

        xsun=xsun+dt*vxsun;
        ysun=ysun+dt*vysun;

        rx=xearth-xsun;
        ry=yearth-ysun;

        r=sqrt((rx*rx+ry*ry));


        fxearth= -G*ms*ma*(rx)/pow(r,3);
        fyearth= -G*ms*ma*(ry)/pow(r,3);

        fxsun= -G*ms*ma*(-rx)/pow(r,3);
        fysun= -G*ms*ma*(-ry)/pow(r,3);


        vxearth=vxearth+dt/ma*fxearth;
        vyearth=vyearth+dt/ma*fyearth;

        vxsun=vxsun+dt/ms*fxsun;
        vysun=vysun+dt/ms*fysun;


        t=t+dt;


        fprintf(pFile,"%f\t %f\t %f\t %f\t %f\n",xearth,yearth,xsun,ysun,t);
    }

    fclose (pFile);

    return 0;

}

Answer 1

这是因为初始条件给了系统一个非零的净动量。你可以通过计算系统的初始平均速度并从所有物体的速度中减去它来修正这个问题：

double vxavg = (vxsun*ms + vxearth*ma) / (ms + ma);
double vyavg = (vysun*ms + vyearth*ma) / (ms + ma);

vxsun -= vxavg;
vysun -= vyavg;
vxearth -= vxavg;
vyearth -= vyavg;