Wie kann man die steife Gleichung in diesem Restricted Three Body Problem numerisch lösen?

Ich bin bei der Lösung des Circular Restricted Three Body-Problems auf eine steife Gleichung gestoßen. [Ein Objekt bewegt sich unter Berücksichtigung des Effekts der Gravitationskräfte, die durch zwei in einem 2D-Raum fixierte Gravitationsquellen verursacht werden.]

Die Gleichungen sind folgende:

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

Weder die Euler-Methode noch Runge Kutta funktionieren, da die Eigenschaft in der Nähe von $(x_1, 0)$ oder $(x_1, 0)$ nicht gut ist. Die Derivate ändern sich zu schnell. Die Simulation kann nicht richtig gelöst werden. Das Objekt ist zu leicht auf die Gravitationsquelle zu treffen.

Wie kann ich das beheben?

Vielen Dank!

Herkömmliche Integratoren behalten die "Form" des Phasenraums nicht bei, was zu einem systematischen Energiegewinn oder -verlust führt. Daher sollten Sie "sympelektische" Integratoren in Betracht ziehen. Für enge Begegnungen sollten Sie die Methode von Chambers (1999) in Betracht ziehen. Ein hybrider symplektischer Integrator, der enge Begegnungen zwischen massiven Körpern ermöglicht .

Eine gut lesbare Einführung in moderne ODE-Löser finden Sie in Kapitel 7 von Cleve Molers Buch Numerical Computing With MATLAB, das online hier verfügbar ist: http://www.mathworks.com/moler/odes.pdf Zu den Themen, die er behandelt, gehört Stabilität , wie man eine vorgeschriebene Genauigkeit mit variablen Zeitschrittalgorithmen und Anwendung auf das Zweikörperproblem erhält.