Einführung in die Programmierung für Studierende der Physik

(Introduction to Programming for Physicists)

Vorlesung gehalten an der J.W.Goethe-Universität in Frankfurt am Main

(Sommersemester 2026)

von Dr.phil.nat. Dr.rer.pol. Matthias Hanauske

Frankfurt am Main 26.02.2026

Das Dreikörperproblem

Die Bewegung eines Körpers in einem gravitativem Zentralkraftfeld ist ein oft behandeltes System in der klassischen Mechanik (siehe z.B. Walter Greiner, 'Mechanik Teil 1' [5. Auflage, 1989, Kapitel 26. Seite 266]). Die Bewegungsgleichung eines Planeten (z.B. Erde oder Venus) im Gravitationsfeld der Sonne wird in vielen Lehrbüchern analytisch behandelt, indem man den Planeten als einen Massenpunkt mit der Masse $m$ behandelt. Die Bewegungsgleichungen des betrachteten Systems formuliert mittels der Newtonschen Gravitationstheorie lauten dann:

\begin{equation} m \cdot \frac{d^2 \vec{r}(t)}{dt^2} = \frac{G \cdot m \cdot M_\odot}{r^3} \cdot \vec{r}(t) \quad \Rightarrow \quad \frac{d^2 \vec{r}(t)}{dt^2} = \frac{G \cdot M_\odot}{r^3} \cdot \vec{r}(t) \quad , \end{equation}

wobei $\vec{r}(t)$ der Ortsvektor des betrachteten Planeten in Bezug auf unser Zentralgestirn, die Sonne (Masse $M_\odot$) und $G$ die Gravitationskonstante ist. In drei Dimensionen ($\vec{r} = ( x,y,z )$) stellt die Bewegungsgleichung ein System von drei Differenzialgleichungen zweiter Ordnung dar, dessen Lösungen man analytisch angeben kann (siehe z.B. Walter Greiner, 'Mechanik Teil 1' [5. Auflage, 1989, Kapitel 26. Seite 266]). Betrachtet man jedoch drei Himmelskörper (z.B. Sonne, Erde und Mond), die sich gegenseitig durch ihre Schwerkraft beeinflussen, so ist das entstehende System von Bewegungsgleichungen im Allgemeinen nicht mehr analytisch lösbar. Da es keine allgemeine Lösung des Dreikörperproblems gibt, haben Forscherinnen und Forscher schon früh begonnen, Spezialfälle zu untersuchen. Zum Beispiel, wenn eines der drei Objekte deutlich schwerer ist als die anderen beiden – wie es beim System Erde-Mond-Sonne der Fall ist. Des Weiteren existieren für bestimmte Anfangsbedingungen spezielle Lösungen (z.B. die Lagrange-Dreieckslösungen, siehe später).

Das Dreikörperproblem beschreibt die mathematische und physikalische Herausforderung, die Bewegungen von drei Himmelskörpern (z. B. drei Sterne oder Sonne, Erde und Mond) vorherzusagen, die sich gegenseitig durch ihre Schwerkraft beeinflussen.

Das Thema wurde vor einigen Jahren sogar in der Netflix-Serie "3 Body Problem" thematisiert, welche auf dem berühmtem Bestseller Die drei Sonnen des chinesischen Sciencefiction-Autor Liu Cixin basiert. In dieser Trisolaris-Trilogie wird eine außerirdische Zivilisation beschrieben, die auf einem Planeten lebt, der drei Sonnen umkreist. Der Bahnverlauf der Sonnen lässt sich nicht genau vorhersagen, wodurch die Aliens immer wieder vernichtet werden, wenn eine der Sonnen zu nah an ihren Planeten heranreicht.

Die zugrundeliegende Bewegungsgleichung wird durch das folgende System von gekoppelten Differenzialgleichungen zweiter Ordnung beschrieben:

\begin{eqnarray} \frac{d^2 \vec{r_1}}{dt^2} &=& G \cdot m_2 \cdot \frac{\vec{r_2}-\vec{r_1}}{\left| \, \vec{r_2}-\vec{r_1} \right|^3} + G \cdot m_3 \cdot \frac{\vec{r_3}-\vec{r_1}}{\left| \, \vec{r_3}-\vec{r_1} \right|^3} \nonumber \\ \frac{d^2 \vec{r_2}}{dt^2} &=& G \cdot m_1 \cdot \frac{\vec{r_1}-\vec{r_2}}{\left| \, \vec{r_1}-\vec{r_2} \right|^3} + G \cdot m_3 \cdot \frac{\vec{r_3}-\vec{r_2}}{\left| \, \vec{r_3}-\vec{r_2} \right|^3} \nonumber \\ \frac{d^2 \vec{r_3}}{dt^2} &=& G \cdot m_1 \cdot \frac{\vec{r_1}-\vec{r_3}}{\left| \, \vec{r_1}-\vec{r_3} \right|^3} + G \cdot m_2 \cdot \frac{\vec{r_2}-\vec{r_3}}{\left| \, \vec{r_2}-\vec{r_3} \right|^3} \quad , \nonumber \end{eqnarray}

wobei $\vec{r}_i$ und $m_i, i \in {1,2,3}$ der Ortsvektor und die Masse des Körpers $i$ und $G$ ist.

Betrachtet man die Bewegung der drei Körper in allen drei Raumdimensionen ($\vec{r_1}=\left( x_{1}, y_{1}, z_{1} \right)$), so hat gelangt man zu neun gekoppelten Differenzialgleichungen zweiter Ordnung, die man dann in ein System von 18 gekoppelten DGLs erster Ordnung umschreiben muss um sie numerisch lösen zu können.

In diesem Jupyter Notebook werden wir das Dreikörperproblem numerisch lösen. Zunächst befassen wir uns mit einigen bekannten regulären Bewegungen der drei Körper und betrachten im Speziellen die elliptischen Lagrangeschen homografischen Lösungen und die Achterschleife. Während die erst im Jahre 2000 gefundene Lösung der Figur 8 stabil ist, sind die elliptischen Lagrangeschen Lösungen, bei gleichen Massenwerten, instabil und zeigen nach einiger Zeit chaotisches Verhalten. Danach betrachten wir weitere, erst in den letzten Jahren gefundenen, reguläre planare und auch dreidimensionale Bewegungen. Am Ende des Jupyter-Notebooks simulieren als Beispiel die Sonnenfinsternis 2026/27, als auch das n-Körperproblem der Planetenbewegungen unseres Sonnensystems.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from matplotlib import rcParams
import matplotlib.gridspec as gridspec

Wir definieren nun das System der gekoppelten Differenzialgleichung. Die Werte der Massenparameter $m_1$, $m_2$ und $m_3$ und die Gravitationskonstante $G$ legen wir noch nicht fest; diese werden als zusätzliche Argumente in der DGL-Funktion definiert.

def DGLsys(t,u_vec,G,m):
    # Entpacke u_vec in 9 Koordinaten und 9 Geschwindigkeiten
    r = u_vec[:9].reshape((3, 3))
    v = u_vec[9:].reshape((3, 3))
    # DGL-System erster Ordnung (18 Gleichungen)
    drdt = v
    dvdt = np.zeros_like(v)
    for i in range(3):
        for j in range(3):
            if i != j:
                dvdt[i] += G * m[j] * (r[j] - r[i]) / np.linalg.norm(r[j] - r[i])**3            
    return np.concatenate([drdt.flatten(), dvdt.flatten()])

Elliptische Lagrangesche homografischen Lösungen

Die elliptischen Lagrangeschen homografischen Lösungen im Dreikörperproblem wurden bereits im Jahre 1772 von Joseph-Louis Lagrange entdeckt. In diesen Lösungen bilden die drei Körper zu jedem Zeitpunkt die Eckpunkte eines gleichseitigen Dreiecks, wobei sich die Größe und Ausrichtung des Dreiecks im Laufe der Zeit ändern kann. Wir betrachten im Folgenden den Spezialfall gleicher Massen ($m_1 = m_2 = m_3 =1 $) und setzen die Gravitationskonstante auf Eins ($G=1$).

# Elliptische Lagrangesche homografischen Lösungen (Exzentrizität e)
set_G = 1
set_m = [1,1,1]

def u0_Lagrange(set_e = 0.5):
    r_1_0 = np.array([1/np.sqrt(3), 0, 0])
    r_2_0 = np.array([-1/(2*np.sqrt(3)), 1/2, 0])
    r_3_0 = np.array([-1/(2*np.sqrt(3)), -1/2, 0])
    v_1_0 = np.array([0, np.sqrt(1-set_e), 0])
    v_2_0 = np.array([-np.sqrt(3*(1-set_e))/2, -np.sqrt(1-set_e)/2, 0])
    v_3_0 = np.array([np.sqrt(3*(1-set_e))/2, -np.sqrt(1-set_e)/2, 0])

    u_0 = np.concatenate([r_1_0, r_2_0, r_3_0, v_1_0, v_2_0, v_3_0])
    return u_0

Mittels der oben definierten Funktion initialisieren wir den 18 Komponentigen Zustandsvektor $\vec{u}(0)$ mit ($u_1(0)=x_{1}(0)$, $u_2(0)=y_{1}(0)$, ..., $u_{18}(0)=v_{3z}(0)$) und betrachten speziell eine Anfangsbedingung mit einer Exzentrizität von $e=0.7$.

u_0 = u0_Lagrange(0.7)

Zum Lösen des gekoppelten Systems von Differenzialgleichungen benutzen wir das Python-Modul SciPy, welches eine breite Kollektion von mathematischen Algorithmen und Funktionen bereitstellt. Im Speziellen werden wir die Funktion 'solve_ivp(...)' verwenden, die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differenzialgleichungen bereitstellt.

Bei dem numerischen Lösen von Differenzialgleichungen in Python kann der Benutzer zusätzlich die relativen und absoluten Fehler-Toleranzen der berechneten numerischen Werte festlegen, wobei der relative Fehler mit der Zusatzoption 'rtol' und der absolute Fehler mit der Zusatzoption 'atol' kontrolliert wird (näheres siehe scipy.integrate.solve_ivp). Im Folgenden setzen wir die relativen und absoluten Fehler-Toleranzen (rtol und atol) auf $10^{-9}$.

Wir lösen das Anfangswertproblem im Bereich $t \in [0,4]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten (mesh points).

N = 10000
t_span = [0,4]
tval = np.linspace(0, t_span[1], N)
fehler = 1e-9

Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)

Wir veranschaulichen uns nun die numerischen Ergebnisse in unterschiedlichen Diagrammen.

rcParams.update({
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
})

plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel("x-Koordinate")
plt.plot(Loes.t,Loes.y[0],color="red", label=r"$\rm x_{1}(t)$")
plt.plot(Loes.t,Loes.y[3],color="blue", label=r"$\rm x_{2}(t)$")
plt.plot(Loes.t,Loes.y[6],color="green", label=r"$\rm x_{3}(t)$")
plt.legend(loc='upper right',fontsize=14);

Die obere Abbildung zeigt, wie sich die x-Koordinaten der drei Körper $x_{1}$, $x_{2}$ und $x_{3}$ im Laufe der Zeit verändern und die untere Abbildung zeigt das zeitliche Verhalten der y-Koordinaten. Man erkennt eine sehr gleichförmigen und periodischen Verlauf.

plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel("y-Koordinate")
plt.plot(Loes.t,Loes.y[1],color="red", label=r"$\rm y_{1}(t)$")
plt.plot(Loes.t,Loes.y[4],color="blue", label=r"$\rm y_{2}(t)$")
plt.plot(Loes.t,Loes.y[7],color="green", label=r"$\rm y_{3}(t)$")
plt.legend(loc='upper right',fontsize=14);

Die folgende Abbildung zeigt die Trajektorien der Bewegungen der drei Körper im Konfigurationsraum $(x,y)$. Die Punkte auf den Kurven zeigen die Anfangskonfigurationen der drei Körper. Anhand der grauen Verbindungslinien erkennt man, dass die Anfangskonfiguration die Form eines gleichseitigen Dreiecks hat.

plt.figure(figsize=(5, 5))
plt.xlabel("x-Koordinate")
plt.ylabel("y-Koordinate")
plt.plot(Loes.y[0],Loes.y[1],color="red", label=r"$\rm \vec{r}_1(t)$")
plt.plot(Loes.y[3],Loes.y[4],color="blue", label=r"$\rm \vec{r}_2(t)$")
plt.plot(Loes.y[6],Loes.y[7],color="green", label=r"$\rm \vec{r}_3(t)$")
plt.scatter(Loes.y[0][0],Loes.y[1][0],color="red")
plt.scatter(Loes.y[3][0],Loes.y[4][0],color="blue")
plt.scatter(Loes.y[6][0],Loes.y[7][0],color="green")
plt.plot([Loes.y[0][0],Loes.y[3][0]], [Loes.y[1][0],Loes.y[4][0]], color='grey', alpha=0.2)
plt.plot([Loes.y[0][0],Loes.y[6][0]], [Loes.y[1][0],Loes.y[7][0]], color='grey', alpha=0.2)
plt.plot([Loes.y[6][0],Loes.y[3][0]], [Loes.y[7][0],Loes.y[4][0]], color='grey', alpha=0.2)
plt.legend(loc='upper right',fontsize=14);
plt.axis('equal');

Man erkennt, dass sich die drei Körper ellipsenförmig um ihren gemeinsamen Schwerpunkt bewegen. Während der Bewegung bleibt die Gesamtenergie und der gesamte Drehimpuls des Systems erhalten. Die kinetische und potentielle Energie des Gesamtsystems berechnet sich wie folgt:

\begin{equation} E = E_{kin} + E_{pot} \quad , \quad E_{kin} = \frac{1}{2} \sum_{i=1}^3 m_i \vec{v}_i \cdot \vec{v}_i \quad , \quad E_{pot} = -\frac{G}{2} \sum_{i=1}^3 \sum_{j=1\\i \neq j}^3 \frac{m_i \, m_j}{\left| \, \vec{r}_j - \vec{r}_i \right|} \end{equation}

def E(sol, G, m):
    E_kin, E_pot = 0, 0
    for i in range(3):
        v = np.sqrt(sol.y[i*3+9]**2 + sol.y[i*3+10]**2 + sol.y[i*3+11]**2)
        E_kin += 1/2 * m[i] * v**2
        for j in range(3):
            if i != j:
                rji = np.sqrt((sol.y[j*3]-sol.y[i*3])**2 + (sol.y[j*3+1]-sol.y[i*3+1])**2 + (sol.y[j*3+2]-sol.y[i*3+2])**2)
                E_pot += -G/2 * m[i] * m[j] / rji
    return E_kin, E_pot

Der Drehimpuls des Gesamtsystems berechnet sich wie folgt

\begin{equation} \vec{L} = \sum_{i=1}^3 \vec{L}_i = \sum_{i=1}^3 m_i \vec{r}_i \times \vec{v}_i \end{equation}

def L_j(sol, m, j):
    Li = []
    L = 0
    for i in range(3):
        Li.append(m[i] * np.cross([sol.y[i*3][j],sol.y[i*3+1][j],sol.y[i*3+2][j]], [sol.y[i*3+9][j],sol.y[i*3+10][j],sol.y[i*3+11][j]]))
        L += Li[i]
    return Li, L
def L_g(sol, m):
    Lig = []
    Lg = []
    for j in range(len(sol.t)):
        Li, L = L_j(sol, m, j)
        Lg.append(L)
        Lig.append(Li)
    return np.array(Lig), np.array(Lg)

Die untere Funktion stellt eine optimierte Version der oberen Python-Funktion dar. Beide Funktionen berechnen jedoch das gleiche und geben

