Wie kann das Zweikörperproblem im ECI-Frame durch numerische Integration gelöst werden?

Dazu gibt es mehrere Möglichkeiten. Am einfachsten und unkompliziertesten ist es, es in zwei Sätze aufzuteilen, indem Sie die Geschwindigkeit als Variable einbeziehen und gemeinsam lösen.

Statt einer einzigen Differentialgleichung zweiter Ordnung

$$\ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r}$$

Wir können das folgende Paar von Differentialgleichungen erster Ordnung parallel lösen

$$\dot{\mathbf{v}} = -\frac{\mu}{r^3}\mathbf{r}$$ $$\dot{\mathbf{r}} = \mathbf{v}$$

unter Verwendung verschiedener einfacher Methoden, einschließlich Standardbibliotheken oder selbst entwickelter Implementierungen von Runge-Kutta, einschließlich meines bevorzugten einfachen RK4/5 mit variabler Schrittgröße .

Es ist wunderbar und wirklich lehrreich, das einmal für sich selbst zu codieren und die Aufgabe zuerst zu schätzen, bevor Sie von da an Standardbibliotheken verwenden.

Weitere Informationen zu Fehlern bei der Verwendung von Differentialgleichungslösern und zu einer Möglichkeit, sie selbst zu testen, finden Sie in meiner Frage in Math SE: Notwendigkeit, eine bessere Genauigkeit der ODE-Lösung im Vergleich zur numerischen Genauigkeit zu verstehen (was ich in Computational Science SE hätte fragen sollen )

Für Programmierzwecke können Sie schreiben$\mathbf{r} / r^3$alsr * (r**2).sum()**-1.5

Aus dieser Antwort können Sie eine 2D-Implementierung entnehmen, die nicht nur den Monopol verwendet$1/r^2$Gravitationsterm, sondern der zusätzliche Quadrupol$J_2$Bezeichnung für die abgeflachte Form und das Feld der Erde. Weitere Informationen dazu finden Sie in dieser Antwort auf Probleme beim Ableiten rechteckiger Komponenten der Beschleunigung eines Satelliten im Orbit um die Erde unter Berücksichtigung von J2.

Siehe auch eine ähnliche Lösung in dieser Antwort auf Wie wird die Umlaufbahn eines Satelliten für eine Kollisionserkennung bestimmt? .

Sie können auch in Betracht ziehen, einheitenlos zu werden, indem Sie verwenden$\mu = 1$und Periode$T = 2 \pi$und große Halbachse$a=1$. Einige Integratoren werden das etwas besser handhaben.

Die blaue Linie ist korrekt, die rote Linie hat J2 das Zehnfache seines tatsächlichen Wertes, um zum Spaß eine übertriebene Apsidenpräzession zu zeigen.

def deriv(X, t):

    x, v = X.reshape(2, -1)

    acc0  = -GMe * x * ((x**2).sum())**-1.5
    acc2  = -1.5 * GMe * J2 * Re**2 * x * ((x**2).sum())**-2.5

    return np.hstack([v, acc0 + acc2])


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

# David Hammen's nice table https://physics.stackexchange.com/a/141981/83380
# See http://www.iag-aig.org/attach/e354a3264d1e420ea0a9920fe762f2a0/51-groten.pdf
# https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere

GMe = 3.98600418E+14  # m^3 s^-2
J2e = 1.08262545E-03  # unitless
Re  = 6378136.3 # meters

X0 = np.hstack([6778000.0, 0.0, 0.0, 10000.])  # x, y, vx, vy

time = np.arange(0, 300001, 100)

J2 = J2e  # correct J2
answerJ2, info = ODEint(deriv, X0, time, full_output=True)

J2 = 10*J2e # 10x larger J2
answer10xJ2, info = ODEint(deriv, X0, time, full_output=True)

if True:
    plt.figure()
    x, y = answerJ2.T[:2]
    plt.plot(x, y, '-b')
    x, y = answer10xJ2.T[:2]
    plt.plot(x, y, '-r')
    plt.plot([0], [0], 'or')
    plt.gca().set_aspect('equal')
    plt.show()