¿Cómo resolver la ecuación rígida en este problema restringido de tres cuerpos numéricamente?

Me he encontrado con una ecuación rígida para resolver el problema circular restringido de tres cuerpos. [Un objeto se mueve considerando el efecto de las fuerzas gravitacionales causadas por dos fuentes gravitacionales fijadas en un espacio 2D].

Las ecuaciones son estas:

$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}$

Ni el Método Euler ni Runge Kutta funcionarán ya que la propiedad cerca de $(x_1, 0)$ o $(x_1, 0)$ no es buena. Los derivados cambian demasiado rápido. La simulación no se puede resolver correctamente. El objeto es demasiado fácil de alcanzar en la fuente gravitacional.

¿Cómo puedo arreglar esto?

¡Gracias!

Los integradores convencionales no conservan la "forma" del espacio de fases, lo que lleva a una pérdida o ganancia de energía sistemática, por lo tanto, debe considerar integradores "simécticos". Para encuentros cercanos, debe considerar el método de Chambers (1999) Un integrador simpléctico híbrido que permite encuentros cercanos entre cuerpos masivos .

Puede encontrar una introducción muy legible a los solucionadores modernos de EDO en el capítulo 7 del libro de Cleve Moler, Numerical Computing With MATLAB, disponible en línea aquí: http://www.mathworks.com/moler/odes.pdf Entre los temas que analiza están la estabilidad , cómo obtener una precisión prescrita con algoritmos de paso de tiempo variable y la aplicación al problema de dos cuerpos.