def L_g_short(sol, m):
    r = sol.y[:9].reshape(3, 3, -1).transpose(0, 2, 1)
    v = sol.y[9:18].reshape(3, 3, -1).transpose(0, 2, 1)
    m = np.array(m)
    Li = m[:, None, None] * np.cross(r, v)
    L_total = Li.sum(axis=0)
    return Li.transpose(1, 0, 2), L_total

Mittels der oben definierten Python-Funktionen berechnen wir für die gesamten Simulationswerte die Energiekomponenten $E$, $E_{kin}$ und $E_{pot}$, den gesamten Drehimpuls $\vec{L}$ und die Drehimpulse der einzelnen Körper $\vec{L}_1$, $\vec{L}_2$ und $\vec{L}_3$.

E_kin, E_pot = E(Loes, set_G, set_m)
Li, L = L_g_short(Loes, set_m)

Wir visualisieren nun die berechneten Werte.

plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel(r"Energie $E$")
plt.plot(Loes.t, E_kin + E_pot,c="black", label=r"Gesamtenergie $E$")
plt.plot(Loes.t, E_kin,c="brown", label=r"kinetische Energie $E_{kin}$")
plt.plot(Loes.t, E_pot,c="orange", label=r"potentielle Energie $E_{pot}$")
plt.legend(loc='upper right',fontsize=12);

In der oberen Abbildung erkennt man, dass sich die potentielle und kinetische Energie periodisch so verändern, dass die Gesamtenergie des Systems erhalten bleibt.

plt.figure(figsize=(11, 5))
plt.xlabel(r"t")
plt.ylabel(r"$L_x, L_y, L_z$")
plt.plot(Loes.t, L[:,0],c="brown", label=r"Gesamter Drehimpuls $L_x $")
plt.plot(Loes.t, L[:,1],c="orange", label=r"Gesamter Drehimpuls $L_y $")
plt.plot(Loes.t, L[:,2],c="black", label=r"Gesamter Drehimpuls $L_z$")
plt.plot(Loes.t, Li[:,0,2],c="red", label=r"Drehimpuls $L_{1z}$")
plt.plot(Loes.t, Li[:,1,2],c="blue", label=r"Drehimpuls $L_{2z}$")
plt.plot(Loes.t, Li[:,2,2],c="green", label=r"Drehimpuls $L_{3z}$")
plt.legend(loc='upper right',fontsize=10);

In der oberen Abbildung sind die Drehimpulskomponenten dargestellt. Der gesamte Drehimpuls besitzt nur eine nichtverschwindende Komponente in z-Richtung $\vec{L}=(\left(0, 0, L_z \right)$, mit $L_z \approx 0.94868$. Dieser Wert setzt sich aus gleichen Teilen aus den zeitlich konstanten Beiträgen der Drehimpulse der einzelnen Körper zusammen ($L_{iz} \approx 0.3162278$), wobei sich die jeweiligen Kurven in der oberen Abbildung überlagern.

Neben der Visualisierung der Bewegung der einzelnen Körper im Konfigurationsraum $(x,y)$ ist es oft sinnvoll die zeitliche Entwicklung auch als Phasenraumtrajektorien darzustellen. Für jeden der Körper lassen sich z.B. die kartesischen Trajektorien in den Phasenräumen $(v_x,x)$, $(v_y,y)$ und $(v_z,z)$ darstellen:

fig = plt.figure(figsize=(10, 5))
gs = gridspec.GridSpec(2, 3, width_ratios=[1, 1, 1], hspace=0.4, wspace=0.4)

ax1 = fig.add_subplot(gs[0, 0]) # Oben links
ax2 = fig.add_subplot(gs[0, 1]) # Oben mitte
ax3 = fig.add_subplot(gs[0, 2]) # Oben rechts
ax4 = fig.add_subplot(gs[1, 0]) # Unten links
ax5 = fig.add_subplot(gs[1, 1]) # Unten mitte
ax6 = fig.add_subplot(gs[1, 2]) # Unten rechts

ax1.set_xlabel(r"$r_{1x}$")
ax1.set_ylabel(r"$v_{1x}$")
ax1.plot(Loes.y[0],Loes.y[9],color="red")
ax1.scatter(Loes.y[0][0],Loes.y[9][0],color="red")

ax2.set_xlabel(r"$r_{2x}$")
ax2.set_ylabel(r"$v_{2x}$")
ax2.plot(Loes.y[3],Loes.y[12],color="blue")
ax2.scatter(Loes.y[3][0],Loes.y[12][0],color="blue")

ax3.set_xlabel(r"$r_{3x}$")
ax3.set_ylabel(r"$v_{3x}$")
ax3.plot(Loes.y[6],Loes.y[15],color="green")
ax3.scatter(Loes.y[6][0],Loes.y[15][0],color="green")


ax4.set_xlabel(r"$r_{1y}$")
ax4.set_ylabel(r"$v_{1y}$")
ax4.plot(Loes.y[1],Loes.y[10],color="red")
ax4.scatter(Loes.y[1][0],Loes.y[10][0],color="red")

ax5.set_xlabel(r"$r_{2y}$")
ax5.set_ylabel(r"$v_{2y}$")
ax5.plot(Loes.y[4],Loes.y[13],color="blue")
ax5.scatter(Loes.y[4][0],Loes.y[13][0],color="blue")

ax6.set_xlabel(r"$r_{3y}$")
ax6.set_ylabel(r"$v_{3y}$")
ax6.plot(Loes.y[7],Loes.y[16],color="green")
ax6.scatter(Loes.y[7][0],Loes.y[16][0],color="green");

Diese zeitliche Entwicklung der Körper im Phasenraum kann man auch in radialen Zentrumskoordinaten visualisieren. Dazu transformieren wir die kartesischen Trajektorien der Körper in Koordinaten, die die radiale Entfernung der Körper vom gemeinsamen Zentrum beschreiben ($r_1, r_2$ und $r_3$). Zusätzlich berechnen wir auch die radialen Geschwindigkeiten der Körper relativ zum gemeinsamen Zentrum ($v_1, v_2$ und $v_3$).

Der Schwerpunkt $\vec{r}_S=\left( x_S, y_S, z_S \right)$ und die Geschwindigkeit des Schwerpunkts $\vec{v}_S=\left( v_{xS}, v_{yS}, v_{zS} \right)$ berechnen sich wie folgt:

\begin{equation} \vec{r}_S = \frac{\sum_{i=1}^3 m_i \vec{r}_i}{\sum_{i=1}^3 m_i} \quad \hbox{und} \quad \vec{v}_S = \frac{\sum_{i=1}^3 m_i \vec{v}_i}{\sum_{i=1}^3 m_i} \end{equation}

Der Abstand jedes Körpers zum gemeinsamen Schwerpunkt ($R_i$) ist die euklidische Norm der relativen Positionsvektoren:

\begin{equation} R_i = \sqrt{ \left( r_{ix} - x_S \right)^2 + \left( r_{iy} - y_S \right)^2 + \left( r_{iz} - z_S \right)^2} \end{equation}

Die radiale Geschwindigkeit $V_i$ des Körpers $i$ ist die Projektion des relativen Geschwindigkeitsvektors auf den Position-Einheitsvektor:

\begin{equation} V_i = \left( \vec{v}_{i} - \vec{v}_S \right) \cdot \frac{\left( \vec{r}_{i} - \vec{r}_S \right)}{r_i} \end{equation}

Die zwei folgenden Python-Funktionen berechnen die gleichen oben definierten Größen, wobei die zweite Version mittels der Einsteinschen Summenkonvention das Skalarprodukt der Projektion berechnet (elegantere Variante).

# Koordinaten-Transformation: Schwerpunkt, r, v_r
def Transf_S(Loes, m):
    # Schwerpunkt
    x_S = (m[0]*Loes.y[0] + m[1]*Loes.y[3] + m[2]*Loes.y[6]) / sum(m)
    y_S = (m[0]*Loes.y[1] + m[1]*Loes.y[4] + m[2]*Loes.y[7]) / sum(m)
    z_S = (m[0]*Loes.y[2] + m[1]*Loes.y[5] + m[2]*Loes.y[8]) / sum(m)
    # Geschwindigkeit Schwerpunkt
    v_x_S = (m[0]*Loes.y[9] + m[1]*Loes.y[12] + m[2]*Loes.y[15]) / sum(m)
    v_y_S = (m[0]*Loes.y[10] + m[1]*Loes.y[13] + m[2]*Loes.y[16]) / sum(m)
    v_z_S = (m[0]*Loes.y[11] + m[1]*Loes.y[14] + m[2]*Loes.y[17]) / sum(m)
    # Relativer Radius
    r_1 = np.sqrt((Loes.y[0]-x_S)**2 + (Loes.y[1]-y_S)**2 + (Loes.y[2]-z_S)**2)
    r_2 = np.sqrt((Loes.y[3]-x_S)**2 + (Loes.y[4]-y_S)**2 + (Loes.y[5]-z_S)**2)
    r_3 = np.sqrt((Loes.y[6]-x_S)**2 + (Loes.y[7]-y_S)**2 + (Loes.y[8]-z_S)**2)
    # Radiale Geschwindigkeit 
    vr_1 = ((Loes.y[0]-x_S)*(Loes.y[9]-v_x_S) + (Loes.y[1]-y_S)*(Loes.y[10]-v_y_S) + (Loes.y[2]-z_S)*(Loes.y[11]-v_z_S)) / r_1
    vr_2 = ((Loes.y[3]-x_S)*(Loes.y[12]-v_x_S) + (Loes.y[4]-y_S)*(Loes.y[13]-v_y_S) + (Loes.y[5]-z_S)*(Loes.y[14]-v_z_S)) / r_2
    vr_3 = ((Loes.y[6]-x_S)*(Loes.y[15]-v_x_S) + (Loes.y[7]-y_S)*(Loes.y[16]-v_y_S) + (Loes.y[8]-z_S)*(Loes.y[17]-v_z_S)) / r_3
    return [x_S,y_S,z_S],[r_1,r_2,r_3],[vr_1,vr_2,vr_3]

# Koordinaten-Transformation: Schwerpunkt, r, v_r (elegantere Variante)
def Transf_S(Loes, m):
    # Daten restrukturieren: (3 Körper, 3 Dimensionen)
    # y[0:9] sind Positionen, y[9:18] sind Geschwindigkeiten
    pos = Loes.y[0:9].reshape(3, 3, -1)
    vel = Loes.y[9:18].reshape(3, 3, -1)
    # Massen für Multiplikation (3 Körper, 1 Achse, 1 Zeit)
    m_arr = np.array(m).reshape(3, 1, 1)
    
    # Schwerpunkt (Zentrum der Masse) berechnen
    pos_S = np.sum(m_arr * pos, axis=0) / np.sum(m)
    vel_S = np.sum(m_arr * vel, axis=0) / np.sum(m) # Geschwindigkeit des Schwerpunkts
    
    # Relative Positionen und Geschwindigkeiten zum Schwerpunkt
    rel_pos = pos - pos_S
    rel_vel = vel - vel_S 
    
    # Radien Positionsvektoren der drei Körper
    r = np.linalg.norm(rel_pos, axis=1)
    
    # Radiale Geschwindigkeit (Skalarprodukt / r)
    # np.einsum: "Einsteinsche Summenkonvention" für das Skalarprodukt
    # v_r mit np.einsum (über Achse 1 summieren, Achse 0 und 2 behalten)
    # 'ijk,ijk->ik' bedeutet: (Körper i, Koordinate j, Zeit k)
    vr = np.einsum('ijk,ijk->ik', rel_pos, rel_vel) / (r + 1e-15)
    
    return pos_S, r, vr

Mittels einer der oben definierten Python-Funktionen berechnen wir für die gesamten Simulationswerte die Schwerpunktskoordinaten $\vec{r}_S=\left( x_S, y_S, z_S \right)$, die relativen Abstände $R_i$ der Körper zum gemeinsamen Schwerpunkt und die radialen Geschwindigkeiten $V_i$ der Körper.

r_S, R, Vr = Transf_S(Loes,set_m)

Wir Visualisieren nun die Bewegung der Körper und stellen die Trajektorien auch im Phasenraum der radialen Zentrumskoordinaten dar.

from matplotlib.animation import FuncAnimation
from IPython.display import HTML

rcParams.update({
    'axes.titlesize' : 10,
    'axes.labelsize' : 10,  
    'xtick.labelsize' : 10 ,
    'ytick.labelsize' : 10 ,
    'legend.fontsize' : 10
})

Die folgende Python-Funktion stellt die Bewegung der Körper in einer Animation dar.

def create_animation(Loes, R, Vr, r_S, N, frames_N=100):
    # Berechnung der zu Animation verwendeten Indices
    indices = [int(f * (N - 1) / (frames_N - 1)) for f in range(frames_N)]
    
    fig = plt.figure(figsize=(11, 5))
    gs = gridspec.GridSpec(2, 3, width_ratios=[2, 1, 1], hspace=0.3, wspace=0.4)
    
    ax1 = fig.add_subplot(gs[:, 0]) # Links (erstreckt sich über beide Reihen)
    ax2 = fig.add_subplot(gs[0, 1]) # Oben mitte
    ax3 = fig.add_subplot(gs[0, 2]) # Oben rechts
    ax4 = fig.add_subplot(gs[1, 1]) # Unten mitte
    ax5 = fig.add_subplot(gs[1, 2]) # Unten rechts
    
    # Statische Plots
    ax1.axis('equal')
    ax1.set_xlabel("x-Koordinate")
    ax1.set_ylabel("y-Koordinate")
    ax1.plot(Loes.y[0], Loes.y[1], color="red", label=r"$\rm \vec{r}_1(t)$")
    ax1.plot(Loes.y[3], Loes.y[4], color="blue", label=r"$\rm \vec{r}_2(t)$")
    ax1.plot(Loes.y[6], Loes.y[7], color="green", label=r"$\rm \vec{r}_3(t)$")
    ax1.legend(loc='upper right')
    ax1.axis('equal')
    
    ax2.set_xlabel(r"Radius $r_1$")
    ax2.set_ylabel(r"radiale Geschwind. $v_{r_1}$")
    ax2.plot(R[0], Vr[0], color="red")
    
    ax3.set_xlabel(r"Radius $r_2$")
    ax3.set_ylabel(r"radiale Geschwind. $v_{r_2}$")
    ax3.plot(R[1], Vr[1], color="blue")
    
    ax4.set_xlabel(r"Radius $r_3$")
    ax4.set_ylabel(r"radiale Geschwind. $v_{r_3}$")
    ax4.plot(R[2], Vr[2], color="green")
    
    ax5.set_xlabel("x-Koordinate")
    ax5.set_ylabel("y-Koordinate")
    ax5.plot(r_S[0], r_S[1], color="black", label=r"$\rm Schwerpunkt$")
    ax5.legend(loc='upper left')
    
    # Animationsobjekte
    dot_s = 30
    dot_r1 = ax1.scatter([], [], color='red', s=dot_s)
    dot_r2 = ax1.scatter([], [], color='blue', s=dot_s)
    dot_r3 = ax1.scatter([], [], color='green', s=dot_s)
    dot_r1a = ax2.scatter([], [], color='red', s=dot_s)
    dot_r2a = ax3.scatter([], [], color='blue', s=dot_s)
    dot_r3a = ax4.scatter([], [], color='green', s=dot_s)
    
    line_S1, = ax1.plot([], [], color='grey', alpha=0.2)
    line_S2, = ax1.plot([], [], color='grey', alpha=0.2)
    line_S3, = ax1.plot([], [], color='grey', alpha=0.2)
    
    def update(frame):
        i = indices[frame]
        dot_r1.set_offsets([[Loes.y[0][i], Loes.y[1][i]]])
        dot_r2.set_offsets([[Loes.y[3][i], Loes.y[4][i]]])
        dot_r3.set_offsets([[Loes.y[6][i], Loes.y[7][i]]])
        dot_r1a.set_offsets([[R[0][i], Vr[0][i]]])
        dot_r2a.set_offsets([[R[1][i], Vr[1][i]]])
        dot_r3a.set_offsets([[R[2][i], Vr[2][i]]])
        
        line_S1.set_data([r_S[0][i], Loes.y[0][i]], [r_S[1][i], Loes.y[1][i]])
        line_S2.set_data([r_S[0][i], Loes.y[3][i]], [r_S[1][i], Loes.y[4][i]])
        line_S3.set_data([r_S[0][i], Loes.y[6][i]], [r_S[1][i], Loes.y[7][i]])
        
        return dot_r1, dot_r2, dot_r3, dot_r1a, dot_r2a, dot_r3a, line_S1, line_S2, line_S3

    ani = FuncAnimation(fig, update, frames=frames_N, interval=70, blit=True)
    plt.close(fig) 
    return ani

Die untere Version stellt die gleiche Animation dar, ist jedoch kompakter geschrieben.

def create_animation_short(Loes, R, Vr, r_S, N, frames_N=100):
    # Berechnung der zu Animation verwendeten Indices
    indices = [int(f * (N - 1) / (frames_N - 1)) for f in range(frames_N)]
    
    fig = plt.figure(figsize=(11, 5))
    gs = gridspec.GridSpec(2, 3, width_ratios=[2, 1, 1], hspace=0.3, wspace=0.4)
    
    # Achsen-Setup
    ax_main = fig.add_subplot(gs[:, 0])
    axes_sub = [fig.add_subplot(gs[0, 1]), fig.add_subplot(gs[0, 2]), fig.add_subplot(gs[1, 1])]
    ax_cm = fig.add_subplot(gs[1, 2])

    colors = ['red', 'blue', 'green']
    labels = [r"$\rm \vec{r}_1(t)$", r"$\rm \vec{r}_2(t)$", r"$\rm \vec{r}_3(t)$"]    
    dots_main, dots_sub, lines_cm = [], [], []

    # Hauptplot (Trajektorien im im Konfigurationsraum (x,y))
    ax_main.axis('equal')
    for j in range(3):
        idx_x, idx_y = j * 3, j * 3 + 1
        ax_main.plot(Loes.y[idx_x], Loes.y[idx_y], color=colors[j], label=labels[j], alpha=0.3)
        dots_main.append(ax_main.scatter([], [], color=colors[j], s=30))
        lines_cm.append(ax_main.plot([], [], color='grey', alpha=0.2)[0])
    ax_main.legend(loc='upper right')

    # Phasenraum-Plots (R vs Vr)
    for j, ax in enumerate(axes_sub):
        ax.plot(R[j], Vr[j], color=colors[j], alpha=0.3)
        dots_sub.append(ax.scatter([], [], color=colors[j], s=30))
        ax.set_xlabel(rf"$r_{j+1}$")
        ax.set_ylabel(rf"$v_{{r_{j+1}}}$")

    # Schwerpunktplot
    ax_cm.plot(r_S[0], r_S[1], color="black", label="Schwerpunkt")
    ax_cm.legend(loc='upper left')

    def update(frame):
        i = indices[frame]
        changed_objects = []
        
        for j in range(3):
            # Update Punkte der Körper im Hauptplot
            x, y = Loes.y[j*3][i], Loes.y[j*3+1][i]
            dots_main[j].set_offsets([[x, y]])
            
            # Update Verbindungslinien zum Schwerpunkt
            lines_cm[j].set_data([r_S[0][i], x], [r_S[1][i], y])
            
            # Update Punkte in den Phasenraum-Plots
            dots_sub[j].set_offsets([[R[j][i], Vr[j][i]]])
            
            changed_objects.extend([dots_main[j], lines_cm[j], dots_sub[j]])
            
        return changed_objects

    ani = FuncAnimation(fig, update, frames=frames_N, interval=80, blit=True)
    plt.close(fig)
    return ani

ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Der Schwerpunkt der Körper ist im Ursprung und verändert sich mit der Zeit nicht (siehe Diagramm rechts unten, beachte x-y-Skala).

Wir möchten uns nun die Bewegung einer weiteren elliptischen Lagrangeschen Lösung des Dreikörperproblems visualisieren und betrachten eine Anfangsbedingung mit einer größeren Exzentrizität ($e=0.9$).

u_0 = u0_Lagrange(0.9)

t_span = [0,4]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)

r_S, R, Vr = Transf_S(Loes,set_m)

ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Im Folgenden betrachten wir noch die Spezialfälle $e=0$ und $e=1$.

u_0 = u0_Lagrange(0)
t_span = [0,4]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

u_0 = u0_Lagrange(1)
t_span = [0,0.641]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Fehlende Stabilität und chaotisches Verhalten

Die dargestellten elliptischen und kreisförmigen Lagrange-Lösung sind jedoch instabil gegen kleine Änderungen. Die numerischen Lösungen zeigen somit alle, ab einem gewissen Zeitpunkt chaotisches Verhalten. Wir möchten dies exemplarisch an der Lösung mit $e=0.7$ zeigen.

u_0 = u0_Lagrange(0.7)

Wir lösen das Anfangswertproblem nun im Bereich $t \in [0,12.5]$ unter Verwendung von den Fehler-Toleranzen rtol und atol gleich $10^{-13}$.

N = 10000
t_span = [0,12.5]
tval = np.linspace(0, t_span[1], N)
fehler = 1e-13
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)
r_S, R, Vr = Transf_S(Loes,set_m)

ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Man erkennt, dass die periodische reguläre Bewegung am Ende der Simulation instabil wird und sich der blaue Körper zum Zeitpunkt $t \approx 12$ von den beiden anderen separiert. Die weitere Entwicklung folgt dann einem chaotischen Verhalten. Ob reguläre Lösungen stabil gegen Störungen sind, hängt unter anderem vom Massenverhältnis der Körper und von der Exzentrizität der Ellipse ab. Bei den Lagrange-Lösungen ist Stabilität nur möglich, wenn eine Masse deutlich größer als die anderen ist. Wir werden chaotisches Verhalten später näher betrachten.

Die Achterschleife (Figure-8)

Die zuvor dargestellten elliptischen Lagrange-Lösungen von drei Körpern gleicher Masse sind reguläre Sonderfälle, die jedoch relativ schnell in numerischen Simulationen in chaotisches Verhalten übergehen, da ihre Bewegung nur bei kleinsten Abweichungen von ihrer Idealbahn instabil wird. Die im Folgenden dargestellte Lösung besitzt hingegen numerische Stabilität und stellt eine wunderschöne, hochsymmetrische und gleichzeitig robuste Inseln der Ordnung inmitten der sonst oft chaotischen Bewegungen des Dreikörperproblems dar.

Wir betrachten nun eine weitere reguläre planare Bewegung, die in der Literatur als die "Figur Acht" bekannt ist, da alle drei Körper exakt der gleichen (Figur-8)-Flugbahn folgen, jedoch zeitlich versetzt (eine sogenannte choreografische Bewegung). Diese stabile Dreikörper-Bahn wurde 1993 von dem Physiker Cris Moore durch numerische Simulationen gefunden und im Jahre 2000 bewiesen die Mathematiker Alain Chenciner und Richard Montgomery die formale Existenz dieser Lösung. Die Figur-8-Lösung war eine mathematische Sensation, da sie einerseits die erste echte 3-Körper-Choreografie war und andererseits bei gleichen Massenwerten Stabilität zeigte. Zusätzlich war schwingen die drei Körper so umeinander, dass das System als Ganzes nicht rotiert (Gesamter Drehimpuls ist Null).

Wir setzen die Gravitationskonstante auf Eins ($G=1$) und die Massenwerte auf $m_1 = m_2 = m_3 =1 $.

# Die reguläre Bewegung der "Chenciner-Montgomery" Acht ("Figur Acht")
set_G = 1
set_m = [1,1,1]

def u0_Figur8():
    r_1_0 = np.array([ 0.97000436, -0.24308753, 0])
    r_2_0 = np.array([-r_1_0[0], -r_1_0[1], 0])
    r_3_0 = np.array([0, 0, 0])
    v_1_0 = np.array([0.466203685, 0.43236573, 0])
    v_2_0 = np.array([v_1_0[0], v_1_0[1], 0])
    v_3_0 = np.array([-2*v_1_0[0], -2*v_1_0[1], 0])
    u_0 = np.concatenate([r_1_0, r_2_0, r_3_0, v_1_0, v_2_0, v_3_0])
    return u_0

Die Bewegung der Körper stellen wir zunächst in einer Animation dar.

u_0 = u0_Figur8()
t_span = [0,6.32591398]
tval = np.linspace(0, t_span[1], N)
fehler = 1e-9
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Mittels der weiter oben definierten Python-Funktionen berechnen wir wider für die gesamten Simulationswerte die Energiekomponenten $E$, $E_{kin}$ und $E_{pot}$, den gesamten Drehimpuls $\vec{L}$ und die Drehimpulse der einzelnen Körper $\vec{L}_1$, $\vec{L}_2$ und $\vec{L}_3$.

E_kin, E_pot = E(Loes, set_G, set_m)
Lj, L = L_g_short(Loes, set_m)

Wir stellen diese Größen in den beiden folgenden Abbildungen dar.

plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel(r"Energie $E$")
plt.plot(Loes.t, E_kin + E_pot,c="black", label=r"Gesamtenergie $E$")
plt.plot(Loes.t, E_kin,c="brown", label=r"kinetische Energie $E_{kin}$")
plt.plot(Loes.t, E_pot,c="orange", label=r"potentielle Energie $E_{pot}$")
plt.legend(loc='upper right',fontsize=12);

plt.figure(figsize=(11, 5))
plt.xlabel(r"t")
plt.ylabel(r"$L_x, L_y, L_z$")
plt.plot(Loes.t, L[:,0],c="brown", label=r"Gesamter Drehimpuls $L_x $")
plt.plot(Loes.t, L[:,1],c="orange", label=r"Gesamter Drehimpuls $L_y $")
plt.plot(Loes.t, L[:,2],c="black", label=r"Gesamter Drehimpuls $L_z$")
plt.plot(Loes.t, Lj[:,0,2],c="red", label=r"Drehimpuls $L_{1z}$")
plt.plot(Loes.t, Lj[:,1,2],c="blue", label=r"Drehimpuls $L_{2z}$")
plt.plot(Loes.t, Lj[:,2,2],c="green", label=r"Drehimpuls $L_{3z}$")
plt.legend(loc='upper right',fontsize=10);

Wie schon oben erwähnt wurde, summieren sich die Drehimpuls der einzelnen Körper zu Null.

Weitere reguläre planare Bewegung

Gerade in den letzten Jahren wurden viele neue reguläre Bewegungen im Dreikörperproblem gefunden. So konnten z.B. die Wissenschaftler Xiaoming Li und Shijun Liao von der Shanghai Jiao Tong Universität in China mittels Methoden der künstlichen Intelligenz (maschinellem Lernen) viele neue periodische Umlaufbahnen im bekannten Dreikörperproblem identifizieren (siehe Literatur weiter unten). Im Gegensatz zu den bisher dargestellten Lösungen betrachten wir nun Systeme mit ungleichen Massenwerten, bei denen Stabilität im Allgemeinen einfacher zu finden ist. Die Anfangsbedingungen der nun dargestellten regulären planaren Bewegungen und noch weitere sind von der folgenden Internetseite erreichbar: Periodic orbits for Newtonian planar three-body problem with unequal mass

Literatur:

Liao, Shijun, Xiaoming Li, and Yu Yang. "Three-body problem - From Newton to supercomputer plus machine learning." New Astronomy 96 (2022): 101850.

Liao, Shijun, and Xiaoming Li. "On the periodic solutions of the three-body problem." National Science Review 6.6 (2019): 1070-1071.

Mittels der folgenden Python-Funktion kann der Benutzer eine von sieben Anfangsbedingungen des oben definierten Artikels auswählen. Die Funktion gibt die Liste der Massen der drei Körper und die Anfangsbedingungen zurück.

# Anfangsbedingungen der planaren Bewegungen von
# https://numericaltank.sjtu.edu.cn/three-body/three-body-unequal-mass.htm
set_G = 1

def u0_Article_1(num):
    match num:
        case 1:
            m1, m2, m3, x1, v1, v2 = 0.1000, 0.2300, 1.0000, -1.371949198659973e+00, -8.862881132488540e-01, -8.279658998898554e-01
        case 2:
            m1, m2, m3, x1, v1, v2 = 0.3300, 3.1500, 1.0000, -8.525335038973433e-01, -1.633511563063199e+00, -3.704384551900162e-01
        case 3:
            m1, m2, m3, x1, v1, v2 = 0.6700, 3.5900, 1.0000, -8.797938089126017e-01, -1.630970837238929e+00, -2.441737423450410e-01
        case 4:
            m1, m2, m3, x1, v1, v2 = 0.1000, 0.1000, 1.0000, -1.349417050308199e+00, -7.762940193673036e-01, -4.028136981897485e-01
        case 5:
            m1, m2, m3, x1, v1, v2 = 2.0000, 0.9000, 1.0000, -1.882679774436631e+00, -5.399521505705154e-01, 1.982242573538665e-02
        case 6:
            m1, m2, m3, x1, v1, v2 = 6.5000, 23.2000, 1.0000, -1.382950469514427e+00, -2.054567296895344e+00, 3.972991955422785e-01
        case 7:
            m1, m2, m3, x1, v1, v2 = 27.4000, 3.9000, 1.0000, -2.941703082411299e+00, -4.188649992386887e-01, 1.908291421590212e+00
        case _:
            m1, m2, m3, x1, v1, v2 = 1.1900, 1.6400, 1.0000, -1.202078235710311e+00, -1.045592652499898e+00, -1.744038482791141e-01 
            
    r_1_0 = np.array([x1, 0, 0])
    r_2_0 = np.array([1, 0, 0])
    r_3_0 = np.array([0, 0, 0])
    v_1_0 = np.array([0, v1, 0])
    v_2_0 = np.array([0, v2, 0])
    v_3_0 = np.array([0, -(m1*v1 + m2*v2)/m3, 0])
    
    set_m = [m1, m2, m3]
    u_0 = np.concatenate([r_1_0, r_2_0, r_3_0, v_1_0, v_2_0, v_3_0])
    return set_m, u_0

In der folgenden Abbildung sind die Trajektorien von sechs unterschiedlichen reguläre planare Bewegung dargestellt.

sims = [1, 3, 4, 5, 6, 7] 
t_span = [0, 12]
tval = np.linspace(0, t_span[1], N)

Loes_i = {}
for i in sims:
    set_m, u_0 = u0_Article_1(i)
    Loes_i[i] = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G, set_m), rtol=fehler, atol=fehler)

fig, axes = plt.subplots(2, 3, figsize=(14, 8), constrained_layout=True)
axes = axes.flatten()
colors = ["red", "blue", "green"]
labels = [r"$\rm \vec{r}_1$", r"$\rm \vec{r}_2$", r"$\rm \vec{r}_3$"]

for k, i in enumerate(sims):
    ax = axes[k]
    Loes = Loes_i[i]
    
    # Die 3 Körper plotten
    for j in range(3):
        ax.plot(Loes.y[j*3], Loes.y[j*3+1], color=colors[j], alpha=0.8)
    
    ax.set_aspect('equal')
    ax.set_title(f"Konfiguration {i}", fontsize=14)
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.grid(True, linestyle=':', alpha=0.5)

