N-body simulation- objects accelerate *past* each other

[RossH]

Homework Statement

I am making a standard N-body simulation for a computer science class and, for some reason, the objects that I am simulating accelerate past each other, and then engage in spring-like motion, getting closer and then farther apart, etc. Obviously incorrect.

Homework Equations

If you can't read C++, this is how I am calculating the forces (in x-y plane):

1) Calculate distance with Pythagorean theorem
2) set total force on the hypotenuse equal to Gm1m2/r^2
3) Find the distance on x and y-axis between the two bodies
4) The x component of the force is equal to the totalforce*x_distance/hypotenuse length
5) Same for y component
6) Update the x and y components by subtracting the forces found from the total force on the object at this time

Is there anything really obvious wrong with this? I'm not quite sure about the subtraction in step (6), but it seems to work. Bodies push away from each other when I have it as addition; they attract when I have it as subtraction. Anyway, I just can't figure out why mo bodies accelerate past each other and exhibit this strange motion where they refuse to stick together. Any help would be much appreciated.


If you can read C++, the code that I am using to calculate forces between two objects is this (should be pretty self-explanatory code)

for(int i=0; i<bodies.size(); i++) //body being acted upon
{
bodies.setforcex(0); //reset all forces to zero
bodies.setforcey(0);
for(int j=0; j<bodies.size(); j++) //acting body
{
if(i!=j) //update forces
{
double dist=distance(bodies.getxloc() , bodies.getyloc() , bodies[j].getxloc() , bodies[j].getyloc());
if(dist>(bodies.getradius()+bodies[j].getradius())) //basic collision detection: don't do this if the bodies are too close
//prevents asymptotic behavior of the force graph
{
double totalforce=G*bodies.getmass()*bodies[j].getmass()/(dist*dist); //force vector, now I need to split it into x, y components
double diff_x=bodies.getxloc()-bodies[j].getxloc();
double diff_y=bodies.getyloc()-bodies[j].getyloc();
//first solve for x,y components
double x_component=totalforce*(diff_x/dist); //cosine triangle times total force
double y_component=totalforce*(diff_y/dist); //sine triangle times total force

//Now add forces to objects
if(bodies[j].getneg())
{
bodies.setforcex(bodies.getforcex()+x_component);
bodies.setforcey(bodies.getforcey()+y_component);
}
else
{
bodies.setforcex(bodies.getforcex()-x_component);
bodies.setforcey(bodies.getforcey()-y_component);
}
}
}
}
}

 

[cepheid]

Yikes, this code is really hard to read. You should really enclose code in [noparse]

[/noparse] tags like below. It preserves whitespace.

for(int i=0; i<bodies.size(); i++) { // body being acted upon 

    bodies[i].setforcex(0); //reset all forces to zero
    bodies[i].setforcey(0);

    for(int j=0; j<bodies.size(); j++) { // acting body

        if(i!=j) { // update forces

            double dist=distance(bodies[i].getxloc(), bodies[i].getyloc(), 
                             bodies[j].getxloc(), bodies[j].getyloc());

            if(dist>(bodies[i].getradius()+bodies[j].getradius())) { 

             // basic collision detection: don't do this if the 
             // bodies are too close
             // prevents asymptotic behavior of the force graph

                double totalforce=G*bodies[i].getmass()
                *bodies[j].getmass()/(dist*dist); 

             // force vector, now I need to split it into x, y components
                double diff_x=bodies[i].getxloc()-bodies[j].getxloc(); 
                double diff_y=bodies[i].getyloc()-bodies[j].getyloc();

             // first solve for x,y components

             // cosine triangle times total force
                double x_component=totalforce*(diff_x/dist); 
             // sine triangle times total force 
                double y_component=totalforce*(diff_y/dist);
 
             // Now add forces to objects
                if(bodies[j].getneg()) {
            
                    bodies[i].setforcex(bodies[i].getforcex()+x_component);
                    bodies[i].setforcey(bodies[i].getforcey()+y_component);

                }   else {
                    bodies[i].setforcex(bodies[i].getforcex()-x_component);
                    bodies[i].setforcey(bodies[i].getforcey()-y_component);
                }
            }
        }
    }
}

I also used my own formatting/indentation etc., because I do not know what yours was (this information was lost when the code got posted).

 

[cepheid]

Homework Statement

I am making a standard N-body simulation for a computer science class and, for some reason, the objects that I am simulating accelerate past each other, and then engage in spring-like motion, getting closer and then farther apart, etc. Obviously incorrect.

Homework Equations

If you can't read C++, this is how I am calculating the forces (in x-y plane):

1) Calculate distance with Pythagorean theorem
2) set total force on the hypotenuse equal to Gm1m2/r^2
3) Find the distance on x and y-axis between the two bodies
4) The x component of the force is equal to the totalforce*x_distance/hypotenuse length
5) Same for y component
6) Update the x and y components by subtracting the forces found from the total force on the object at this time

