1

I'm currently trying to do an assignment where I have to write a simulation for the restricted 3 body gravitational problem, with two fixed masses and one test mass. I have been given the equations I need to use in order to do so, but either I'm not understanding them correctly or I'm not implementing correctly. I would be very grateful is someone could help to push me in the right direction.

The information I've been given is as follows:

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

basically I've been trying to test for a single mass,but my program is giving me a straight line for the movement of the test mass, over small times it looks as though the program is functioning correctly but then it doesn't work as you get higher. (see images, inputs for these were 1 0.4 0.5 0 -1)

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

I originally wrote my program quite faithfully to the equations given as follows:

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

I then tried a new method following the advice given here :Simulate the gravitational pull of a star? , so that i now have: 

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

however, this doesn't give me any better results.

I really don't know where to go from here, I think I'm using the method right but I cant see any issues in the actual programming so I really don't know whats going on, so any advice or pointers in the right direction would be very very much appreciated!

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つの複雑な式にまとめないでください。最新のコンパイラは、コードを最適化できる場合、冗長な一時式を排除するのに非常に優れています。上記のコードは、乗算が通常除算よりも少し速いという事実と、それを使用してい1.0/pow(sqrt(x), 3.0) == pow(x, -1.5)ます。

代わりに、VerletインテグレーターをvelocityVerletに置き換えることをお勧めします。粒子速度を明示的に保存するため、各ステップでシステムの総エネルギー（運動エネルギーと位置エネルギーの合計）を計算できます。実行全体を通してほぼ同じままである必要があります（丸めおよび離散化エラーを与えるか、または取ります）。それが大きく逸脱している場合は、力を正しく計算していないことがわかります。

于 2012-12-13T13:42:36.740 に答える