fig.legend(labels, loc='upper center', ncol=3, bbox_to_anchor=(0.5, 1.05), frameon=False, fontsize=14)

plt.show()

Wir betrachten uns die Simulation der oberen, mittleren Abbildung genauer und stellen diese in einer Animation dar.

idx = 3
Loes = Loes_i[idx]
set_m, u_0 = u0_Article_1(idx)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Wir betrachten uns nun die Simulation der unteren, linken Abbildung genauer (siehe voriges Übersichtsdiagramm) und stellen diese ebenfalls in einer Animation dar.

idx = 5
Loes = Loes_i[idx]
set_m, u_0 = u0_Article_1(idx)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Der Verlauf der Trajektorien, sowohl im Konfigurationsraum, als auch in den radialen Phasenräumen der Körper, sieht bei dieser Bewegung komplizierter aus. Wir betrachten zusätzlich noch den Verlauf der Energie- und Drehimpulskomponenten.

E_kin, E_pot = E(Loes, set_G, set_m)
Lj, L = L_g_short(Loes, set_m)

plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel(r"Energie $E$")
plt.plot(Loes.t, E_kin + E_pot,c="black", label=r"Gesamtenergie $E$")
plt.plot(Loes.t, E_kin,c="brown", label=r"kinetische Energie $E_{kin}$")
plt.plot(Loes.t, E_pot,c="orange", label=r"potentielle Energie $E_{pot}$")
plt.legend(loc='upper right',fontsize=12);

plt.figure(figsize=(11, 5))
plt.xlabel(r"t")
plt.ylabel(r"$L_x, L_y, L_z$")
plt.plot(Loes.t, L[:,0],c="brown", label=r"Gesamter Drehimpuls $L_x $")
plt.plot(Loes.t, L[:,1],c="orange", label=r"Gesamter Drehimpuls $L_y $")
plt.plot(Loes.t, L[:,2],c="black", label=r"Gesamter Drehimpuls $L_z$")
plt.plot(Loes.t, Lj[:,0,2],c="red", label=r"Drehimpuls $L_{1z}$")
plt.plot(Loes.t, Lj[:,1,2],c="blue", label=r"Drehimpuls $L_{2z}$")
plt.plot(Loes.t, Lj[:,2,2],c="green", label=r"Drehimpuls $L_{3z}$")
plt.legend(loc='upper right',fontsize=10);

Wir betrachten nun abschließend die Simulation der unteren, rechten Abbildung genauer (siehe voriges Übersichtsdiagramm) und stellen diese ebenfalls in einer Animation dar.

idx = 7
Loes = Loes_i[idx]
set_m, u_0 = u0_Article_1(idx)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Das Dreikörperproblem und chaotische Bewegungen

Alle bisher dargestellten Bewegungen waren aufgrund der speziellen Wahl der Anfangswerte von periodischer, regulärer Struktur, wobei einige (z.B. die dargestellten elliptischen Lagrange-Lösungen) nach einiger Zeit chaotisches Verhalten zeigen. Bewegungen von Dreikörper-Systemen mit gleichen Massenwerten sind in der Regel chaotisch und wir möchten nun ein weiteres Beispiel einer chaotischen Bewegung im Detail vorstellen und verwenden dazu die folgenden Anfangswerte:

# Anfangswerte einer chaotischen Bewegung
set_m = [1,1,1]
set_G = 1
r_1_0 = np.array([1, 0, 0])
r_2_0 = np.array([-0.5, 0.8, 0])
r_3_0 = np.array([-0.5, -0.8, 0])
v_1_0 = np.array([0, 0.5, 0])
v_2_0 = np.array([0.2, -0.1, 0])
v_3_0 = np.array([-0.2, -0.4, 0])
u_0 = np.concatenate([r_1_0, r_2_0, r_3_0, v_1_0, v_2_0, v_3_0])

Wir lösen das Anfangswertproblem im Bereich $t \in [0,25]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten, setzen den relativen und absoluten Fehler (rtol und atol) auf $10^{-9}$ und stellen uns die Bewegung wieder in einer Animation dar.

fehler = 1e-9
t_span = [0,25]
N = 10000
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Im Gegensatz zu den bisher betrachteten Bewegungen, sehen nun die Trajektorien der Körper schon am Anfang chaotisch aus und es ist keine einfache Periodizität erkennbar. Zusätzlich separiert sich zum Zeitpunkt $t \approx 20$ der blaue Körper von den beiden anderen. Um die weitere zeitliche Entwicklung des Systems zu untersuchen, stellen wir Im Folgenden die Trajektorien der Bewegung für einen längeren Zeitraum ($t \in [0,45]$), unter Verwendung von $N=100000$ zeitlichen Gitterpunkten dar.

fehler = 1e-9
t_span = [0,45]
N = 100000
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_short(Loes, R, Vr, r_S, N, 150)
HTML(ani.to_html5_video())

Wie erwartet kehrt der blaue Körper zu einem späteren Zeitpunkt zu den anderen beiden zurück, separiert sich dann jedoch erneut. Wir möchten im Folgenden den chaotischen Charakter der Bewegung genauer studieren und fertigen dazu eine weitere numerische Lösung mit den gleichen Anfangsbedingungen an, jedoch setzen wir nun den relativen und absoluten Fehler (rtol und atol) auf $f=10^{-13}$.

fehler = 1e-13
t_span = [0,45]
N = 100000
tval = np.linspace(0, t_span[1], N)
Loes1 = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)

N1 = 0
plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel("x-Koordinate")
plt.plot(Loes.t[N1:N],Loes.y[1][N1:N],color="red", alpha=0.2, label=r"$r_{1x}(t), f=10^{-9}$")
plt.plot(Loes1.t[N1:N],Loes1.y[1][N1:N],color="red", linestyle=":", label=r"$r_{1x}(t), f=10^{-13}$")
plt.plot(Loes.t[N1:N],Loes.y[3][N1:N],color="blue", alpha=0.2, label=r"$\rm r_{2x}(t), f=10^{-9}$")
plt.plot(Loes1.t[N1:N],Loes1.y[3][N1:N],color="blue", linestyle=":", label=r"$\rm r_{2x}(t), f=10^{-13}$")
plt.plot(Loes.t[N1:N],Loes.y[6][N1:N],color="green", alpha=0.2, label=r"$\rm r_{3x}(t), f=10^{-9}$")
plt.plot(Loes1.t[N1:N],Loes1.y[6][N1:N],color="green", linestyle=":", label=r"$\rm r_{3x}(t), f=10^{-13}$")
plt.legend(loc='upper left',fontsize=12);

In der oberen Abbildung ist die Bewegung der drei Körper in x-Richtung ($r_{1x}, r_{2x}, r_{3x},$) für die beiden unterschiedlich genauen Simulationen dargestellt, wobei die gepunkteten Kurven die Ergebnisse der Simulation mit rtol=atol=$f=0^{-13}$ und die durchgezogenen, transparenten Kurven die Ergebnisse der Simulation mit rtol=atol=$f=10^{-9}$ darstellen. Man erkennt, dass sich die beiden Simulationen nach dem zweiten Entfernen des blauen Körpers unterscheiden. Die folgende Abbildung verdeutlicht dies nochmals, indem sie diesen Zeitbereich vergrößert.

N1 = int(0.88*N)
plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel("x-Koordinate")
plt.plot(Loes.t[N1:N],Loes.y[1][N1:N],color="red", alpha=0.2, label=r"$r_{1x}(t), f=10^{-9}$")
plt.plot(Loes1.t[N1:N],Loes1.y[1][N1:N],color="red", linestyle=":", label=r"$r_{1x}(t), f=10^{-13}$")
plt.plot(Loes.t[N1:N],Loes.y[3][N1:N],color="blue", alpha=0.2, label=r"$\rm r_{2x}(t), f=10^{-9}$")
plt.plot(Loes1.t[N1:N],Loes1.y[3][N1:N],color="blue", linestyle=":", label=r"$\rm r_{2x}(t), f=10^{-13}$")
plt.plot(Loes.t[N1:N],Loes.y[6][N1:N],color="green", alpha=0.2, label=r"$\rm r_{3x}(t), f=10^{-9}$")
plt.plot(Loes1.t[N1:N],Loes1.y[6][N1:N],color="green", linestyle=":", label=r"$\rm r_{3x}(t), f=10^{-13}$")
plt.legend(loc='upper left',fontsize=12);

Das chaotische Verhalten der Trajektorien bedeutet auch, dass sich kleine, $\epsilon$-ähnliche Störungen nach einiger Zeit so aufschaukeln, dass die Bewegung nicht mehr adäquat vorhergesagt werden kann. Dies wollen wir im Folgenden numerisch betrachten, indem wir eine Vergleichrechnung mit rtol=atol=$f=0^{-13}$ machen, bei der der $r_{1x}(t=0)$-Wert der Anfangsbedingungen um ein $\epsilon$ (hier $\epsilon=10^{-8}$) verschoben ist. Die untere Abbildung zeigt den Unterschied der simulierten Bewegungen.

fehler = 1e-13
t_span = [0,45]
N = 100000
tval = np.linspace(0, t_span[1], N)

u_0_eps = list(u_0)
u_0_eps[0] = u_0[0] + 1e-8

Loes2 = solve_ivp(DGLsys, t_span, u_0_eps, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)

N1 = 0
plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel("x-Koordinate")
plt.plot(Loes2.t[N1:N],Loes2.y[1][N1:N],color="red", alpha=0.2, label=r"$r_{1x}(t), \epsilon=10^{-8}$")
plt.plot(Loes1.t[N1:N],Loes1.y[1][N1:N],color="red", linestyle=":", label=r"$r_{1x}(t)$")
plt.plot(Loes2.t[N1:N],Loes2.y[3][N1:N],color="blue", alpha=0.2, label=r"$\rm r_{2x}(t), \epsilon=10^{-8}$")
plt.plot(Loes1.t[N1:N],Loes1.y[3][N1:N],color="blue", linestyle=":", label=r"$\rm r_{2x}(t)$")
plt.plot(Loes2.t[N1:N],Loes2.y[6][N1:N],color="green", alpha=0.2, label=r"$\rm r_{3x}(t), \epsilon=10^{-8}$")
plt.plot(Loes1.t[N1:N],Loes1.y[6][N1:N],color="green", linestyle=":", label=r"$\rm r_{3x}(t)$")
plt.legend(loc='upper left',fontsize=12);

N1 = int(0.85*N)
plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel("x-Koordinate")
plt.plot(Loes2.t[N1:N],Loes2.y[1][N1:N],color="red", alpha=0.2, label=r"$r_{1x}(t), \epsilon=10^{-8}$")
plt.plot(Loes1.t[N1:N],Loes1.y[1][N1:N],color="red", linestyle=":", label=r"$r_{1x}(t)$")
plt.plot(Loes2.t[N1:N],Loes2.y[3][N1:N],color="blue", alpha=0.2, label=r"$\rm r_{2x}(t), \epsilon=10^{-8}$")
plt.plot(Loes1.t[N1:N],Loes1.y[3][N1:N],color="blue", linestyle=":", label=r"$\rm r_{2x}(t)$")
plt.plot(Loes2.t[N1:N],Loes2.y[6][N1:N],color="green", alpha=0.2, label=r"$\rm r_{3x}(t), \epsilon=10^{-8}$")
plt.plot(Loes1.t[N1:N],Loes1.y[6][N1:N],color="green", linestyle=":", label=r"$\rm r_{3x}(t)$")
plt.legend(loc='upper left',fontsize=14);

Die obere Abbildung zeigt wie empfindlich das dynamische System auf kleinste Änderungen seiner Anfangsbedingungen reagiert. Die Stabilität von zeitabhängiger Trajektorien wurde bereits im Jahre 1892 von dem russischen Mathematiker Lyapunov untersucht. Eine stabile Bahn liegt im Sinne von Lyapunov dann vor, wenn ein Punkt auf einer benachbarten Trajektorie stets in der Nähe der Referenztrajektorie bleibt. Mittels des von ihm definierten Lyapunov-Exponenten $\lambda$ kann man quantifizieren, wie schnell zwei fast identische Zustände eines Systems im Laufe der Zeit auseinanderdriften.

Das Kernkonzept des Lyapunov-Exponenten gründet auf der Annahme eines näherungsweise exponentiellen Wachstums der Fehler im Laufe der Zeit. Ist $d_0$ der Fehler zur Zeit $t=0$ (Abstand zweier benachbarter Trajektorien), dann entwickelt sich die Entfernung zwischen zwei Trajektorien zum Zeitpunkt $t$ näherungsweise wie folgt $d(t) \approx d_0 \, e^{\lambda t}$. Für $\lambda \leq 0$ ist das System stabil, wohingegen für $\lambda > 0$ langfristige Vorhersagen unmöglich werden, da benachbarte Trajektorien sich exponentiell voneinander entfernen - ein eindeutiges Indiz für Chaos. Im Folgenden stellen wir den Logarithmus von $d(t)$ dar:

log($d(t)$ = log($d_0$)) + $\lambda \, t$

diff = np.abs(Loes1.y[3] - Loes2.y[3])
log_dists = np.log(diff + 1e-20)

plt.figure(figsize=(11, 5))
plt.xlabel("t")
plt.ylabel("log(d(t))")
plt.plot(Loes1.t[N1:N],log_dists[N1:N],color="orange")
plt.plot([Loes1.t[N1],Loes1.t[-1]],[0,0],color="black",linestyle=":");

