Restricted 3 body simulation not working correctly

Question

I am 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. The information i have been given on the problem is: Check out this link and here is my program so far:

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

My program is working to certain extent but i am getting some weird problems that i hope someone can help me with.

The main issue is that when the test mass gets to a point in its trajectory when it should go off and start to orbit about the other mass it instead just shoots off on a straight line to infinity! at first i thought it was that the masses were colliding so i put in the radius check, but in some cases this does work, in some cases it doesn't, and in some cases the masses collide earlier on before the trajectory goes wrong anyway so this clearly isn't the issue. I am not sure if i have explained that all too well so here is a picture to show you what i mean. (the simulation on the right is from here)

However, this is not always the case, sometimes instead of going in a straight line, the trajectory just goes crazy when it should go over to the other mass, like this:

I really have absolutely no idea whats going on i have spent days trying to figure this out but just cant seem to get anywhere, so any help in identifying where my problem is would be very much appreciated.

Answer 1

This is just too long to fit in a comment and I also might be of use to future visitors.

The proper choice of timestep for a given computation is not an easy task. The family of Verlet integrators are symplectic, which means that they preserve the phase space volume and hence should preserve the total energy of the system given infinite precision and an infinitely small timestep, but unfortunately real computers operate with finite precision and the human life is too short in order for us to wait for an infinite number of timesteps.

The Verlet integrator, like the one that you have implemented and the velocity Verlet scheme, have global error which is O(Δt²). It means that the algorithm is only quadratically sensitive to the timestep and in order to greatly improve the precision, one has to decrease the timestep accordingly by as many times as the square root of the desired precision improvement factor. Click on the into button of the Flash applet that you compare your trajectories with and you'll see that it uses a completely different integrator - the Euler-Cromer algorithm (also known as semi-implicit Euler method). It has different precision given the same timestep (actually it is worse than that of the Verlet scheme given the same timestep) and hence you cannot and should not directly compare both trajectories, rather only their statistical properties (e.g. mean energy, mean velocity, etc.)

My point was that you have to decrease the timestep, because it is too large to handle the cases when the test body comes too close to one of the gravitational centres. There is another problem hidden here and it is the finite numerical precision. Observe this term:

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

Whenever you subtract two closely valued floating point numbers (x[n] and 1.0), a very unfortunate event happens, known as precision loss as most of the higher significant bits in their mantissas cancel each other and in the end, after the normalisation step, you get a number with much less significant bits than the original two numbers. This precision loss gets even bigger as the result is then squared and used as a denominator. Note that this mostly happens near the axis of symmetry of the system where y[n] comes close to 0. Otherwise y[n] might be big enough so that dxp*dxp is only a tiny correction to the value of y[n]*y[n]. The net result is that the force would come out totally wrong near each fixed mass and would usually be greater in magnitude than the actual force. This is prevented in your case as you test for the point being outside the prescribed radius.

Greater forces lead to greater displacements given a fixed timestep. This also leads to an artificial increase in the total energy of the system, i.e. the test mass would tend to move faster than in a finer simulation. It also might happen that the test body ends up so close to the gravitational centre, that the huge force times the square of the timestep might give so huge a displacement, that your test mass would end up much far away, but this time with the increased total energy it would result in high kinetic energy and the body would practically be ejected from the simulation volume. This could also happen even if you compute the force with infinite precision - simply the displacement between two timesteps might be so large (because of the large timestep) that the system would make an unrealistic jump in the phase space to a completely different energy isosurface. And with gravity (as well as with electrostatics) it is so easy to get to such a case as the force increases as 1/r^2 and near the radius it is many orders of magnitude stronger than in the initial state.

One might come up with different rules of a thumb to estimate the size of the timestep given the largest expected force value, but in general the higher the maximum force, the lower the timestep should be. These kind of things can usually be roughly estimated given the initial conditions, which saves a lot of failed simulations due to "ejection" effects.

Now as the Verlet schemes are symplectic, the best way to control the correctness of the simulation is to observe the total energy of the system. Note that the velocity Verlet integrator is a bit better as it is numerically more stable (but still it has the same dependence of the accuracy on the square of the timestep). With the standard Verlet scheme you can get an approximation of the speed v[i] by taking (x[i+1] - x[i-1])/(2*dt). With the velocity Verlet, the speed is included explicitly in the equations.

Either way, it would make sense to take the speed and to compute the total energy of the system at each timestep and to observe the value. If the timestep is Just Right (tm), then the total should be almost conserved with relatively small oscillations around the mean value. If it goes crazy upwards, then your timestep is too big and should be decreased.

Decreasing the timestep increases the run-time of the simulation accordingly. One could also observe that in the far field the forces are small, the point moves slowly and a long timestep is just fine. Shorter timestep would not improve the solution there but would only increase the run-time. That's why people have invented the multi-timestep algorithms and also the adaptive timestep algorithms that automatically refine the solutions in the near field. Also a different method to compute the forces might be applied there by transforming the equations so as to not include subtraction of closely valued variables.

(well, this came out way larger than even several comments)

Answer 2

I found it difficult to understand your code, so I will try to be as helpful as I can and apologize, if I am telling you things you already know.

The best way I have found to calculate the physics for simulations involving gravity is to use Newtons Law of Universal Gravitation. This is given by the formula:

F = ( ( -G * M1 * M2 ) / ( R * R ) ) * r_unit_vector;

Where:

G ~= 6.67e-11, M1 is the mass of the first object, M2 is the mass of the second object, R is the distance between the two objects: sqrt(pow(X2 - X1, 2) + pow(Y2 - Y1, 2))

X1 is the X-coordinate of object 1, X2 is the X-coordinate of object 2, Y1 is the Y-coordinate of object 1, Y2 is the Y-coordinate of object 2.

r_unit_vector is a unit vector pointing from object 2 to object 1

struct r_unit_vector_struct{
    double x, y;
}r_unit_vector;

r_unit_vector has a x component which is object 2's x-coordinate - object 1's x-coordinate, r_unit_vector has a y component which is object 2's y-coordinate - object 1's y-coordinate.

To make r_unit_vector a unit vector, you must divide (both the x aand y components separately) by its length, which is given by sqrt(pow(r_unit_vector.x, 2) + pow(r_unit_vector.y - Y1, 2))

And you should all be ready to go! Hopefully this makes sense. If not, I will write you a class to do this stuff or explain it further if I can!