c - 受限3人体模拟无法正常工作

我现在正在做一个作业，我必须为受限的3体引力问题写一个模拟，有两个固定质量和一个测试质量。我得到的关于这个问题的信息是：Check out this link到目前为止，这是我的程序：

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

int main (int argc, char* argv[])
{
    double dt=0.005, x[20000],y[20000],xv,yv,ax,ay,mneg,mpos,time,radius=0.01;
    int n,validation=0;
    FILE* output=fopen("proj1.out", "w");

    printf("\n");


    if((argv[1]==NULL) || (argv[2]==NULL) || (argv[3]==NULL) || (argv[4]==NULL) ||     (argv[5]==NULL) || (argv[6]==NULL))
    {
        printf("************************ ERROR! ***********************\n");
        printf("**     Not enough comand line arguments input.       **\n");
        printf("**     Please run again with the correct amount (6). **\n");
        printf("*******************************************************\n");
        validation=1;
        goto VALIDATIONFAIL;
    }
if((sscanf(argv[1], "%lf", &mneg)==NULL) || (sscanf(argv[2], "%lf", &mpos)==NULL) || (sscanf(argv[3], "%lf", &x[0])==NULL) ||
(sscanf(argv[4], "%lf", &y[0])==NULL) || (sscanf(argv[5], "%lf", &xv)==NULL) || (sscanf(argv[6], "%lf", &yv)==NULL) )
{
    printf("************************* ERROR! ************************\n");
    printf("** Input values must be numbers. Please run again with **\n");
    printf("** with numerical inputs (6).                          **\n");
    printf("*********************************************************\n");
    validation=1;
  goto VALIDATIONFAIL;
}

sscanf(argv[1], "%lf", &mneg);
sscanf(argv[2], "%lf", &mpos);
sscanf(argv[3], "%lf", &x[0]);
sscanf(argv[4], "%lf", &y[0]);
sscanf(argv[5], "%lf", &xv);
sscanf(argv[6], "%lf", &yv);


    x[1]=x[0]+(xv*dt);
    y[1]=y[0]+(yv*dt);

for(n=1;n<10000;n++)
    {
        if(x[n-1]>=(1-radius) && x[n-1]<=(1+radius) && y[n-1]>=(0-radius) && y[n-1]<=(0+radius))
        {
            printf("Test mass has collided with M+ at (1,0), Exiting...\n");
            goto EXIT;
        }
        else if(x[n-1]>=(-1-radius) && x[n-1]<=(-1+radius) && y[n-1]>=(0-radius) && y[n-1]<=(0+radius))
        {
            printf("Test mass has collided with M- at (-1,0), Exiting...\n");
            goto EXIT;
        }
        else
        {
            double dxn = x[n] + 1;
            double dxp = x[n] - 1;
            double mnegdist = pow((dxn*dxn + (y[n]*y[n])), -1.5);
            double mposdist = pow((dxp*dxp + (y[n]*y[n])), -1.5);

            ax =  -(mpos*dxp*mposdist+mneg*dxn*mnegdist);
            ay =  -(mpos*y[n]*mposdist+mneg*y[n]*mnegdist);


            x[n+1]=((2*x[n])-x[n-1] +(dt*dt*ax));
            y[n+1]=((2*y[n])-y[n-1]+(dt*dt*ay));


            fprintf(output, "%lf %lf\n",x[n-1], y[n-1]);
        }
    }
VALIDATIONFAIL:
printf("\n");
return(EXIT_FAILURE);

EXIT:
return(EXIT_SUCCESS);
}

我的程序在一定程度上起了作用，但我遇到了一些奇怪的问题，我希望有人能帮助我。
主要问题是，当测试质量到达其轨道的某一点时，它应该离开并开始围绕另一个质量旋转，而不是沿着直线发射到无穷远！一开始我以为是质量碰撞，所以我做了半径检查，但在某些情况下，这确实有效，在某些情况下，它不起作用，在某些情况下，质量碰撞之前，轨道出错无论如何，所以这显然不是问题。我不确定我是否解释得太好了，所以这里有一张图片向你展示我的意思。（右边的模拟来自here）
然而，情况并非总是这样，有时不是沿着直线走，而是在轨道应该转到另一个质量时变得疯狂，比如：
我真的完全不知道发生了什么事，我花了好几天的时间想弄清楚，但似乎什么也找不到，所以如果能帮我找出问题所在，我将不胜感激。

