Solving N-Body Simulation Problems with Gravitational Equations

[Lijnk]

I'm having problems with an n-body simulator. I tried to use the gravitational equation for vectors, however that didn't exactly work, so I split up the forces into their respective X and Y components and used two separate gravitational equations. The problem is the entire body is exploding rather than coming in. I'm obviously doing something wrong with the gravity, I'm just not sure what it is.

void forceInit() //Initializes/Updates the forces acting on each particle from the other particles
{
  int i, j;
  double rx=0, ry=0, r1, r2;
  for(i=0;i<N;i++) //N is the number of particles
  {
    for(j=0;j<N;j++)
    {
      if(j==i) continue; //so we don't check the particle against itself
      rx += pow(part[0][j] - part[0][i], 2);
      ry += pow(part[1][j] - part[1][i], 2);
    }
    r1 = pow(rx, -1); //converts to 1/r since all particle masses are =1
    r2 = pow(ry, -1);
    part[2][i]=G*r1; //G is gravitational constant
    part[3][i]=G*r2;
    rx = 0;
    ry = 0;
  }
}

void forceNext() //Converts the force into the next location for the particle
{
  int i;
  forceInit();
  for(i=0;i<N;i++)
  {
    part[4][i] += part[2][i] * T; //T is the change in time
    part[5][i] += part[3][i] * T;
    part[0][i] += part[4][i] * T;
    part[1][i] += part[5][i] * T;
  }
}

/*part[x][y] is an array listing all the particle's properties
x is the property type and y is the particle number

property types:
0 - X coordinate
1 - Y coordinate
2 - X force
3 - Y force
4 - X Velocity
5 - Y Velocity
6 - Theta (used for initial velocity and position)
*/

I can answer any questions about the code or make it more verbose if needed.

The equation I'm basing it off of is two of these. One for the X value and the other for the Y.

(ignore the last right bracket looking thing. The original equation was the vector one, but the one I'm using doesn't need the unit vector.)

 

[Lijnk]

I've narrowed the problem down to the fact that the forces are always adding and never subtracting. I'm not completely sure why that is. I'm going to take the question elsewhere since this seems to be more programming rather than physics at the moment (I've fixed an error that was physics related in the code though).

 

[Born2bwire]

You still have some problems here. The force on a particle is

\mathbf{F}_i = -G m_i \sum_{j=0, \ j\neq i}^{N-1} \frac{m_j}{\left| \mathbf{r}_i-\mathbf{r}_j \right|^2} \hat{r}_{ij}

You can't split it up r into discrete x and y components as you have done so. Because here, \hat{r}_{ij}, is the directional vector that points from your particle i to the jth particle. Your x and y components are taken from the directional vector \hat{r}_{ij}. So for each term you need:

A = -Gm_i

B_j = \frac{m_j}{\left| \mathbf{r}_i-\mathbf{r}_j \right|^2}

C_x = \cos \phi_{ij}

C_y = \sin \phi_{ij}

where \phi_{ij} is the angle between the location of the ith particle and the jth particle. So the x component of the force is

F_{ix} = A \sum B_jC_x

and likewise the y component is

F_{iy} = A \sum B_jC_y

Again, r_{ij} is the distance between the ith and jth particles. The direction of this force is directed from the ith particle to the jth particle (since it is an attractive force). You need to find this direction and decompose it into the x and y components and scale it by the force between the two particles. Of course another way of doing it, if you do not wish to convert to cylinderical coordinates you note that

-\hat{r}_{ij} = \frac{\mathbf{r}_j-\mathbf{r}_i}{\left| \mathbf{r}_i-\mathbf{r}_j \right|}

Well it appears that I may have been a bit loose on the direction because I forgot about the minus sign in front of the equation but that can easily be cleared up by a bit of thought before implementing.

You also seem to have assumed that all the particles have unit mass. It would be trivial and more correct for you to account for their mass. In addition, you may want to take a look into a better integrator. You are using what is known as an Euler integrator. A Verlet integrator is only trivially more complex but it allows for an excellent error bound and is much more stable.

 