Man erkennt, dass sich die Größe log($d(t)$, nach dem zweiten Entfernen des blauen Körpers ($t \approx 42.5$), mit positiver Steigung weiter entwickelt $\rightarrow$ positiver Lyapunov-Exponent.

Dreidimensionale Bewegungen

Die bisher betrachteten Bewegungen der drei Körper fanden alle auf der ($z=0$)-Ebene statt (planare Bewegungen). Wir werden im Folgenden diese Beschränkung aufgeben und allgemeine, dreidimensionale Bewegungen simulieren.

In den letzten Jahren wurden viele neue reguläre dreidimensionale Bewegungen im Dreikörperproblem gefunden. Im Jahre 2025 publizierten z.B. die Wissenschaftler Xiaoming Li und Shijun Liao von der Shanghai Jiao Tong Universität in China eine große Anzahl von neuen dreidimensionalen periodischen Bewegungen, von denen wir einige hier vorstellen möchten.

Literatur:

Li, Xiaoming, and Shijun Liao. "Discovery of 10,059 new three-dimensional periodic orbits of general three-body problem." arXiv preprint arXiv:2508.08568 (2025).

Mittels der folgenden Python-Funktion kann der Benutzer eine von 13 Anfangsbedingungen des oben definierten Artikels auswählen. Die Funktion gibt die Periode der Bewegung $T$, die Liste der Massen der drei Körper und die Anfangsbedingungen zurück.

# 3-D Anfangsdaten von https://github.com/sjtu-liao/three-body
set_G = 1

def u0_Article_2(num):
    m1, m2 = 1, 1
    T = 1
    
    match num:
        case 0:
            m3, z0, vx, vy, vz, T = 0.4, 2.02618227737279e-01, 2.39490733660267e-01, 1.54345430161257e-02, -1.59048053673444e-02, 1.51755245595921e+02
        case 1:
            m3, z0, vx, vy, vz, T = 0.6, 7.23794540561745e-01, 2.17453097136500e-01, -2.40235097321237e-01, 4.80721353013616e-01, 8.25135312975706e+00
        case 2:
            m3, z0, vx, vy, vz, T = 0.6, 6.73603712265433e-01, 2.66956499583551e-01, 1.45463973879742e-01, -6.42523240035411e-03, 3.73530217936097e+01
        case 3:
            m3, z0, vx, vy, vz, T = 0.9, -2.12165808428858e-02, 4.25238996598787e-01, 1.85105818590150e-01, 1.05183534476409e-02, 2.81797546926098e+01
        case 4:
            m3, z0, vx, vy, vz, T = 1.0, 2.33457651476844e-01, 2.55036783150535e-01, 1.25020205187277e-01, -1.37193142484893e-03, 7.14338842054417e+01
        case 5:
            m3, z0, vx, vy, vz, T = 1.0, 4.76878264280312e-01, 4.02136910074724e-01, 1.80356951286259e-01, 2.10445128137873e-01, 6.83162203628444e+00
        case 6:
            m3, z0, vx, vy, vz, T = 1.2, 6.46976570434943e-01, 4.15972648554844e-01, 2.72801109394611e-01, 3.34433078432372e-01, 7.86113798386169e+00
        case 7:
            m3, z0, vx, vy, vz, T = 1.7, 7.01144163926328e-01, 4.58924075805285e-01, 2.86482398733522e-01, 1.09684411424347e-03, 6.54566561497015e+01
        case 8:
            m3, z0, vx, vy, vz, T = 0.6, 6.14604358837056e-01, 1.03810482319459e-01, 7.24946332304501e-02, 5.89669051037992e-02, 7.97434883999123e+00
        case 9:
            m3, z0, vx, vy, vz, T = 1.1, 6.17172796472047e-01, 5.05639109864963e-01, 2.89239073968664e-01, 6.34784314640099e-03, 2.73815024124891e+01
        case 10:
            m3, z0, vx, vy, vz, T = 1.7, 5.68896058936806e-01, 1.36733707131425e-01, 5.42122376227224e-01, -1.44122788112534e-02, 1.19165128972095e+01
        case 11:
            m3, z0, vx, vy, vz, T = 1.2, 5.99263503038488e-01, 5.19875708236299e-01, 3.51407008027847e-01, -6.09898484609544e-03, 6.91026168710127e+01
        case 12:
            m3, z0, vx, vy, vz, T = 1.0, 1.42545600211905e-01, 3.49739525795482e-01, 6.02459695802159e-01, 5.62585581220896e-03, 4.44132581514811e+01 
        case _:
            m3, z0, vx, vy, vz, T = 0.6, 7.23794540561745e-01, 2.17453097136500e-01, -2.40235097321237e-01, 4.80721353013616e-01, 8.25135312975706e+00
    
    r_1_0 = np.array([-1, 0, 0])
    r_2_0 = np.array([1, 0, 0])
    r_3_0 = np.array([0, 0, z0])
    v_1_0 = np.array([vx, vy, vz])
    v_2_0 = np.array([vx, vy, -vz])
    v_3_0 = np.array([-(m1 + m2)*vx/m3, -(m1 + m2)*vy/m3, 0])
    
    set_m = [m1, m2, m3]
    u_0 = np.concatenate([r_1_0, r_2_0, r_3_0, v_1_0, v_2_0, v_3_0])
    return T, set_m, u_0

In der folgenden Abbildung sind die Trajektorien von sechs unterschiedlichen regulären dreidimensionalen periodischen Bewegungen dargestellt.

sims = [1, 6, 8, 9, 10, 12] 
t_span = [0,20]
tval = np.linspace(0, t_span[1], N)

loesungen = []
for i in sims:
    T, set_m, u_0 = u0_Article_2(i)
    loesungen.append(solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G, set_m), rtol=fehler, atol=fehler))

from mpl_toolkits.mplot3d import Axes3D

colors = ["red", "blue", "green"]
labels = [r"$\rm \vec{r}_1(t)$", r"$\rm \vec{r}_2(t)$", r"$\rm \vec{r}_3(t)$"]

fig, axes = plt.subplots(2, 3, figsize=(12, 9), subplot_kw={'projection': '3d'}, 
                         gridspec_kw={'hspace': -0.2, 'wspace': 0.2})

for ax, Loes in zip(axes.flat, loesungen):
    ax.view_init(azim=-75, elev=20)
    ax.set(xlabel="x-Koordinate", ylabel="y-Koordinate", zlabel="z-Koordinate")
    
    for i in range(3):
        ax.plot(Loes.y[3*i], Loes.y[3*i+1], Loes.y[3*i+2], 
                color=colors[i], label=labels[i], lw=0.5)
        ax.scatter(Loes.y[3*i][0], Loes.y[3*i+1][0], Loes.y[3*i+2][0], 
                   color=colors[i], s=20)

Die Bewegung der drei Körper im dreidimensionalen ($x,y,z$)-Konfigurationsraum ist in den Abbildungen oben-links und oben-mitte gut erkennbar. In den Abbildungen oben-rechts, unten--links und unten-mitte überlagern sich die Trajektorien des roten und blauen Körpers; sogenannte "piano-trio orbits" (ein Trio für zwei Violinen und ein Klavier). In der Abbildung unten-rechts überlagern sich schließlich alle drei Trajektorien (choreografische periodische 3D-Orbits des Dreikörperproblems).

Wir betrachten uns die Simulation der oberen, linken Abbildung genauer und stellen diese in einer Animation dar.

def create_animation_3d(Loes, R, Vr, r_S, N, frames_N=100):
    dN = int(N / (frames_N + 1))
    
    fig = plt.figure(figsize=(11, 5))
    gs = gridspec.GridSpec(2, 3, width_ratios=[2, 1, 1], hspace=0.3, wspace=0.4)
    
    ax1 = fig.add_subplot(gs[:, 0], projection='3d')
    ax2 = fig.add_subplot(gs[0, 1])
    ax3 = fig.add_subplot(gs[0, 2])
    ax4 = fig.add_subplot(gs[1, 1])
    ax5 = fig.add_subplot(gs[1, 2])
    
    colors = ["red", "blue", "green"]
    labels = [r"$\rm\vec{r}_1(t)$", r"$\rm \vec{r}_2(t)$", r"$\rm \vec{r}_3(t)$"]
    # Statische Plots
    ax1.view_init(azim=-75, elev=20)
    ax1.set(xlabel="x-Koordinate", ylabel="y-Koordinate", zlabel="z-Koordinate")
    for i in range(3):
        ax1.plot(Loes.y[3*i], Loes.y[3*i+1], Loes.y[3*i+2], 
                color=colors[i], label=labels[i], lw=0.5)
    ax1.legend(loc='upper left')
    
    ax2.set_xlabel(r"Radius $r_1$")
    ax2.set_ylabel(r"radiale Geschwind. $v_{r_1}$")
    ax2.plot(R[0], Vr[0], color="red")
    
    ax3.set_xlabel(r"Radius $r_2$")
    ax3.set_ylabel(r"radiale Geschwind. $v_{r_2}$")
    ax3.plot(R[1], Vr[1], color="blue")
    
    ax4.set_xlabel(r"Radius $r_3$")
    ax4.set_ylabel(r"radiale Geschwind. $v_{r_3}$")
    ax4.plot(R[2], Vr[2], color="green")
    
    ax5.set_xlabel("x-Koordinate")
    ax5.set_ylabel("y-Koordinate")
    ax5.plot(r_S[0], r_S[1], color="black", label=r"$\rm Schwerpunkt$")
    ax5.legend(loc='upper left')
    
    # Animationsobjekte
    dot_s = 30
    dot_r1 = ax1.scatter([], [], [], color='red', s=dot_s)
    dot_r2 = ax1.scatter([], [], [], color='blue', s=dot_s)
    dot_r3 = ax1.scatter([], [], [], color='green', s=dot_s)
    dot_r1a = ax2.scatter([], [], color='red', s=dot_s)
    dot_r2a = ax3.scatter([], [], color='blue', s=dot_s)
    dot_r3a = ax4.scatter([], [], color='green', s=dot_s)
    
    line_S1, = ax1.plot([], [], [], color='grey', alpha=0.2)
    line_S2, = ax1.plot([], [], [], color='grey', alpha=0.2)
    line_S3, = ax1.plot([], [], [], color='grey', alpha=0.2)
    
    def update(frame):
        i = frame * dN + 1
        dot_r1._offsets3d = ([Loes.y[0][i]], [Loes.y[1][i]], [Loes.y[2][i]])
        dot_r2._offsets3d = ([Loes.y[3][i]], [Loes.y[4][i]], [Loes.y[5][i]])
        dot_r3._offsets3d = ([Loes.y[6][i]], [Loes.y[7][i]], [Loes.y[8][i]])
        
        dot_r1a.set_offsets([[R[0][i], Vr[0][i]]])
        dot_r2a.set_offsets([[R[1][i], Vr[1][i]]])
        dot_r3a.set_offsets([[R[2][i], Vr[2][i]]])
        
        line_S1.set_data([r_S[0][i], Loes.y[0][i]], [r_S[1][i], Loes.y[1][i]])
        line_S1.set_3d_properties([r_S[2][i], Loes.y[2][i]])        
        line_S2.set_data([r_S[0][i], Loes.y[3][i]], [r_S[1][i], Loes.y[4][i]])
        line_S2.set_3d_properties([r_S[2][i], Loes.y[5][i]])
        line_S3.set_data([r_S[0][i], Loes.y[6][i]], [r_S[1][i], Loes.y[7][i]])
        line_S3.set_3d_properties([r_S[2][i], Loes.y[8][i]])
        
        return dot_r1, dot_r2, dot_r3, dot_r1a, dot_r2a, dot_r3a, line_S1, line_S2, line_S3

    ani = FuncAnimation(fig, update, frames=frames_N, interval=100, blit=True)
    plt.close(fig) 
    return ani

def create_animation_3d_short(Loes, R, Vr, r_S, N, frames_N=100):
    # Berechnung der zu Animation verwendeten Indices
    indices = [int(f * (N - 1) / (frames_N - 1)) for f in range(frames_N)]
    
    fig = plt.figure(figsize=(11, 5))
    gs = gridspec.GridSpec(2, 3, width_ratios=[2, 1, 1], hspace=0.3, wspace=0.4)
    
    # Achsen-Setup
    ax_main = fig.add_subplot(gs[:, 0], projection='3d')
    axes_sub = [fig.add_subplot(gs[0, 1]), fig.add_subplot(gs[0, 2]), fig.add_subplot(gs[1, 1])]
    ax_cm = fig.add_subplot(gs[1, 2])

    colors = ['red', 'blue', 'green']
    labels = [r"$\rm \vec{r}_1(t)$", r"$\rm \vec{r}_2(t)$", r"$\rm \vec{r}_3(t)$"]    
    dots_main, dots_sub, lines_cm = [], [], []

    # Hauptplot (Trajektorien im Konfigurationsraum (x,y,z))
    ax_main.view_init(azim=-75, elev=20)
    for j in range(3):
        idx_x, idx_y, idx_z = j * 3, j * 3 + 1, j * 3 + 2
        ax_main.plot(Loes.y[idx_x], Loes.y[idx_y], Loes.y[idx_z], 
                     color=colors[j], label=labels[j], lw=0.6, alpha=0.8)
        dots_main.append(ax_main.scatter([], [], [], color=colors[j], s=20))
        lines_cm.append(ax_main.plot([], [], [], color='grey', alpha=0.2)[0])
    ax_main.legend(loc='upper right')

    # Phasenraum-Plots (R vs Vr)
    for j, ax in enumerate(axes_sub):
        ax.plot(R[j], Vr[j], color=colors[j], alpha=0.3)
        dots_sub.append(ax.scatter([], [], color=colors[j], s=25))
        ax.set_xlabel(rf"$r_{j+1}$")
        ax.set_ylabel(rf"$v_{{r_{j+1}}}$")

    # Schwerpunktplot
    ax_cm.plot(r_S[0], r_S[1], color="black", label="Schwerpunkt")
    ax_cm.legend(loc='upper left')

    def update(frame):
        i = indices[frame]
        changed_objects = []
        
        for j in range(3):
            # Update Punkte der Körper im Hauptplot
            x, y, z = Loes.y[j*3][i], Loes.y[j*3+1][i], Loes.y[j*3+2][i]
            dots_main[j]._offsets3d = ([x],[y], [z])
            
            # Update Verbindungslinien zum Schwerpunkt (Hauptplot)
            lines_cm[j].set_data([r_S[0][i], x], [r_S[1][i], y])
            lines_cm[j].set_3d_properties([r_S[2][i], z])
            
            # Update Punkte in den Phasenraum-Plots
            dots_sub[j].set_offsets([[R[j][i], Vr[j][i]]])
            
            changed_objects.extend([dots_main[j], lines_cm[j], dots_sub[j]])
            
        return changed_objects

    ani = FuncAnimation(fig, update, frames=frames_N, interval=120, blit=True)
    plt.close(fig)
    return ani

fehler = 1e-9
N = 10000
T, set_m, u_0 = u0_Article_2(1)
t_span = [0,T]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_3d_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Wir betrachten uns auch die Simulation der oberen, rechten Abbildung genauer.

T, set_m, u_0 = u0_Article_2(8)
t_span = [0,T]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_3d_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Man erkennt im Hauptplot, dass die Trajektorien des roten und blauen Körpers im Konfigurationsraum (x,y,z) sich überlagern ("piano-trio orbit").

Nun betrachten wir uns die Simulation der unteren, rechten Abbildung genauer.

T, set_m, u_0 = u0_Article_2(12)
t_span = [0,T]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)
r_S, R, Vr = Transf_S(Loes,set_m)
ani = create_animation_3d_short(Loes, R, Vr, r_S, N)
HTML(ani.to_html5_video())

Man erkennt, dass sich bei dieser Simulation alle Trajektorien der Körper im Konfigurationsraum (x,y,z) überlagern ("choreographischer periodischer 3D-Orbit").

Im folgenden Unterpunkt werden wir nun das N-Körperproblem der Planetenbewegungen unseres Sonnensystems und die Sonnenfinsternis 2026/27 behandeln.

Planetenbewegungen und das n-Körper-Problem

Um die Bewegung der Planeten in unserem Sonnensystem mit absoluter Präzision vorherzusagen, muss man alle Himmelskörper gleichzeitig als n-Körper-Problem behandeln. Am Ende dieses Jupyter-Notebooks werden wir eine solche Simulation auch durchführen. Zunächst betrachten wir jedoch die drei Himmelskörper, Sonne-Erde-Mond, als ein klassisches Dreikörper-Problem und vernachlässigen dadurch z.B. den Einfluss von Jupiter und Venus, die messbare, aber langsame Veränderungen der Erdbahn verursachen.

Wir verwenden nun den folgenden Wert für die Gravitationskonstante $G = 6.6743015 \cdot 10^{-20} \frac{\rm km^3}{\rm kg \, s^2}$ und setzen die Massen der drei Körper auf: Masse der Sonne: $M_\odot =1.9884 \cdot 10^{30} {\rm kg}$, Masse der Erde: $M_{\!E} = 5.9722 \cdot 10^{24} {\rm kg}$ und Masse des Mondes: $M_{\!M} = 7.346 \cdot 10^{22} {\rm kg}$

set_G = 6.6743015*10**(-20)

m_S = 1.9884*10**30
m_E = 5.9722*10**24
M_M = 7.346*10**22
set_m = [m_S, m_E, M_M]

Die Anfangswerte der Himmelskörper können wir uns mittels der Python-Bibliothek Skyfield generieren. Dies ist ein Python-Modul für astronomische Berechnungen, mit dem man hochpräzise Positionen und Geschwindigkeiten für Himmelskörper wie Planeten, Sterne und Erdsatelliten erzeugen kann. Das Modul basiert auf der Ephemeriden-Datei 'de421.bsp', die von der NASA entwickelt wurde.