最佳答案

这太长了，不适合发表评论，我也可能对未来的访问者有用。
正确选择给定计算的时间步长不是一件容易的事情。Verlet积分器族是辛的，这意味着它们保留相空间体积，因此在给定无穷精度和无穷小时间步长的情况下，它们应该保留系统的总能量，但不幸的是，真正的计算机以有限的精度运行，而人类的寿命太短，我们无法等待无限的时间步。
Verlet积分器，像你已经实现的和velocity-Verlet方案一样，有全局误差O（Δt2）。这意味着该算法只对时间步长二次敏感，为了大大提高精度，必须相应地将时间步长减小到所需精度提高因子平方根的数倍。点击Flash小程序的into按钮，比较你的轨迹，你会发现它使用了一个完全不同的积分器-Euler-Cromer算法（也称为半隐式Euler方法）。给定相同的时间步长（实际上比给定相同时间步长的Verlet方案的精度差），因此不能也不应该直接比较两个轨迹，而只比较它们的统计特性（如平均能量、平均速度等）
我的观点是，你必须减少时间步长，因为它太大，当测试物体离其中一个引力中心太近时，无法处理这些情况。这里还隐藏着另一个问题，那就是有限的数值精度。遵守本条款：

double dxp = x[n] - 1;
double mposdist = pow((dxp*dxp + (y[n]*y[n])), -1.5);

每当您减去两个值相近的浮点数（x[n]和1.0）时，就会发生一个非常不幸的事件，称为精度损失，因为它们尾数中的大部分高有效位互相抵消，最后，在标准化步骤之后，您会得到一个比原来两个数的有效位低得多的数字。当结果被平方并用作分母时，这种精度损失会变得更大。注意，这主要发生在系统的对称轴附近，其中y[n]接近0。否则y[n]可能足够大，因此dxp*dxp只是对y[n]*y[n]值的一个微小修正。最终的结果是，在每个固定质量附近，力会完全错误地出来，并且通常比实际力的大小更大。在您的情况下，当您测试的点不在指定的radius范围内时，这是可以避免的。
在给定固定时间步长的情况下，较大的力会导致较大的位移。这也会导致系统总能量的人为增加，即测试质量的移动速度往往比精细模拟中的快。也可能会发生这样的情况：试验体最终离引力中心太近，巨大的力乘以时间步的平方可能会产生巨大的位移，你的试验质量最终会远得多，但这次随着总能量的增加，它会产生高动能，而且实际上会被弹射出来从模拟卷。这也可能发生，即使你以无限精度计算力-简单地说，两个时间步之间的位移可能太大（因为时间步太大），系统会在相空间中不切实际地跳到一个完全不同的能量等值面。在重力作用下（以及在静电作用下），很容易得到这样一种情况，即力的增加小于cc>，在半径附近，它比初始状态强很多个数量级。
人们可能会提出不同的拇指法则来估计时间步长的大小，因为最大的期望力值，但一般来说，最大力越大，时间步长越低。在给定初始条件的情况下，这类情况通常可以粗略估计，从而节省了许多由于“喷射”效应而失败的模拟。
由于Verlet格式是辛格式，所以控制仿真正确性的最佳方法是观察系统的总能量。请注意，速度Verlet积分器在数值上更稳定，因此它更好一些（但它对时间步长平方的精度依赖性仍然相同）。用标准的Verlet格式，可以通过使用1/r^2来获得速度v[i]的近似值。对于速度Verlet，速度被显式地包含在方程中。
不管是哪种方式，都有必要计算速度，计算系统在每个时间步的总能量，并观察其值。如果时间步正好（tm），那么总的时间步应该几乎是守恒的，在平均值附近有相对较小的振荡。如果它变得疯狂向上，那么你的时间步太大，应该减少。
减少时间步长会相应地增加模拟的运行时间。我们还可以观察到，在远场中，力很小，点移动缓慢，长时间的步进就可以了。较短的时间步不会改进那里的解决方案，但只会增加运行时间。这就是为什么人们发明了多时间步算法，也发明了自适应时间步算法，可以在近场自动优化解。此外，还可以通过转换方程来应用不同的计算力的方法，从而不包括近值变量的减法。
（好吧，这比几条评论都要大得多）