In nbody simulations, often gravitational softening is included. This is a small value added to the distance in the denominator of the force computation: F = G*m_j*m_2 normalize(r_j - r) * 1/(|r_j -r|^2+epsilon)
This rounds off the singularity at |r_j - r| = 0, and is in fact also physically meaningful since in practice bodies are not point-like particles anyway, so in the real world there is never a singularity.
Numerically, this has two advantages:
When two particles have a close encounter, avoid forming bound pairs orbiting each other at extremely close distances and with velocities which may not be physically meaningful; this is particularly important when adaptive time stepping is involved (which we don't do) since this phenomenon can force time steps to be extremely small
Avoid energy errors exploding in a close encounter in case of fixed time stepping as in our case
In nbody simulations, often gravitational softening is included. This is a small value added to the distance in the denominator of the force computation:
F = G*m_j*m_2 normalize(r_j - r) * 1/(|r_j -r|^2+epsilon)This rounds off the singularity at |r_j - r| = 0, and is in fact also physically meaningful since in practice bodies are not point-like particles anyway, so in the real world there is never a singularity.
Numerically, this has two advantages: