Does this look accurate to you?


    That's the output of my home-made NBody solver. I'm wondering if it looks right to you all. Two things to be aware of:

    1) The particles shoot away from the singularity because when they are solved in on timestep with r approching 0, the force is enormous and they fly too far for the next timestep to catch them. I haven't figured out a way to fix this yet.

    2) It appears to go slower at the corners.... is this because they're more far away?

    Thanks a lot,

    -Jack C
  2. jcsd
  3. I also wrote a tutorial/paper/something on the subject that, let me tell you, is not anywhere close to as good as most of you could likely write. If you'd like to take a look though (suggestions?), you can see it (and the source to my program) at

    Thanks again,

    -Jack C
  4. tony873004

    tony873004 1,555
    Science Advisor
    Gold Member

    Do you have any Z axis in this? Some of the particles seem to hover a little too long, but a Z position could explain this.

    Yes, the corners contain the furthest test particles. They don't accelerate as fast at first.

    To fix the r approaching 0 problem, you have to prevent r from ever reaching 0. In real life, you can not have a distance of 0 from anything. The closest you can come is the combined radii of the objects. That's why you can lift a ball off the ground. If the distance were 0, then the gravitational force would be infinity and you couldn't lift it. The ball's distance from the Earth is the ball's radius + the Earth's radius.

    One way to prevent the test particles from coming too close is collision detection. Whatever that object in the middle is (Star, planet) has a size. When a test particle's distance (r) is equal to or less than the middle object's radius + the radius of the test particle (if any), destroy the test particle. If it has mass, transfer its mass and momentum to the parent object.

    Another way to prevent the test particles from coming too close is to give them some tangental velocity which should cause them to orbit the center object instead of crash to its surface. That would make a more interesting simulation too.
  5. No, just x and y.

    But the issue is not that the distance is evaulated at zero, it's that the distance is CLOSE to zero, giving that large force. The radii idea sounds interesting, but I was going for more of a point-mass simulation.

    Hmm I'll think about that.


    -Jack C
  6. tony873004

    tony873004 1,555
    Science Advisor
    Gold Member

    You can't evaluate it at exactly 0. You'll get a division by 0 error or overflow. But even close to 0 gives you unrealistic results. It still is a point-mass simulation even if you consider radius. Point mass simply means that in your formula for gravitational acceleration you consider it to be a point mass rather than a flattened sphere like Jupiter. Earth's shape is not perfectly spherical and the mass is not evenly distributed. This has effects on low earth orbit satellites. But the further you get from Earth, the more Earth starts behaving like a point mass. So you can use radii and still have point mass.

    The only way to correct the error if you don't give your objects tangental velocity or incorporate collision detection is to take a slower time step. But without any tangental velocity your objects will still be drawn towards the 0 radius point no matter how slow you run it. They are dropping like rocks, just like the Space Shuttle would drop like a rock if it had no tangental velocity.
  7. OK... could you explain this tangental velocity concept a little deeper?


    -Jack C
  8. tony873004

    tony873004 1,555
    Science Advisor
    Gold Member

    Consider two objects M and m. Big M is the heavier object. Draw a line between the two objects. Now draw another line perpendicular to this line, intersecting object m.

    Whatever direction object m moves, you can look at its velocity as a combination of radial velocity and tangental velocity. In the case of your objects, they all seem to start out with a velocity of 0. The center object pulls on them directly towards itself. All their velocity is contained in the radial component.

    Now consider a satellite in a circular low Earth orbit. It would need to have a velocity of about 7 km / s to maintain orbit. Since its orbit is circular, its radial velocity (towards or away from Earth is 0) and its tangental velocity 7 km/s (perpendicular to the direction of the center of the Earth)

    If you add your radial and tangental velocity with the pythagorean theorem you get your total velocity. Total velocity = sqrt (radial velocity squared + tangental velocity squared).

    When you give your objects some initial tangental velocity, it prevents them from falling directly towards the center mass. Instead, it makes a closest approach and then starts to receed away.

    In the following lines of your code:
    Code (Text):
        Particles[i].vi = 0;
        Particles[i].vj = 0;
    try some values other than 0. I don't know what units you're using, so just play around with them.

    Maybe put one of them outside of your loop that generates random positions so you can specify its position as well.

    You'll probably find some values that give you orbits rather than freefall trajectories onto your singularity.
  9. Aha... that's smart.

    But at somepoint won't I still have to use some method to deal with the singularity, because if they orbit but still get closer and closer.... eventually some particle is bound to get very close.

    How do "real" (professional) NBody codes solve this problem?


    -Jack C
    jack {[@]}
Know someone interested in this topic? Share a link to this question via email, Google+, Twitter, or Facebook