c - C 中的 Velocity Verlet 实现困难(太阳漂移)

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;