Wir simulieren nun die Bewegung der drei Himmelskörper für ca. ein Jahr (365 Tage) und wählen als Anfangszeitpunkt das heutige Datum. Um später unsere Simulation mit den wirklichen Werten vergleichen zu können, laden wir uns die Positionen und Geschwindigkeiten für ein gesamtes Jahr in $N$-äquidistanten Zeitschritten.

from skyfield.api import load

# Lade Ephemeriden-Datei
eph = load('de421.bsp')
ts = load.timescale()

# Beispiel: 14. April 2026 um 15:30 Uhr UTC
#zeitpunkt = ts.utc(2026, 4, 14, 15, 30)
start_time = ts.now()
N = 10001
days = np.linspace(0, 365, N) 
times = start_time + days

r_1 = eph['sun'].at(times).position.km
v_1 = eph['sun'].at(times).velocity.km_per_s
r_2 = eph['earth'].at(times).position.km
v_2 = eph['earth'].at(times).velocity.km_per_s
r_3 = eph['moon'].at(times).position.km
v_3 = eph['moon'].at(times).velocity.km_per_s

u_0 = np.concatenate([r_1[:,0], r_2[:,0], r_3[:,0], v_1[:,0], v_2[:,0], v_3[:,0]])

# Beispiel-Ausgabe der Anfangsdaten der Erde
print(f"Anfangszeitpunkt: {start_time.utc_strftime('%Y-%m-%d %H:%M UTC')}")
print(f"Erdposition (km): {r_2[:,0]}")
print(f"Erdgeschwindigkeit (km/s): {v_2[:,0]}")

Anfangszeitpunkt: 2026-04-04 14:46 UTC
Erdposition (km): [-1.45225990e+08 -3.50768580e+07 -1.51845824e+07]
Erdgeschwindigkeit (km/s): [  6.96779692 -26.55583522 -11.51125392]

Wir simulieren nun die Bewegung der drei Himmelskörper für ein Jahr. Zunächst benutzen wir dabei die Einheiten Kilometer und Sekunde. Wir werden jedoch bei den späteren numerischen Berechnungen geeignetere Einheiten wählen.

fehler = 1e-9
t_span = [0,365*24*60*60]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys, t_span, u_0, t_eval=tval, args=(set_G,set_m), rtol=fehler, atol=fehler)

Wir stellen uns im Folgenden die zeitliche Entwicklung der x- und y-Koordinate der Erde im Vergleich zu den wirklichen Skyfield-Daten dar.

plt.figure(figsize=(11, 4))
plt.xlabel("t")
plt.ylabel("x-, y-Koordinate")
plt.plot(times.utc_datetime(),Loes.y[3],color="blue", label=r"$\rm x_{E}(t)=x_{2}(t), Simulation$", alpha = 0.8)
plt.plot(times.utc_datetime(),r_2[0,:],color="black", linestyle = ":", label=r"$\rm x_{E}(t)=x_{2}(t), \, de421-Daten$")
plt.plot(times.utc_datetime(),Loes.y[4],color="blue", label=r"$\rm y_{E}(t)=y_{2}(t), Simulation$", alpha = 0.2)
plt.plot(times.utc_datetime(),r_2[1,:],color="black", linestyle = ":", label=r"$\rm y_{E}(t)=y_{2}(t), \, de421-Daten$")
plt.grid(True)
plt.legend(loc='upper left',fontsize=11);

Die obere Abbildung zeigt, dass die Simulation gut mit den Skyfield-Daten übereinstimmt. Bei dem genauen Vergleich der Daten (siehe folgende print-Ausgabe) erkennt man jedoch Abweichungen; dies werden wir später diskutieren.

print(f"Endzeitpunkt: {times[-1].utc_strftime('%Y-%m-%d %H:%M UTC')}")
print(f"Erdposition (km): {r_2[:,-1]}")
print(f"Erdgeschwindigkeit (km/s): {v_2[:,-1]}")
print("Simulation:")
print(f"Erdposition (km): {Loes.y[3:6,-1]}")
print(f"Erdgeschwindigkeit (km/s): {Loes.y[12:15,-1]}")

Endzeitpunkt: 2027-04-04 14:46 UTC
Erdposition (km): [-1.45054902e+08 -3.43713441e+07 -1.48868301e+07]
Erdgeschwindigkeit (km/s): [  6.84425812 -26.59988896 -11.53174225]
Simulation:
Erdposition (km): [-1.45005519e+08 -3.44449112e+07 -1.49196515e+07]
Erdgeschwindigkeit (km/s): [  6.84842988 -26.60369689 -11.53354938]

Bevor wir die Trajektorien der Bewegungen in drei Dimensionen darstellen, wollen wir uns zuächst die kommenden Sonnenfinsternisse berechnen. Eine Sonnenfinsternis findet statt, wenn der Richtungsvektor Erde-Sonne ($\vec{w}_1 := \vec{r}_S - \vec{r}_E$) annähernd gleich dem Richtungsvektor Erde-Mond ($\vec{w}_2 :=\vec{r}_M - \vec{r}_E$) ist (Sonne, Mond und Erde müssen nahezu auf einer Geraden liegen). Der Winkel $\alpha$ zwischen diesen Vektoren ist dann sehr klein, wobei $\alpha = \hbox{arcos} \left( \frac{\vec{w}_1 \cdot \vec{w}_2}{\left| \vec{w}_1 \right| \, \left| \vec{w}_2 \right|} \right)$. Damit eine Sonnenfinsternis für einen Beobachter auf der Erde stattfindet, muss der Winkel $\alpha$ kleiner sein als die Summe der scheinbaren Radien von Sonne und Mond. Mittels der simulierten Daten berechnen wir uns im Folgenden die zeitliche Entwicklung des Winkels ($\alpha(t)$).

def winkel_sonne_mond(pos_sun, pos_earth, pos_moon):
    vec_earth_to_sun  = pos_sun - pos_earth
    vec_earth_to_moon = pos_moon - pos_earth
    
    # Einheitsvektoren 
    v1_e = vec_earth_to_sun / np.linalg.norm(vec_earth_to_sun)
    v2_e = vec_earth_to_moon / np.linalg.norm(vec_earth_to_moon)
    
    # Skalarprodukt
    dot_product = np.dot(v1_e, v2_e)
    
    # Winkel Sonne-Mond in Radiant
    winkel_rad = np.arccos(dot_product)
    
    # Rückgabewert in Grad 
    return np.degrees(winkel_rad)

winkel = []
for i in range(N):
    pos_sun   = Loes.y[0:3,i]
    pos_earth = Loes.y[3:6,i]
    pos_moon  = Loes.y[6:9,i]
    winkel.append(winkel_sonne_mond(pos_sun, pos_earth, pos_moon))

Eine kürzere Variante der Winkelberechnung stellt die untere Python-Funktion dar:

def winkel_sonne_mond_short(pos_sun, pos_earth, pos_moon):
    earth_to_sun = pos_sun - pos_earth
    earth_to_moon = pos_moon - pos_earth

    # Einheitsvektoren
    v1_e = earth_to_sun / np.linalg.norm(earth_to_sun, axis=0)
    v2_e = earth_to_moon / np.linalg.norm(earth_to_moon, axis=0)

    # Skalarprodukt für jeden Zeitschritt 
    dot_product = np.sum(v1_e * v2_e, axis=0)

    # Winkel Sonne-Mond in Radiant
    winkel_rad = np.arccos(dot_product)
    
    # Rückgabewert in Grad 
    return np.degrees(winkel_rad)

winkel = winkel_sonne_mond_short(Loes.y[0:3], Loes.y[3:6], Loes.y[6:9])

Die untere Abbildung zeigt die zeitliche Entwicklung des Winkels $\alpha(t)$ im betrachteten Zeitraum.

F_1 = ts.utc(2026, 8, 12, 20, 30)
F_2 = ts.utc(2027, 2, 6, 12, 0)
plt.figure(figsize=(11, 4))
plt.xlabel("t")
plt.ylabel(r"Winkel $\alpha(t)$")
plt.plot(times.utc_datetime(),winkel, color="orange")
plt.scatter(F_1.utc_datetime(),0.1, color="black", s=30, label="Totale Finsternis am 12. August 2026")
plt.scatter(F_2.utc_datetime(),0.1, color="brown", s=30, label="Partielle Finsternis am 6. Februar 2027")
plt.grid(True)
plt.ylim(0, 7.5)
plt.legend(loc='lower left',fontsize=10);

In der oberen Abbildung sind auch die zwei kommenden Sonnenfinsternisse markiert.

• Am 12. August 2026 findet eine totale Sonnenfinsternis statt, wobei der Kernschatten über Grönland, Island und Nordspanien (u. a. Mallorca und Ibiza) verläuft.
• Am 6. Februar 2027 findet eine ringförmige Sonnenfinsternis statt, die unter anderem in Argentinien und Brasilien zu sehen sein wird.

Die simulierten Berechnungen des Winkel $\alpha$ stimmen dabei sehr gut mit den Daten der Sonnenfinsternisse überein. Gerade an den Zeitpunkten der Sonnenfinsternisse ist der Winkel $\alpha$ sehr klein.

Wir stellen uns nun die Trajektorien der Bewegungen in drei Dimensionen dar. Da man aufgrund der sehr unterschiedlichen Abstandskalen die Bewegung des Systems Sonne-Erde-Mond nicht adäquat in einer Abbildung darstellen kann, wurden die Systeme Sonne-Erde und Erde-Mond in zwei getrennten Abbildungen dargestellt.

def create_animation_eclipse(Loes, N, frames_N=100):
    dN = int(N / (frames_N + 1))
    dot_s = 30
    fac = 0.005
    
    fig = plt.figure(figsize=(10, 6))
    gs = gridspec.GridSpec(2, 2, width_ratios=[1,1], height_ratios=[3,1], hspace=0.3, wspace=0)
    
    ax1 = fig.add_subplot(gs[0, 0], projection='3d')
    ax2 = fig.add_subplot(gs[0, 1], projection='3d')
    ax3 = fig.add_subplot(gs[1, :])
    
    colors = ["red", "blue", "green"]
    labels = [r"$\rm \vec{r}_1(t)$", r"$\rm \vec{r}_2(t)$", r"$\rm \vec{r}_3(t)$"]
    labels = [r"$\rm \vec{r}_{sun}$", r"$\rm \vec{r}_{earth}$", r"$\rm \vec{r}_{moon}$"]
    
    # Statische Plots
    ax1.view_init(azim=10, elev=25)
    ax1.set(xlabel="x[km]", ylabel="y[km]", zlabel="z[km]")
    for i in [0,1]:
        ax1.plot(Loes.y[3*i], Loes.y[3*i+1], Loes.y[3*i+2], 
                color=colors[i], label=labels[i], lw=0.5)
    ax1.legend(loc='upper left')
    
    ax2.view_init(azim=10, elev=25)
    ax2.set(xlabel="x", ylabel="y", zlabel="z")
    t_end = int(N/12)
    for i in [1,2]:
        ax2.plot(Loes.y[3*i][0:t_end]-Loes.y[3][0:t_end], Loes.y[3*i+1][0:t_end]-Loes.y[4][0:t_end], 
                 Loes.y[3*i+2][0:t_end]-Loes.y[5][0:t_end], color=colors[i], label=labels[i], lw=0.5)
    ax2.scatter(0, 0, 0, color='blue', s=dot_s)
    ax2.legend(loc='upper left')
    
    ax3.set_xlabel("t")
    ax3.set_ylabel("Winkel")
    ax3.plot(times.utc_datetime(),winkel, color="orange")
    ax3.scatter(F_1.utc_datetime(),0.1, color="black", s=10)
    ax3.scatter(F_2.utc_datetime(),0.1, color="brown", s=10)
    ax3.set_ylim(0, 7.5)
    
    
    # Animationsobjekte
    dot_r1 = ax1.scatter([], [], [], color='red', s=dot_s)
    dot_r2 = ax1.scatter([], [], [], color='blue', s=dot_s)
    line_S1, = ax1.plot([], [], [], color='red', alpha=0.2)
    
    dot_r3 = ax2.scatter([], [], [], color='green', s=dot_s)
    line_S2, = ax2.plot([], [], [], color='blue', alpha=0.2)
    line_S3, = ax2.plot([], [], [], color='red', alpha=0.2)
    
    hline, = ax3.plot([], [], color='grey', linewidth=3, alpha=0.8)
    
    def update(frame):
        i = frame * dN + 1
        time = times.utc_datetime()[i]
        dot_r1._offsets3d = ([Loes.y[0][i]], [Loes.y[1][i]], [Loes.y[2][i]])
        dot_r2._offsets3d = ([Loes.y[3][i]], [Loes.y[4][i]], [Loes.y[5][i]])
        line_S1.set_data([Loes.y[3][i],Loes.y[0][i]], [Loes.y[4][i],Loes.y[1][i]])
        line_S1.set_3d_properties([Loes.y[5][i],Loes.y[2][i]])

        dot_r3._offsets3d = ([Loes.y[6][i]-Loes.y[3][i]], [Loes.y[7][i]-Loes.y[4][i]], [Loes.y[8][i]-Loes.y[5][i]])
        line_S2.set_data([0,Loes.y[6][i]-Loes.y[3][i]], [0,Loes.y[7][i]-Loes.y[4][i]])
        line_S2.set_3d_properties([0,Loes.y[8][i]-Loes.y[5][i]])
        line_S3.set_data([0,fac*(Loes.y[0][i]-Loes.y[3][i])], [0,fac*(Loes.y[1][i]-Loes.y[4][i])])
        line_S3.set_3d_properties([0,fac*(Loes.y[2][i]-Loes.y[5][i])])

        hline.set_data([time, time], [0, 7.5])
        
        return dot_r1, dot_r2, dot_r3, line_S1, line_S2, line_S3, hline

    ani = FuncAnimation(fig, update, frames=frames_N, interval=100, blit=True)
    plt.close(fig) 
    return ani

ani = create_animation_eclipse(Loes, N)
HTML(ani.to_html5_video())

Die obere linke Animation stellt die Bewegung der Erde um die Sonne dar und die obere rechte Abbildung zeigt die Bewegung des Mondes um die Erde, wobei die rote Linie den Verbindungsvektor Erde-Sonne und die blaue Linie den Verbindungsvektor Erde-Mond skizziert. Die untere Abbildung zeigt die zeitliche Entwicklung des Winkel $\alpha(t)$.

Wir möchten uns nun auch die Bewegung der anderen Planeten visualisieren, sodass wir nun ein n-Körperproblem zu lösen haben, wobei n die Zahl der simulierten Planeten darstellt. Die zugrundeliegende Bewegungsgleichung wird nun durch das folgende System von gekoppelten Differenzialgleichungen zweiter Ordnung beschrieben:

\begin{equation} \frac{d^2 \vec{r}_i}{dt^2} = \sum_{j=1\\i \neq j}^n G \cdot m_j \cdot \frac{\vec{r}_j-\vec{r}_i}{\left| \, \vec{r}_j-\vec{r}_i \right|^3} \nonumber \end{equation}

wobei $\vec{r}_i$ und $m_i, i \in {1,2,...,n}$ der Ortsvektor und die Masse des Körpers $i$ und $G$ die Gravitationskonstante ist.

Betrachtet man die Bewegung der n-Körper in allen drei Raumdimensionen ($\vec{r_i}=\left( x_{i}, y_{i}, z_{i} \right)$), so hat gelangt man zu 3n gekoppelten Differenzialgleichungen zweiter Ordnung, die man dann in ein System von 6n gekoppelten DGLs erster Ordnung umschreiben muss um sie numerisch lösen zu können.