Is there anything really obvious wrong with this? I'm not quite sure about the subtraction in step (6), but it seems to work. Bodies push away from each other when I have it as addition; they attract when I have it as subtraction. Anyway, I just can't figure out why mo bodies accelerate past each other and exhibit this strange motion where they refuse to stick together. Any help would be much appreciated.If you can read C++, the code that I am using to calculate forces between two objects is this (should be pretty self-explanatory code)

You asked about the step in red above. Well, it seems to me that force is a vector, and so whether or not you add or subtract x_component to/from the cumulative sum that is being stored** depends on where the two bodies are in relation to each other. If body j is to the right of body i, and if right is defined as the "positive" x-direction, then the x-component of the force on body i due to body j is positive. If body j is to the left of body i, then the opposite is true. The same goes for the y-component. Is body j above or below body i?

Maybe you are already doing this with getneg(). I don't know, because it's not clear to me what getneg() does (EDIT, but I don't think that that is what it does, because it seems to change sign for both x and y).

If taking this into account properly still doesn't sort out your problem, then I would say that maybe the problem lies with the code that takes the x and y accelerations of each body that you deduce from the forces, and numerically integrates it to get updated velocities and positions for each body. (I presume that this code must exist as part of the program, even if it is not in the excerpt that you have provided us).

**presumably the cumulative force sum that is being stored can be retrieved using getforce(). If it can't, then I don't understand how the code works.

EDIT: I just realized that x_component and y_component have intrinsic signs, because diff_x and diff_y have intrinsic signs. (The diff would be positive if body i's x_loc were greater than body j's x_loc, and negative otherwise). But it seems like you really want the diff to have the opposite sign (i.e. subtract i's position from j's, not the other way around). On the other hand, the negative sign you put in front of the force component when you add it to the cumulative sum should take care of this reversal in the sign of diff. So...I don't know.

Maybe the behaviour is correct. Any two-body system within the simulation will presumably have SOME non-zero angular momentum to it, suggesting that the bodies will not head straight towards each other and all collapse to a point, but rather "swing past" each other on various elliptical, parabolic, or hyperbolic trajectories.

EDIT 2: try simulating it with total number of bodies = 2.
- if they start at rest, do they simply accelerate towards each other?
- if they start with some non-zero tangential speed, do they "orbit" each other?
- Does the shape of this orbit (elliptical and bound vs parabolic or hyperbolic and unbound) depend on the initial speeds in a sensible way?

EDIT 3: I finally figured it out: it's your collision "veto" that is causing the problem. If the bodies are intitially at rest, they'll accelerate towards each other. But once they get within 2 object radii of each other, the collision veto kicks in, the force between them goes to 0, and fails to get updated. As a result of this (and Newton's 1st law), they coast past each other until their distance once again gets large enough for the collision veto to get turned off. Suddenly, a gravitational force between them "turns back on". They then begin to accelerate back towards each other. This repeats ad infinitum, hence the oscillatory behavoiur.

You need to add code that realistically deals with the collision once the two objects get within two object radii of each other. Maybe you could assume the collisions were elastic and give them "after collision" velocities that depend on the ratio of their masses.

EDIT 4: On the other hand, if you are trying to simulate a collisionless system, then maybe this behaviour is just fine and is what you want. Two particles can just fall "through" each other and come out the other side. Maybe they are dark matter particles or something. I think this is actually a pretty good description of something like dark matter particles in a dark matter halo, or stars in a globular cluster, or maybe even galaxies in a galaxy cluster (although the latter do tend to collide and interact with each other). The individual bodies in the system don't have enough energy to escape, so they remain bound to the system. They individually fall into the potential well on some trajectory, get accelerated, pass through the centre, and come back out the other side (eventually decelerating as they climb back out of the potential well). This sort of thing happens over and over again, but on a different trajectory each time. So the motions, rather than being ordered, are chaotic. On the other hand, dark matter does clump, in spite of being collisionless. This is something that I don't fully understand (but it is also off topic).

 

[RossH]

Haha, sorry, I didn't realize that this site had code tags.

Anyway, I've kind of come to the same conclusion as you. I don't actually do numerical integration; I just use the deltaT variable as a timestep, and then sort of do a Riemann sum (actually, I guess that is numerical integration of a sort). I think that the bodies were going past each other since the deltaT was too large, and at the next recalculation, the bodies were too far from each other to overcome their already existing velocity. And then, if I wait long enough, it goes into endless oscilations. Yeah, actually, I don't think the physics is wrong in my program. Oh, and I figured out that I had to reverse the negative and positive since I'm using the vector form of the universal gravitation equation and forgot to include the negative sign before G.

Thanks for the help!

 

[D H]

Oh, and I figured out that I had to reverse the negative and positive since I'm using the vector form of the universal gravitation equation and forgot to include the negative sign before G.

But you are doing that, right here:

             // cosine triangle times total force
                double x_component=totalforce*(diff_x/dist); 
             // sine triangle times total force 
                double y_component=totalforce*(diff_y/dist);

Adding getneg() to the mix appears to be undoing this.