Gravity Simulation failing at short distances

Question

I have a 2D space simulation with gravity acting on three bodies (e.g. sun and 2x planets). The simulation seems pretty stable with RK4 (was unstable with Euler). The problem arises if the orbit of one of the planets get VERY close to the sun (or other planet), it is liable to shoot off towards alpha-centuri.

I think the problem is due to the asymmetric nature of numerical integration in this situation; if a planet is VERY close to a sun in one time step, it will experience huge acceleration and end up (next update) on the far side of the sun, much further away with much weaker gravity (1/r^2), which won't counteract the huge acceleration experienced in the previous update.

Implementing RK4 had improved the situation over Euler, but not eliminated the problem. The game I am working on requires bodies to fly very close to other bodies and still behave correctly.

Any potential solutions?

Edit

One the the answerers, @Thelvyn suggested using fixed elliptical orbits, but that is not an option here as the planets can interact gravitationally. Also, yes I am using a fixed timestep.

Answer 1

You could consider breaking down a single physical frame based on the absolute magnitude of the acceleration vector for any body in the simulation. Essentially, you could map the number of digits in the acceleration magnitude to the number of iterations per each physical frame:

timeSlicesRequired = ceil(log10(maximumMeasuredAcceleration))
adjustedTimeStep = 1 / timeSlicesRequired
for (i = 0; i < timeSlicesRequired; i++) {
    advanceSimulationBy(adjustedTimeStep)
}

That way, you're still keeping real time, but calculating more smaller steps inside a single physical frame, hopefully mitigating most of your problems (but not all, extremely close orbits will still cause problems due to precision).

Answer 2

For stuff like this I would use elliptic orbits, this guaranties stability by definition.

Due to limited precision of floating point numbers you will reach a point where you will lose stability I guess. This often no big deal in physics engines because you have energy loss at each bounce so you're convergent.

Do you really need integration here ?

Edit: One last remark, if you use integration ensure that you use a fixed time step.