# System der Differenzialgleichungen des n-Körper-Problems
# G_m[i]: Masse des Körpers i mal Gravitationskonstanten G
def DGLsys_n(t,u_vec,G_m,n):
    dim = int(3*n)
    # Entpacke u_vec in dim=3*n Koordinaten und 3*n Geschwindigkeiten
    r = u_vec[:dim].reshape((n, 3))
    v = u_vec[dim:].reshape((n, 3))
    # DGL-System erster Ordnung (3*n*2 Gleichungen)
    drdt = v
    dvdt = np.zeros_like(v)
    for i in range(n):
        for j in range(n):
            if i != j:
                dvdt[i] += G_m[j] * (r[j] - r[i]) / np.linalg.norm(r[j] - r[i])**3            
    return np.concatenate([drdt.flatten(), dvdt.flatten()])

Die obere Python-Funktion der Differenzialgleichung benutzt die kombinierten Massen-Gravitationsparameter $\mu_j = G \cdot m_j$ der Himmelskörper. Die Zahlenwerte dieser Parameter für die einzelnen Himmelskörper sind in dem unteren Python-Dictionary definiert.

# Gravitationsparameter GM in [AU^3 / d^2]
# Basierend auf JPL-Daten (DE421/DE440)
# Farben für Planeten-Plots, Sonnengelb, Merkur-Grau, Venus-Beige (Atmosphäre), Erd-Blau, Mond-Hellgrau, 
# Mars-Rostrot, Jupiter-Bänder-Braun, Saturn-Goldgelb, Uranus-Hellblau/Cyan, Neptun-Tiefblau, Pluto-Graubraun
bodies_data = {
    'sun': {'gm': 2.959122082855911e-04, 'color': '#FFCC00', 'label': r"$\rm \vec{r}_{sun}$"},
    'mercury': {'gm': 4.91248045036476e-11, 'color': '#A5A5A5', 'label': r"$\rm \vec{r}_{mercury}$"},
    'venus': {'gm': 7.243452486161142e-10, 'color': '#E3BB76', 'label': r"$\rm \vec{r}_{venus}$"},
    'earth': {'gm': 8.887692390113509e-10, 'color': '#2E5FAC', 'label': r"$\rm \vec{r}_{earth}$"},
    'moon': {'gm': 1.093189450742374e-11, 'color': '#D3D3D3', 'label': r"$\rm \vec{r}_{moon}$"},
    'mars': {'gm': 9.549535105170911e-11, 'color': '#E27B58', 'label': r"$\rm \vec{r}_{mars}$"},
    'jupiter barycenter': {'gm': 2.825345909524226e-07, 'color': '#D39C7E', 'label': r"$\rm \vec{r}_{jupiter}$"},
    'saturn barycenter': {'gm': 8.459706073245032e-08, 'color': '#C5AB6E', 'label': r"$\rm \vec{r}_{saturn}$"},
    'uranus barycenter': {'gm': 1.292024825783306e-08, 'color': '#B5E2E2', 'label': r"$\rm \vec{r}_{uranus}$"},
    'neptune barycenter': {'gm': 1.524357347885194e-08, 'color': '#4B70DD', 'label': r"$\rm \vec{r}_{neptune}$"},
    'pluto barycenter': {'gm': 2.17844105333832e-12, 'color': '#96745E', 'label': r"$\rm \vec{r}_{pluto}$"}
}

Die Anfangswerte der Himmelskörper können wir uns wieder mittels der Python-Bibliothek Skyfield generieren. Wir extrahieren die Daten nun in den Einheiten [AU] (Astronomische Einheit, 1 AU = 149597870700 m) und Tagen (d) und fügen sie ebenfalls in dem oben definierten Python-Dictionary. Zum späteren Vergleich mit den simulierten Werten, speichern wir uns die Positions- und Geschwindigkeitswerte der gesamten, simulierten Zeitreihe.

# Daten der Positionen und Geschwindigkeiten für alle Zeiten in "times" abrufen und in bodies_data speichern
for name in bodies_data.keys():
    bodies_data[name]['r'] = eph[name].at(times).position.au 
    bodies_data[name]['v'] = eph[name].at(times).velocity.au_per_d 

Wir simulieren nun zunächst die Bewegung der inneren Planeten (Merkur, Venus, Erde und Mars) und des Mondes.

set_bodies = ['sun','mercury','venus','earth','moon','mars']
set_n = len(set_bodies)

set_m = []
r_0 = []
v_0 = []
for name in set_bodies:
    set_m.append(bodies_data[name]['gm'])
    r_0.append(bodies_data[name]['r'][:, 0])
    v_0.append(bodies_data[name]['v'][:, 0])

u_0 = np.concatenate([np.concatenate(r_0), np.concatenate(v_0)])

fehler = 1e-9
t_span = [0,365]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys_n, t_span, u_0, t_eval=tval, args=(set_m,set_n), rtol=fehler, atol=fehler)

Wir stellen die Bewegung wieder in einer Animation dar.

def create_animation_planets_a(Loes, N, bodies, frames_N=100):
    dN = int(N / (frames_N + 1))
    dot_s = 30
    fac = 0.005
    n = len(bodies)
    
    fig = plt.figure(figsize=(10, 6))
    gs = gridspec.GridSpec(2, 2, width_ratios=[1,1], height_ratios=[3,1], hspace=0.3, wspace=0)
    
    ax1 = fig.add_subplot(gs[0, 0], projection='3d')
    ax2 = fig.add_subplot(gs[0, 1], projection='3d')
    ax3 = fig.add_subplot(gs[1, :])
    
    
    i_moon = bodies.index('moon')
    i_earth = bodies.index('earth')
    # Statische Plots
    ax1.view_init(azim=10, elev=25)
    ax1.set(xlabel="x[AU]", ylabel="y[AU]", zlabel="z[AU]")
    for i in range(n):
        if i != i_moon:
            ax1.plot(Loes.y[3*i], Loes.y[3*i+1], Loes.y[3*i+2], 
                color=bodies_data[bodies[i]]['color'], lw=0.5, label=bodies_data[bodies[i]]['label'])
    ax1.legend(bbox_to_anchor=[0, 0.9])
    
    ax2.view_init(azim=10, elev=25)
    ax2.set(xlabel="x", ylabel="y", zlabel="z")
    t_end = int(N/12)
    for i in range(n):
        if i == i_moon:
            ax2.plot(Loes.y[3*i_moon][0:t_end]-Loes.y[3*i_earth][0:t_end], Loes.y[3*i_moon+1][0:t_end]-Loes.y[3*i_earth+1][0:t_end], 
                 Loes.y[3*i_moon+2][0:t_end]-Loes.y[3*i_earth+2][0:t_end], color=bodies_data[bodies[i]]['color'], lw=0.5
                    , label=bodies_data[bodies[i]]['label'])
    ax2.scatter(0, 0, 0, color=bodies_data['earth']['color'], s=dot_s, label=bodies_data['earth']['label'])
    ax2.legend(bbox_to_anchor=[0, 0.7])
    
    ax3.set_xlabel("t")
    ax3.set_ylabel("Winkel")
    ax3.plot(times.utc_datetime(),winkel, color="orange")
    ax3.scatter(F_1.utc_datetime(),0.1, color="black", s=10)
    ax3.scatter(F_2.utc_datetime(),0.1, color="brown", s=10)
    ax3.set_ylim(0, 7.5)
    
    
    # Animationsobjekte
    dots_ax1 = []
    for j, body in enumerate(bodies):
        if body != 'moon':
            dot = ax1.scatter([], [], [], color=bodies_data[body]['color'], s=dot_s)
            dots_ax1.append((j, dot)) 
    
    dot_moon = ax2.scatter([], [], [], color=bodies_data['moon']['color'], s=dot_s)
    line_S2, = ax2.plot([], [], [], color=bodies_data['earth']['color'], alpha=0.7, lw=0.7)
    line_S3, = ax2.plot([], [], [], color=bodies_data['sun']['color'], alpha=0.7, lw=0.7)
    
    hline, = ax3.plot([], [], color='grey', linewidth=3, alpha=0.8)
    
    def update(frame):
        i = frame * dN + 1
        time = times.utc_datetime()[i]
        
        changed_objects = []
        for j, dot in dots_ax1:
            dot._offsets3d = ([Loes.y[3*j][i]], [Loes.y[3*j+1][i]], [Loes.y[3*j+2][i]])
            changed_objects.append(dot)

        dot_moon._offsets3d = ([Loes.y[3*i_moon][i]-Loes.y[3*i_earth][i]], 
                               [Loes.y[3*i_moon+1][i]-Loes.y[3*i_earth+1][i]], 
                               [Loes.y[3*i_moon+2][i]-Loes.y[3*i_earth+2][i]])
        changed_objects.append(dot_moon)
        line_S2.set_data([0,Loes.y[3*i_moon][i]-Loes.y[3*i_earth][i]], 
                         [0,Loes.y[3*i_moon+1][i]-Loes.y[3*i_earth+1][i]])
        line_S2.set_3d_properties([0,Loes.y[3*i_moon+2][i]-Loes.y[3*i_earth+2][i]])
        changed_objects.append(line_S2)
        line_S3.set_data([0,fac*(Loes.y[0][i]-Loes.y[3*i_earth][i])], 
                         [0,fac*(Loes.y[1][i]-Loes.y[3*i_earth+1][i])])
        line_S3.set_3d_properties([0,fac*(Loes.y[2][i]-Loes.y[3*i_earth+2][i])])
        changed_objects.append(line_S3)

        hline.set_data([time, time], [0, 7.5])
        changed_objects.append(hline)
        
        return tuple(changed_objects)

    ani = FuncAnimation(fig, update, frames=frames_N, interval=100, blit=True)
    plt.close(fig) 
    return ani

ani = create_animation_planets_a(Loes, N, set_bodies)
HTML(ani.to_html5_video())

Wir stellen nun auch die anderen Planeten in einer Animation dar.

set_bodies = ['sun','mercury', 'venus','earth','moon','mars','jupiter barycenter', 
              'saturn barycenter', 'uranus barycenter', 'neptune barycenter', 'pluto barycenter']
set_n = len(set_bodies)

set_m = []
r_0 = []
v_0 = []
for name in set_bodies:
    set_m.append(bodies_data[name]['gm'])
    r_0.append(bodies_data[name]['r'][:, 0])
    v_0.append(bodies_data[name]['v'][:, 0])

u_0 = np.concatenate([np.concatenate(r_0), np.concatenate(v_0)])

Die Animationsdauer umspannt wieder ein Jahr, wobei wir die Trajektorien der Planeten jedoch für 50 Jahre darstellen, um auch die Bewegung der langsamen äußeren Planeten visualisieren zu können.

fehler = 1e-9
t_span = [0,365]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys_n, t_span, u_0, t_eval=tval, args=(set_m,set_n), rtol=fehler, atol=fehler)

t_span_l = [0,50*365]
tval_l = np.linspace(0, t_span_l[1], N)
Loes_l = solve_ivp(DGLsys_n, t_span_l, u_0, t_eval=tval_l, args=(set_m,set_n), rtol=fehler, atol=fehler)

def create_animation_planets_b(Loes, Loes_L, N, bodies, frames_N=100):
    # Berechnung der zu Animation verwendeten Indices
    indices = [int(f * (N - 1) / (frames_N - 1)) for f in range(frames_N)]
    
    dot_s = 5
    fac = 0.005
    n = len(bodies)
    
    fig = plt.figure(figsize=(10, 6))
    gs = gridspec.GridSpec(2, 2, width_ratios=[1,1], height_ratios=[3,1], hspace=0.3, wspace=0)
    
    ax1 = fig.add_subplot(gs[0, 0], projection='3d')
    ax2 = fig.add_subplot(gs[0, 1], projection='3d')
    ax3 = fig.add_subplot(gs[1, :])
    
    set_bodies_small = ['sun','mercury', 'venus','earth','mars']
    i_moon = bodies.index('moon')
    
    # Statische Plots
    ax1.view_init(azim=10, elev=40)
    ax1.set(xlabel="x", ylabel="y", zlabel="z")
    for i in range(n):
        if i != i_moon:
            ax1.plot(Loes_L.y[3*i], Loes_L.y[3*i+1], Loes_L.y[3*i+2], 
                color=bodies_data[bodies[i]]['color'], 
                     label=bodies_data[bodies[i]]['label'], lw=0.5)
    ax1.legend(bbox_to_anchor=[0, 0.9],fontsize=10)
    
    ax2.view_init(azim=10, elev=40)
    ax2.set(xlabel="x", ylabel="y", zlabel="z")
    for i in range(n):
        if bodies[i] in set_bodies_small:
            ax2.plot(Loes_L.y[3*i], Loes_L.y[3*i+1], Loes_L.y[3*i+2], 
                color=bodies_data[bodies[i]]['color'],
                     label=bodies_data[bodies[i]]['label'], lw=0.3)
    ax2.legend(bbox_to_anchor=[0, 0.7],fontsize=8)
    
    ax3.set_xlabel("t")
    ax3.set_ylabel("Winkel")
    ax3.plot(times.utc_datetime(),winkel, color="orange")
    ax3.scatter(F_1.utc_datetime(),0.1, color="black", s=10)
    ax3.scatter(F_2.utc_datetime(),0.1, color="brown", s=10)
    ax3.set_ylim(0, 7.5)
    
    # Animationsobjekte
    dots_ax1 = []
    dots_ax2 = []
    for j, body in enumerate(bodies):
        if body != 'moon':
            dot = ax1.scatter([], [], [], color=bodies_data[body]['color'], s=dot_s)
            dots_ax1.append((j, dot))
    for j, body in enumerate(bodies):
        if body != 'moon':
            dot = ax2.scatter([], [], [], color=bodies_data[body]['color'], s=2*dot_s)
            dots_ax2.append((j, dot)) 
    
    hline, = ax3.plot([], [], color='grey', linewidth=3, alpha=0.8)
    
    def update(frame):
        i = indices[frame]
        time = times.utc_datetime()[i]
        
        changed_objects = []
        for j, dot in dots_ax1:
            dot._offsets3d = ([Loes.y[3*j][i]], [Loes.y[3*j+1][i]], [Loes.y[3*j+2][i]])
            changed_objects.append(dot)
        for j, dot in dots_ax2:
            dot._offsets3d = ([Loes.y[3*j][i]], [Loes.y[3*j+1][i]], [Loes.y[3*j+2][i]])
            changed_objects.append(dot)
              
        hline.set_data([time, time], [0, 7.5])
        changed_objects.append(hline)
        
        return tuple(changed_objects)

    ani = FuncAnimation(fig, update, frames=frames_N, interval=100, blit=True)
    plt.close(fig) 
    return ani

ani = create_animation_planets_b(Loes, Loes_l, N, set_bodies,150)
HTML(ani.to_html5_video())

Planetenbewegungen: Vergleich Python und C++ Simulation vs. Daten

Wir vergleichen nun die simulierten Daten mit den Ephemeriden-Daten der DE421-Datei. Wir führen dazu eine neue Simulation durch, welche die Anfangsdaten von Weihnachten 2000 (genauer: 24. Dezember 2000 um 20:30 Uhr) verwendet und bis nächste Weinachten (genauer: 24. Dezember 2026 um 20:30 Uhr) simuliert. Da die Erde nicht genau 365 Tage für eine Umrundung der Sonne benötigt, sondern etwa 365 Tage, 5 Stunden, 49 Minuten und 12 Sekunden, verwenden wir einen Wert von 365.2307715 Tagen für ein Jahr in unserer Simulation.

#  24. Dezember 2000 um 20:30 Uhr UTC
start_time = ts.utc(2000, 12, 24, 20, 30)
N = 10001
days = np.linspace(0, 26*365.2307715, N) 
times = start_time + days

