NBody evaluation, is my model correct?

  • Thread starter crepincdotcom
  • Start date
  • Tags
    Model
In summary: +9 for x, y, z int xsign,ysign,zsign; //signs: +ve,+ve,+ve for x, y, z int mx,my,mz; //coordinates of center mass, leftmost particle, rightmost particle double vi,vj,vk; //velocities of center mass, leftmost particle, rightmost particle double fx
  • #1
crepincdotcom
24
0
Hello all,

I've written a very simple nbody solver in c. Basically, I just went through an array of particles and solved F = G(m1)(m2) / r ^2

The output is here: (sorry for the windows format)
http://www.crepinc.com/files/galax.wmv

As you can see, some of the particles don't do what one would expect...

So far the only possible reason I can come up for is that I'm solving the equations at a set dT. Thus if the dT is too large, "stuff happens" between the two different solves. But if dT is too small, the whole thing basically blows up.

Am I overlooking something, or is my math ok and my code is just bad?

Thanks a lot,

-Jack Carrozzo
http://www.crepinc.com/
 
Last edited by a moderator:
Astronomy news on Phys.org
  • #2
crepincdotcom said:
Am I overlooking something, or is my math ok and my code is just bad?

Can't tell without seeing the code, especially the lines that contain the time step variable.
 
  • #3
The time step can intrdouce such problems - you tend to find that energy isn't conserved. I would tend to put in energy conservation 'by hand'. However, there does seem to be another problem with your animation - it takes some time for all of the particles to start moving, and then some of them stop moving. I would guess that this is due to a problem with the code.
 
  • #4
OK, I will post the code at 8:00, sorry I can't get it before then.

As to the particles that stop, It usually doesn't do that. I think I've decreased dT too much or something... I didn't make any hardcore code changes.

But, I don't want to bother you guys with code issues. In terms of theory, does evaluating F and then making a sum of all forces do the trick? Or do I need to find a way to use the differentia l equations and such? :bugeye:

-Jack Carrozzo
http://www.crepinc.com/
 
  • #5
Hello,
Try to work with different DT. Best of Luck. Zaheer A
 
  • #6
OK, at the bottom of this post is the code. I got the model to look FAIRLY accurate by playing with the dT. But if I go higher or lower, the whole system no longer works accuratly. IS there a logical explanation for this?

The other issue is that of the Y axis. You can see in the "MoveParticles" function that when taking the sum of all y-forces and adding them to the velocity, I multiply the acceleration by 25. I have no idea why this needs to happen. All I know is that, whithout that line, the simulation seems to "favor" the X axis, with the particles hardly moving in the Y.

Finally, one should notice that I removed some 0's from the constant G. This is because multiplying something X 10^ -11 with a number (Mass) like 10 ^ 6, there are a lot of extra spaces in there. I was told it was easier and more accurate in the evaluation of that code if I simply lopped off some zero's. So, any mass you see is actually large by a factor of 10^5.

Here's the code: (it will only let me post little bits, the whole thing is at crepinc.com/files/nbody.c)

Code:
#define NUMPARTS   100 //how many particles
#define GRID_X     640
#define GRID_Y     480
//known working values for these are 10, 100000, 0.004, 2
#define MASS       10 //plus 10 zeros
#define CENTERMASS 100000 //plus 10 zeros
#define TIMESTEP   0.004 //seconds per step
#define IMAGESTEP  2 //make bmp every n steps

#define rad(d)     (3.1415926 * d) / 180 //radians from degress
#define DEBUG      0 //0 off, 1 on

struct PartType
{
  int               x; //x coord for body
  int               y; //y coord
	unsigned long int m;
  double            vi;//i-hat of velocity vector
  double            vj;//j-hat "                "
};

struct ForceType
{
  double i;
  double j;
};

struct PartType Particles[NUMPARTS];
struct ForceType ForceSum;
struct ForceType MoveForce[NUMPARTS]; //where forces are stored until move()
int ImageGrid[GRID_X][GRID_Y];
int Steps,ERR;
int ImgSteps=0;