[Lijnk]

I did some experimentation while I was away and came to the exact conclusion. I fixed it up. It's a cluster of a function now .

void forceInit() {
  int i, j; 
  double rx=0, ry=0, r1=0, r2=0;
  for(i=0;i<N;i++) {
    for(j=0;j<N;j++) {
      if(j==i) continue;
      r1 = pow(pow(part[0][i] - part[0][j], 2) + pow(part[1][i] - part[1][j], 2), -1);
      r2 = ((part[0][i] - part[0][j])*pow(sqrt(pow(part[0][i], 2)+pow(part[0][j], 2)), -1)) - ((part[1][i] - part[1][j])*pow(sqrt(pow(part[1][i], 2)+pow(part[1][j], 2)), -1));
      rx += r1*r2;
      r1 = pow(pow(part[0][j] - part[0][i], 2) + pow(part[1][j] - part[1][i], 2), -1);
      r2 = ((part[0][j] - part[0][i])*pow(sqrt(pow(part[0][i], 2)+pow(part[0][j], 2)), -1)) - ((part[1][j] - part[1][i])*pow(sqrt(pow(part[1][i], 2)+pow(part[1][j], 2)), -1));
      ry += r1*r2;
    }
    part[2][i]=-G*rx;
    part[3][i]=-G*ry;
    rx = 0;
    ry = 0;
  }
}

As for the mass, This was basically just a debugging stage to make sure the gravity works properly (it's much easier to test the program parts at a time rather than make the entire program and debug everything at once). Mass is going to be implemented in the future as well as more efficient ways to calculate.

 

[K^2]

You REALLY shouldn't use pow() function in this program. These are done with a call to logarithm and a call to exponent.

Your call to pow(a,2) is computed as exp(2*ln(a)), when all you need is a*a. Similarly, pow(a,-1) is done as exp(-ln(a)), when all you need is 1/a. Your program will run roughly 100 times faster after you fix that.

 

[Born2bwire]

You REALLY shouldn't use pow() function in this program. These are done with a call to logarithm and a call to exponent.

Your call to pow(a,2) is computed as exp(2*ln(a)), when all you need is a*a. Similarly, pow(a,-1) is done as exp(-ln(a)), when all you need is 1/a. Your program will run roughly 100 times faster after you fix that.

Oh yeah wow. You've got a pow embedded inside another pow. If you have any integer powers or half-integer powers you should just write out the multiplication, division and/or square root explicitly.

In addition, your coefficient r1 should be the same for both components (as I showed in the equations above). In fact, I think it is in your calculations since you are squaring the terms. No need to calculate it twice.

This kind of stuff may seem trivial, but you are doing a N^2 calculation here with your nested for loops. These kind of things can quickly stack up in a code that is order N^2 or N^3.

 

[Lijnk]

The reason I put pow(a, -1) is because 1/a did some weird things previously (although it's most likely just me). I'll be fixing a lot of things up with that function such as optimization as well as the Y component still needs work done on my part as well.

I'll probably end up redoing the entire function because of how broken it is.

 

[Born2bwire]

Well, another point about 1/x. Divide takes longer than a multiply. So if you have to divide by a factor multiple times then it is generally faster to first calculate the inverse and then do a multiplication. IE:

A = 1.0/x;
B *= A;
C *= 2.0*A;
...

However, I seem to recall that this can be ever so slightly less accurate than doing the explicit division. I do a lot of numerical codes though and I have not come to experience any problems with that. So I do not think that was your problem but it can be a problem from what I have read.

 

[K^2]

It depends on relative magnitudes of factors going into numerator and denominator. But you have to get creative to it really matter in a simulation like this.

Of course, the difference between division and multiplication in floating point is nowhere near as severe as between multiplication and exponentiation, so that's not exactly an optimization I'd worry about until you need to push the sim for all it got.