Wir definieren das Python-Dictionary der kombinierten Massen-Gravitationsparameter der Himmelskörper neu und fügen diesem die Positions- und Geschwindigkeitswerte der gesamten neuen Zeitreihe zu.

# Gravitationsparameter GM in [AU^3 / d^2]
# Basierend auf JPL-Daten (DE421/DE440)
# Farben für Planeten-Plots, Sonnengelb, Merkur-Grau, Venus-Beige (Atmosphäre), Erd-Blau, Mond-Hellgrau, 
# Mars-Rostrot, Jupiter-Bänder-Braun, Saturn-Goldgelb, Uranus-Hellblau/Cyan, Neptun-Tiefblau, Pluto-Graubraun
bodies_data = {
    'sun': {'gm': 2.959122082855911e-04, 'color': '#FFCC00', 'label': r"$\rm \vec{r}_{sun}$"},
    'mercury': {'gm': 4.91248045036476e-11, 'color': '#A5A5A5', 'label': r"$\rm \vec{r}_{mercury}$"},
    'venus': {'gm': 7.243452486161142e-10, 'color': '#E3BB76', 'label': r"$\rm \vec{r}_{venus}$"},
    'earth': {'gm': 8.887692390113509e-10, 'color': '#2E5FAC', 'label': r"$\rm \vec{r}_{earth}$"},
    'moon': {'gm': 1.093189450742374e-11, 'color': '#D3D3D3', 'label': r"$\rm \vec{r}_{moon}$"},
    'mars': {'gm': 9.549535105170911e-11, 'color': '#E27B58', 'label': r"$\rm \vec{r}_{mars}$"},
    'jupiter barycenter': {'gm': 2.825345909524226e-07, 'color': '#D39C7E', 'label': r"$\rm \vec{r}_{jupiter}$"},
    'saturn barycenter': {'gm': 8.459706073245032e-08, 'color': '#C5AB6E', 'label': r"$\rm \vec{r}_{saturn}$"},
    'uranus barycenter': {'gm': 1.292024825783306e-08, 'color': '#B5E2E2', 'label': r"$\rm \vec{r}_{uranus}$"},
    'neptune barycenter': {'gm': 1.524357347885194e-08, 'color': '#4B70DD', 'label': r"$\rm \vec{r}_{neptune}$"},
    'pluto barycenter': {'gm': 2.17844105333832e-12, 'color': '#96745E', 'label': r"$\rm \vec{r}_{pluto}$"}
}

# Daten der Positionen und Geschwindigkeiten für alle Zeiten in "times" abrufen und in bodies_data speichern
for name in bodies_data.keys():
    bodies_data[name]['r'] = eph[name].at(times).position.au 
    bodies_data[name]['v'] = eph[name].at(times).velocity.au_per_d 

Wir führen nun zwei sehr genaue Python-Simulationen (rtol=atol=$10^{-13}$) durch. Eine simulierte nur die inneren Himmelskörper (Merkur, Venus, Erde, Mond und Mars) und die andere alle im Python-Dictionary definierten Himmelskörper.

set_m = []
r_0 = []
v_0 = []
for name in set_bodies:
    set_m.append(bodies_data[name]['gm'])
    r_0.append(bodies_data[name]['r'][:, 0])
    v_0.append(bodies_data[name]['v'][:, 0])
u_0 = np.concatenate([np.concatenate(r_0), np.concatenate(v_0)])

set_bodies_6 = ['sun','mercury','venus','earth','moon','mars']
set_n_6 = len(set_bodies_6)

set_m_6 = []
r_0 = []
v_0 = []
for name in set_bodies_6:
    set_m_6.append(bodies_data[name]['gm'])
    r_0.append(bodies_data[name]['r'][:, 0])
    v_0.append(bodies_data[name]['v'][:, 0])

u_0_6 = np.concatenate([np.concatenate(r_0), np.concatenate(v_0)])

fehler = 1e-13
t_span = [0,26*365.2307715]
tval = np.linspace(0, t_span[1], N)
Loes = solve_ivp(DGLsys_n, t_span, u_0, t_eval=tval, args=(set_m,set_n), rtol=fehler, atol=fehler)
Loes_6 = solve_ivp(DGLsys_n, t_span, u_0_6, t_eval=tval, args=(set_m_6,set_n_6), rtol=fehler, atol=fehler)

Die untere Abbildung vergleicht nun die Simulation, in der nur die Inneren Himmelskörper verwendet wurden mit den Daten der Ephemeriden-Datei 'de421'.

plt.figure(figsize=(11, 5))
plt.title(r"Abweichungen $\Delta x$ der Python-Lösung von den de421-Ephemeriden-Daten (nur innere Planeten)")
plt.xlabel("t")
plt.ylabel(r"$\Delta$x")
for i in range(len(set_bodies_6)):
    if i != set_bodies.index("moon"):
        plt.plot(times.utc_datetime(),Loes_6.y[3*i]-bodies_data[set_bodies_6[i]]['r'][0,:],color=bodies_data[set_bodies_6[i]]['color'], label=bodies_data[set_bodies_6[i]]['label'])
plt.legend(loc='upper left',fontsize=14);

Man erkennt, dass die simulierten inneren Planeten wohl einem systematischen Fehler unterliegen. Wie wir im Folgenden sehen werden ist der Fehler wirklich nicht von numerischer Natur, sondern die nicht simulierten anderen Planeten wirken auch gravitativ auf die inneren Himmelskörper und verursachen den in der oberen Abbildung gezeigten Fehler.

Die unetere Abbildung zeigt nun die Abweichungen der Simulation in welcher alle Planeten berücksichtigt wurden.

plt.figure(figsize=(11, 5))
plt.title(r"Abweichungen $\Delta x$ der Python-Lösung von den de421-Ephemeriden-Daten (alle Planeten)")
plt.xlabel("t")
plt.ylabel(r"$\Delta$x")
for i in range(len(set_bodies)):
    if i != set_bodies.index("moon"):
        plt.plot(times.utc_datetime(),Loes.y[3*i]-bodies_data[set_bodies[i]]['r'][0,:],color=bodies_data[set_bodies[i]]['color'], label=bodies_data[set_bodies[i]]['label'])
plt.legend(loc='upper left',fontsize=14);

Die Simulation stimmt nun gut mit den Daten der Ephemeriden-Datei überein (siehe obere Abbildung) und auch die simulierte Bewegung des Mondes (siehe untere Abbildung) zeigt nur kleine Abweichungen. Der hauptsächliche Grund, für die Abweichungen liegt hier wiederum nicht in der Numerik, sondern es handelt sich auch hier um einen kleinen systematischen Fehler, der dadurch bedingt ist, dass die Daten in der Ephemeriden-Datei Korrekturen verursacht durch die Allgemeine Relativitätstheorie mit in ihre Berechnungen einbeziehen. Die Allgemeine Relativitätstheorie im Vergleich zur Newtonschen Gravitation sorgt unter anderem für die Periheldrehung (besonders stark bei Merkur) und eine leicht veränderte gravitative Anziehung in Sonnennähe.

plt.figure(figsize=(11, 5))
plt.title(r"Abweichungen $\Delta x$ der Python-Lösung von den de421-Ephemeriden-Daten (nur Mond)")
plt.xlabel("t")
plt.ylabel(r"$\Delta$x")
for i in range(len(set_bodies)):
    if i == set_bodies.index("moon"):
        plt.plot(times.utc_datetime(),Loes.y[3*i]-bodies_data[set_bodies[i]]['r'][0,:],color=bodies_data[set_bodies[i]]['color'], label=bodies_data[set_bodies[i]]['label'])
plt.legend(loc='upper left',fontsize=14);

Wir stellen uns die Simulation wieder in einer Animation dar.

def create_animation_planets_c(Loes, N, bodies, bodies_data, frames_N=100, trail_length=50):
    dN = int(N / (frames_N + 1))
    dot_s = 15
    n = len(bodies)
    inner_bodies = ['sun', 'mercury', 'venus', 'earth', 'moon', 'mars']
    
    fig = plt.figure(figsize=(10, 6), facecolor='white')
    gs = gridspec.GridSpec(1, 2, width_ratios=[1.5, 1], wspace=0.2, top=0.75)
    
    ax1 = fig.add_subplot(gs[0], projection='3d')
    ax2 = fig.add_subplot(gs[1], projection='3d')
    
    dots1, trails1 = [], []
    dots2, trails2 = [], []
    all_legend_handles = [] 
    
    for i, name in enumerate(bodies):
        color = bodies_data[name]['color']
        label = bodies_data[name]['label']
        
        # Hauptplot: Alle Planeten
        ax1.plot(Loes.y[3*i], Loes.y[3*i+1], Loes.y[3*i+2], color=color, lw=0.4, alpha=0.5, linestyle=":")
        d1 = ax1.scatter([], [], [], color=color, s=dot_s, label=label)
        t1, = ax1.plot([], [], [], color=color, lw=1.2, alpha=0.3)
        dots1.append(d1)
        trails1.append(t1)
        all_legend_handles.append(d1)
        
        # Fokusplot: nur innere Planeten
        if name in inner_bodies:
            d2 = ax2.scatter([], [], [], color=color, s=dot_s)
            t2, = ax2.plot([], [], [], color=color, lw=1.2, alpha=0.6)
            dots2.append((i, d2, t2))

    ax2.legend(handles=all_legend_handles, 
               loc='lower center', 
               bbox_to_anchor=(0.5, 0.95), 
               ncol=4, 
               fontsize=10, 
               frameon=True,
               edgecolor='gray')

    # Achsen-Limits 
    limit_main = max(abs(Loes.y[0::3].min()), abs(Loes.y[0::3].max())) * 1.1
    limit_inner = 2.0 

    for ax, lim in zip([ax1, ax2], [limit_main, limit_inner]):
        ax.set_xlim(-lim, lim)
        ax.set_ylim(-lim, lim)
        ax.set_zlim(-lim, lim)

    def update(frame):
        idx = frame * dN
        progress = frame / frames_N
        
        # Änderung im Viewpoint
        ax1.view_init(elev=25 + (90 - 25) * progress, azim=10 + (0 - 10) * progress)
        ax2.view_init(elev=25 + (90 - 25) * progress, azim=10 + (0 - 10) * progress)
        
        changed_objects = []
        
        # Update Hauptplot
        for i in range(n):
            dots1[i]._offsets3d = ([Loes.y[3*i][idx]], [Loes.y[3*i+1][idx]], [Loes.y[3*i+2][idx]])
            s_idx = max(0, idx - trail_length * dN)
            trails1[i].set_data(Loes.y[3*i][s_idx:idx+1], Loes.y[3*i+1][s_idx:idx+1])
            trails1[i].set_3d_properties(Loes.y[3*i+2][s_idx:idx+1])
            changed_objects.extend([dots1[i], trails1[i]])

        # Update Fokusplot
        for i, d2, t2 in dots2:
            d2._offsets3d = ([Loes.y[3*i][idx]], [Loes.y[3*i+1][idx]], [Loes.y[3*i+2][idx]])
            s_idx = max(0, idx - trail_length * dN)
            t2.set_data(Loes.y[3*i][s_idx:idx+1], Loes.y[3*i+1][s_idx:idx+1])
            t2.set_3d_properties(Loes.y[3*i+2][s_idx:idx+1])
            changed_objects.extend([d2, t2])
            
        return tuple(changed_objects)

    ani = FuncAnimation(fig, update, frames=frames_N, interval=150, blit=True)
    plt.close(fig) 
    return ani

ani = create_animation_planets_c(Loes, N, set_bodies,bodies_data,150,50)
HTML(ani.to_html5_video())

Wir vergleichen nun die Simulation mit der C++-Version. Die Simulation benutzte die gleichen Massen und Anfangswerte der Planeten und simulierte auch für einen Zeitbereich von 26 Jahren. Der Zeitbereich wurde in $N_{C++}=500000$ zeitlichen Gitterpunkten unterteilt, wobei wir in der Ausgabedatei ("Planeten.dat") nur jeden 50-ten Punkt ausgaben.

data_C = np.genfromtxt("./Planeten.dat")  # Daten des C++ Simulation einlesen

plt.figure(figsize=(11, 5))
plt.title(r"Abweichungen $\Delta x$ der C++-Lösung von den de421-Ephemeriden-Daten")
plt.xlabel("t")
plt.ylabel(r"$\Delta$x")
for i in range(len(set_bodies)):
    if i != set_bodies.index("moon"):
        plt.plot(times.utc_datetime(),data_C[:,2+i*6]-bodies_data[set_bodies[i]]['r'][0,:],color=bodies_data[set_bodies[i]]['color'], label=bodies_data[set_bodies[i]]['label'])
plt.legend(loc='upper left',fontsize=10);

Die C++-Simulation stimmt auch gut mit den Daten der Ephemeriden-Datei überein (siehe obere Abbildung) und auch die simulierte Bewegung des Mondes (siehe untere Abbildung) zeigt nur kleine Abweichungen. Die Restabweichungen liegen wieder an dem systematischen Fehler der Nichtberücksichtigung der Effekte der Allgemeinen Relativitätstheorie.

plt.figure(figsize=(11, 5))
plt.title(r"Abweichungen $\Delta x$ der C++-Lösung von den de421-Ephemeriden-Daten (nur Mond)")
plt.xlabel("t")
plt.ylabel(r"$\Delta$x")
for i in range(len(set_bodies)):
    if i == set_bodies.index("moon"):
        plt.plot(times.utc_datetime(),data_C[:,2+i*6]-bodies_data[set_bodies[i]]['r'][0,:],color=bodies_data[set_bodies[i]]['color'], label=bodies_data[set_bodies[i]]['label'])
plt.legend(loc='upper left',fontsize=10);

Die untere Abbildung vergleicht nun die C++-Simulation mit der Python-Simulation. Man erkennt einen sehr kleinen Fehler von $\Delta x \approx 10^{-7}$.

plt.figure(figsize=(11, 5))
plt.title(r"Abweichungen $\Delta x$ der Python mit der C++-Lösung")
plt.xlabel("t")
plt.ylabel(r"$\Delta$x")
for i in range(len(set_bodies)):
    if i != set_bodies.index("moon"):
        plt.plot(times.utc_datetime(),data_C[:,2+i*6]-Loes.y[3*i],color=bodies_data[set_bodies[i]]['color'], label=set_bodies[i])
plt.legend(loc='lower left',fontsize=10);

Betrachtet man Merkur und Erde nicht, dann reduziert sich der Fehler sogar auf $\Delta x \approx 10^{-10}$.

plt.figure(figsize=(11, 5))
plt.title(r"Abweichungen $\Delta x$ der Python mit der C++-Lösung (ohne Merkur und Erde)")
plt.xlabel("t")
plt.ylabel(r"$\Delta$x")
for i in range(len(set_bodies)):
    if i not in [set_bodies.index("moon"),set_bodies.index("mercury"),set_bodies.index("earth")]:
        plt.plot(times.utc_datetime(),data_C[:,2+i*6]-Loes.y[3*i],color=bodies_data[set_bodies[i]]['color'], label=set_bodies[i])
plt.legend(loc='lower left',fontsize=10);

Lediglich die Mond-Berechnungen unterscheiden sich noch um $\Delta x \approx 10^{-6}$.

plt.figure(figsize=(11, 5))
plt.title(r"Abweichungen $\Delta x$ der Python mit der C++-Lösung (nur Mond)")
plt.xlabel("t")
plt.ylabel(r"$\Delta$x")
for i in range(len(set_bodies)):
    if i == set_bodies.index("moon"):
        plt.plot(times.utc_datetime(),data_C[:,2+i*6]-Loes.y[3*i],color=bodies_data[set_bodies[i]]['color'], label=set_bodies[i])
plt.legend(loc='lower left',fontsize=10);