Was ist die analytische geschlossene Lösung des Zweikörperproblems, um seine numerischen Integrationsergebnisse zu verifizieren?

Dies ist vorerst eine ergänzende Antwort, denn obwohl wir wissen , dass eine Zwei-Körper-Umlaufbahn auf eine Ein-Körper-Umlaufbahn um ein zentrales Potential reduziert werden kann, wird dies hier etwas ablenkend sein, und ich denke, das Ergebnis für den einen Körper im zentralen Potential sieht aus Reiniger. Siehe auch Antworten zu Können die radialen Schwingungen einer Ellipsenbahn mit einem fiktiven Zentrifugalpotential gelöst werden?

Aufgrund dieses Kommentars weiß ich, dass ich irgendwo auf dieser Seite (oder in Astronomy SE ) eine Diskussion hatte, in der mir zuerst erklärt wurde, dass Kepler-Umlaufbahnen analytische Lösungen haben, die Sie für die Zeit als Funktion der Position aufschreiben können , obwohl wir müssen immer noch numerische Techniken (z. B. die Newton-Methode) verwenden, um die Position als Funktion der Zeit zu lösen. (siehe auch Wie haben Newton und Kepler das (eigentlich) gemacht? )

Wenn jemand es vor mir findet, können Sie gerne einen Link hier hinzufügen, danke!

Gleichung 27 in der Kepler-Umlaufbahn von Wikipedia; Eigenschaften der Bahngleichung ist

Wo$a$ist die große Halbachse,$\mu$ist der Standard-Gravitationsparameter, auch als Produkt bekannt$GM$,$e$ist die Exzentrizität und$E$ist die exzentrische Anomalie .

Die Beziehung zwischen$E$und die wahre Anomalie $\theta = \arctan2(y, x)$Ist

und Lösung für$E$:

Wiedereinsetzen in die erste Gleichung (aber nicht alles ausschreiben):

Versuchen wir eine numerische Überprüfung dieses erstaunlichen Ergebnisses. Beachte das mit$a=1$Und$\mu=1$die Periode ist$2 \pi$.

Das letzte Diagramm unten links zeigt, dass die analytische$t(\theta)$bezogen auf$\theta$aus einer numerisch integrierten Umlaufbahn entspricht der Zeit, die in der numerischen Berechnung für an verwendet wird$e=0.8$elliptische Umlaufbahn. Es wird numerische Störungen oder Singularitäten an den Endpunkten und für geben$e=1$aber es sieht gut aus!

Python-Skript:

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

def deriv(X, t):
    x, v = X.reshape(2, -1)
    acc = -x * ((x**2).sum())**-1.5
    return np.hstack((v, acc))

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

e = 0.8
a = 1.0
mu = 1.0
r_peri, r_apo = a*(1.-e), a*(1.+e)
v_peri, v_apo = [np.sqrt(2./r - 1./a) for r in (r_peri, r_apo)]
T = twopi * np.sqrt(a**3/mu)

X0 = np.array([r_peri, 0, 0, v_peri])
X0 = np.array([-r_apo, 0, 0, -v_apo])
times = np.linspace(-T/2., T/2., 1001)

answer, info = ODEint(deriv, X0, times, full_output=True)

x, y = answer[1:-1].T[:2]
theta = np.arctan2(y, x)
E = 2. * np.arctan(np.sqrt((1.-e)/(1.+e)) * np.tan(theta/2))
t = a * np.sqrt(a/mu) * (E - e * np.sin(E))

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    plt.plot(x, y)
    plt.plot([0], [0], 'ok')
    plt.gca().set_aspect('equal')
    plt.title('y vs. x numerical')
    plt.subplot(2, 1, 2)
    plt.plot(times[1:-1], x)
    plt.plot(times[1:-1], y)
    plt.xlim(-pi, pi)
    plt.title('x(t) and y(t) numerical')
    plt.show()

    plt.subplot(2, 2, 1)
    plt.title('theta(t_numerical)')
    plt.plot(times[1:-1], theta)
    plt.xlim(-pi, pi)
    plt.ylim(-pi, pi)
    plt.gca().set_aspect('equal')
    plt.subplot(2, 2, 2)
    plt.title('E_analytic(theta_numerical)')
    plt.plot(E, theta)
    plt.xlim(-pi, pi)
    plt.ylim(-pi, pi)
    plt.gca().set_aspect('equal')
    plt.subplot(2, 2, 3)
    plt.title('theta(t_analytic)')
    plt.plot(t, theta)
    plt.xlim(-pi, pi)
    plt.ylim(-pi, pi)
    plt.gca().set_aspect('equal')
    plt.subplot(2, 2, 4)
    plt.title('t_analytic(t_numerical)')
    plt.plot(t, times[1:-1])
    plt.xlim(-pi, pi)
    plt.ylim(-pi, pi)
    plt.gca().set_aspect('equal')
    plt.show()

Die Entfernung vom Anziehungspunkt einer Umlaufbahn kann als Funktion der wahren Anomalie (Winkel) ausgedrückt werden, die durch gegeben ist$r(\theta)=a\frac{1-e^2}{1+ecos(\theta)}$, Wo$a$ist die große Halbachse und$e$ist die Exzentrizität.