1

我目前正在尝试做一个任务，我必须为受限的 3 体引力问题编写一个模拟，具有两个固定质量和一个测试质量。我已经获得了为此需要使用的方程式，但要么我没有正确理解它们，要么我没有正确实施。如果有人可以帮助我朝着正确的方向前进，我将不胜感激。

我得到的信息如下：

http://www.flickr.com/photos/91029993@N07/8269806430/in/photostream

基本上我一直在尝试测试单个质量，但是我的程序为测试质量的运动提供了一条直线，在很短的时间内，它看起来好像程序运行正常，但后来它不起作用你变得更高。（见图片，这些输入为 1 0.4 0.5 0 -1）

http://www.flickr.com/photos/91029993@N07/8268764897/in/photostream

我最初非常忠实地按照以下给出的方程式编写了我的程序：

#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[20000],ay[20000],mneg,mpos,time;
    int n;
    FILE* output=fopen("proj1.out", "w");

    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<150;n++)
    {
        ax[n]= (-mneg*(x[n]+1)/(pow((sqrt(pow((x[n]+1),2))),3))) -(mpos*(x[n]-1)/(pow((sqrt(pow((x[n]-1),2))),3)));
        ay[n]= (-mneg*(y[n])/(pow(y[n],3))) -(mpos*(y[n])/(pow(y[n],3)));



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


        fprintf(output, "%lf %lf\n",x[n-1], y[n-1]); 
    }
return(0);
}

然后我按照这里给出的建议尝试了一种新方法：模拟恒星的引力？，所以我现在有：

#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[20000],ay[20000],mneg,time,r,a;
    int n;
    FILE* output=fopen("proj1.out", "w");

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

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



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

        r=sqrt(pow((x[n]+1),2)+pow(y[n],2));
        a=mneg/(r*r);



        ax[n]=a*((x[n]+1)/r);
        ay[n]=a*((y[n])/r);


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


        fprintf(output, "%lf %lf\n",x[n-1], y[n-1]); 
    }
return(0);
}

但是，这并没有给我任何更好的结果。

我真的不知道从哪里开始，我想我使用的方法是正确的，但我在实际编程中看不到任何问题，所以我真的不知道发生了什么，所以任何建议或指针都在右边方向将非常非常感谢！

4

1 回答 1

1

您没有将方程正确地转换为代码，尤其是测试质量和固定质量之间距离的表达式。计算加速度的正确方法是：

double dxm = x[n] + 1.0;
double dxp = x[n] - 1.0;
double dy = y[n];
double denom_minus_inv = pow(dxm*dxm + dy*dy, -1.5);
double denom_plus_inv = pow(dxp*dxp + dy*dy, -1.5);
ax[n] = -mneg*dxn*denom_minus_inv - mpos*dxp*denom_plus_inv;
ay[n] = -mneg*dy*denom_minus_inv - mpos*dy*denom_plus_inv;

请使用临时变量来存储中间表达式 - 不要将它们全部放在一个复杂的表达式中。现代编译器在优化代码时非常擅长消除多余的临时表达式。上面的代码使用了乘法通常比除法快一点的事实1.0/pow(sqrt(x), 3.0) == pow(x, -1.5)。

我建议您将 Verlet 积分器替换为Velocity Verlet。它显式地存储粒子速度，这使您可以在每个步骤中计算系统的总能量（动能和势能的总和）。它应该在整个运行过程中几乎保持不变（给出或接受舍入和离散化误差）。如果它严重偏离，那么你就知道你没有正确计算你的力。

于 2012-12-13T13:42:36.740 回答