Code:
void InitParticles( void )
{
  int i,j;

	srand(time()*getpid());
	Particles[0].x = (int)(GRID_X / 2);
  Particles[0].y = (int)(GRID_Y / 2);
  Particles[0].m = CENTERMASS;
  Particles[0].vi = 0;
  Particles[0].vj = 0;

	for (i=1;i<NUMPARTS;i++) //starts at 1, not 0
  {
    Particles[i].x = ((double)rand()/((double)pow(2,31))) * GRID_X;
    Particles[i].y = ((double)rand()/((double)pow(2,31))) * GRID_Y;
    Particles[i].m = MASS;
    Particles[i].vi = 0;
    Particles[i].vj = 0;
  }

	for (i=0;i<GRID_X;i++)
	{
		for (j=0;j<GRID_Y;j++)
		{
			ImageGrid[i][j]=0;
		}
	}
}

Code:
void ComputeForces( void )
{
  int i,j,fxsign,fysign; //signs: 0 +, 1 -
  double r,F,theta,dx,dy,fx,fy;
  double G=0.667; //has 10 zeros taken off

  for (i=0;i<NUMPARTS;i++)
  {
		ForceSum.i=0;
  	ForceSum.j=0;
    for (j=0;j<NUMPARTS;j++)
    {
      if (j!=i)
      {
				fxsign=0;
				fysign=0;
				if ((Particles[j].x - Particles[i].x) < 0) fxsign++;
				if ((Particles[j].y - Particles[i].y) < 0) fysign++;
				dx = abs(Particles[i].x - Particles[j].x);
        dy = abs(Particles[i].y - Particles[j].y);
        theta = atan(dy/dx); //theta in radians
        r = sqrt((dx * dx) + (dy * dy)); //distance
				F = ((G * (double)(Particles[i].m) * (double)(Particles[j].m)) / (r * r));
				if (DEBUG) printf("F= %e, G= %e m1= %d, m2= %d, r= %f\n",F,G,Particles[i].m,Particles[j].m,r);
				fx = cos(rad(theta)) * F;
        fy = sin(rad(theta)) * F;
				if (DEBUG) if (i==2) printf("compareing to %d, r - %f, force x %f, force y %f, theta %f\n",j,r,fx,fy,theta);
				if (fxsign==0)
					ForceSum.i += fx;
				else
					ForceSum.i -= fx;

				if (fysign==0)
        	ForceSum.j += fy;
				else
					ForceSum.j -= fy;
      }
    }
    MoveForce[i].i = ForceSum.i;
    MoveForce[i].j = ForceSum.j;
  }
}

