How to solve the stiff equation in this Restricted Three Body Problem numerically?

Question

I've come across a stiff equation in solving the Circular Restricted Three Body Problem. [An object is moving considering the effect of the gravitational forces caused by two gravitational sources fixed in a 2D Space.]

The equations are these:

$x''=-\frac{GM_1 (x-x_1)}{\sqrt{(x-x_1)^2+y^2}^3}--\frac{GM_2 (x-x_2)}{\sqrt{(x-x_2)^2+y^2}^3}$

$y''=-\frac{GM_1 y}{\sqrt{(x-x_1)^2+y^2}^3}--\frac{GM_2 y}{\sqrt{(x-x_2)^2+y^2}^3}$

Neither Euler Method or Runge Kutta will work as the property near $(x_1, 0)$ or $(x_1, 0)$ is not good. The derivatives change too fast. The simulation can't be solved out right. The object is too easy to hit on the gravitational source.

How can I fix this?

Thank you!

Answer 1

Conventional integrators do not preserve the "shape" of phase space, leading to systematic energy gain or loss, thus you should consider "sympelectic" integrators. For close encounters, you should consider the method of Chambers (1999) A hybrid symplectic integrator that permits close encounters between massive bodies.