Code:
void MoveParticles( void )
{
  int i;
  double A; //acceleration

  for (i=0;i<NUMPARTS;i++)
  {
    //F=MA, A=F/M
    A = (MoveForce[i].i / (Particles[i].m));
    Particles[i].vi += (A * TIMESTEP);
    A = (MoveForce[i].j / (Particles[i].m));
		A *= 25; //don't know why, but need this
    Particles[i].vj += (A * TIMESTEP);

		Particles[i].x += (int)(Particles[i].vi);
		Particles[i].y += (int)(Particles[i].vj);

		if (DEBUG) if (i==2) printf("part 2, y %f, f %f, %f.\n",Particles[i].vj,MoveForce[i].i,MoveForce[i].j);
  }
  Steps++;
}
("DataOut and bitmap functions left out)
Code:
int main( int argc, char *argv[] )
{
	int i, frames;

	printf("Simple NBODY solver by Jack Carrozzo.\n");

	if (argc!=2)
	{
		printf("Usage: ./nbody <frames>\n");
		exit(1);
	}

	frames = atoi(argv[1]);

	printf("Will go for %d iterations, writing %d frames.\n",(int)(frames*IMAGESTEP),frames);

	InitParticles();

	if (DEBUG)
	{
			printf("\n   [+] Time: %d [+]\n",Steps);
			for (i=0;i<NUMPARTS;i++)
			{
				printf("----------------------\n");
				printf("Particle Num     : %d\n",i);
				printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
				printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
				printf("----------------------\n");
		}
	}
	DataOut();

	for (i=0;i<(frames*IMAGESTEP);i++)
	{
		ComputeForces();
		MoveParticles();
		DataOut();
	}

	if (DEBUG)
	{
		printf("\n   [+] Time: %d [+]\n",Steps);
		for (i=0;i<NUMPARTS;i++)
		{
			printf("----------------------\n");
			printf("Particle Num     : %d\n",i);
			printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
			printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
			printf("----------------------\n");
		}
	}
	printf("Done.\n");

  return 0;
}

Thanks a lot,

-Jack Carrozzo
http://www.crepinc.com/
 
Last edited:
  • #7
OK the source above can by found fully at: crepinc.com/files/nbody.c
The movie it put out is at: crepinc.com/files/galax2.wmv (sorry for WMV)
Notice how it's better than the first movie (crepinc.com/files/galax.wmv), but how some particles still stop. I can find no reason for this. All the particles are evaluated exactly the same way, why are some treated differently?

FYI to compile the code: gcc -O3 -lm -onbody nbody.c

Thanks

-Jack Carrozzo
http://www.crepinc.com/
 
Last edited:
  • #8
Does the code work properly when TIMESTEP = 1?
 
  • #9
It DID work, though everything went so fast it was a little hard to tell.

That's a good point, I'll check it out tomorrow when I can get to my other comp.

Thanks

-Jack Carrozzo
http://www.crepinc.com/
 
  • #10
In fx = cos(rad(theta)) * F; theta is already in radians, so multiplying it by pi/180 is an error
 
  • #11
Aha! I'm an idiot :redface: , quite sorry for that. I suppose that would account for the X-favouring issue. I'm not sure why it would have worked before and not not, but tomorrow I'll fix this and see what I get.

Thanks a lot,

-Jack Carrozzo
http://www.crepinc.com/
 
  • #12
OK, I fixed the Radians problem. It seems to be working now. The only issue seems to be the singularity issue. I would imagine that this is an inherant issue in every code of this sort, am I right?

The pretty-much ok video is at: http://creponc.com/files/galax3.wmv
A funny mess up: http://crepinc.com/files/blowup.wmv

Thanks.

-Jack Carrozzo
http://www.crepinc.com/
 
Last edited by a moderator:
  • #13

FAQ: NBody evaluation, is my model correct?

1. How do you determine the accuracy of an NBody model?

The accuracy of an NBody model is typically determined by comparing its predictions to observed data. This can involve analyzing the trajectory and motion of celestial objects, such as planets or stars, and seeing if they align with the model's predictions. Additionally, running simulations with different initial conditions and comparing the results can also help assess the accuracy of the model.

2. What factors can affect the accuracy of an NBody model?

There are several factors that can affect the accuracy of an NBody model, including the complexity of the model, the precision of the calculations, and the accuracy of the initial conditions and input parameters. Additionally, factors such as external forces, like gravitational pull from other objects in the universe, can also impact the accuracy of the model.

3. How do you incorporate new data into an existing NBody model?

To incorporate new data into an existing NBody model, the model would need to be updated with the new information and the calculations would need to be re-run. This could involve adjusting the initial conditions or input parameters, depending on the nature of the new data. The updated model can then be compared to the previous version to see how the new data affects its accuracy.

4. What techniques can be used to improve the accuracy of an NBody model?

There are several techniques that can be used to improve the accuracy of an NBody model. These include using more precise calculations, incorporating more complex equations or algorithms, and using more accurate initial conditions and input parameters. Additionally, running multiple simulations with different settings and comparing the results can help identify any discrepancies and improve the overall accuracy of the model.

5. How do you account for unknown variables in an NBody model?

Unknown variables, such as the presence of undiscovered celestial objects or the effects of dark matter and energy, can be difficult to account for in an NBody model. One approach is to make educated assumptions and estimations based on existing data and theories. Another approach is to run simulations with different scenarios and see how the model responds to these unknown variables, which can help refine the accuracy of the model over time.

Similar threads

Back
Top