Planetenbewegungen und das n-Körper-Problem
Die Bewegung eines Körpers in einem gravitativen Zentralkraftfeld ist ein oft behandeltes System in der klassischen Mechanik (siehe z. B. Walter Greiner, Mechanik Teil 1 [5. Auflage, 1989, Kapitel 26, S. 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, S. 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ühmten Bestseller Die drei Sonnen des chinesischen Science-Fiction-Autors 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 zugrunde liegende 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$ die Gravitationskonstante 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 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.
Wir werden im Laufe dieses Unterkapitels auch Bewegungen mit mehr als drei Körpern simulieren und stellen deshalb die obere Bewegungsgleichung auch schon an dieser Stelle für $n$-Körper dar. Wir behandeln dann ein sogenanntes n-Körperproblem, wobei $n$ die Zahl der simulierten Körper darstellt. Die zugrunde liegende 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 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.
Numerische Lösung des Dreikörperproblems mit Python
Wir betrachten im Folgenden die numerische Lösung des Dreikörperproblems mittels Python. 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. Beim Klicken auf das nebenstehende Bild gelangen Sie zu dem Jupyter Notebook 3-Body-Problem.ipynb in dem die numerischen Lösungen des Dreikörperproblems und n-Körperproblems in Python programmiert sind.
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 $), nehmen eine planare Bewegung an ($z(t)=0 \,\, \forall \, t$) und setzen die Gravitationskonstante auf Eins ($G=1$).
Zum Lösen des gekoppelten Systems von Differenzialgleichungen benutzen wir die im Python-Modul SciPy definierte Funktion 'solve_ivp(...)', die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differenzialgleichungen bereitstellt. Wir betrachten zunächst eine Anfangsbedingung mit einer Exzentrizität von $e=0.7$ und lösen das Anfangswertproblem im Bereich $t \in [0,4]$.
Wir stellen nun die Bewegung der drei Körper als eine Animation im Konfigurationsraum $(x,y)$ dar. Neben der Visualisierung im Konfigurationsraum ist es oft auch sinnvoll die zeitliche Entwicklung im Phasenraum der radialen Zentrumskoordinaten darzustellen. 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$, wobei der Index den jeweiligen Körper kennzeichnet). Zusätzlich berechnen wir auch die radialen Geschwindigkeiten der Körper relativ zum gemeinsamen Zentrum ($v_1, v_2$ und $v_3$). Details der Berechnungen sind in dem Jupyter Notebook 3-Body-Problem.ipynb dargestellt.
Die Animation oben links zeigt die Bewegung der Körper im Konfigurationsraum $(x,y)$ und stellt zusätzlich die Trajektorien im Phasenraum der radialen Zentrumskoordinaten dar ($t \in [0,4]$). Anhand der grauen Zentrums-Verbindungslinien erkennt man die periodische Veränderung der radialen Entfernung der Körper. Die Körper selbst bilden während der Bewegung zu jedem Zeitpunkt die Form eines gleichseitigen Dreiecks. Der Schwerpunkt der Körper ist im Ursprung und verändert sich mit der Zeit nicht (siehe Diagramm rechts unten, beachte x-y-Skala).
Die dargestellte elliptische Lagrangesche Lösung ist jedoch instabil gegen kleine Änderungen. In der Animation oben rechts wurde hierfür das zu simulierende Zeitfenster auf $t \in [0,12.5]$ erweitert und man sieht, dass ab einem gewissen Zeitpunkt chaotisches Verhalten entsteht. 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.
Die 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 (die 'Figur Acht', siehe Animation unten links) besitzt hingegen numerische Stabilität und stellt eine wunderschöne, hochsymmetrische und gleichzeitig robuste Insel der Ordnung inmitten der sonst oft chaotischen Bewegungen des Dreikörperproblems dar.
Bei dieser Lösung folgen alle drei Körper exakt der gleichen 'Figur Acht'-Flugbahn, jedoch zeitlich versetzt (eine so genannte 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 schwingen die drei Körper so umeinander, dass das System als Ganzes nicht rotiert (gesamter Drehimpuls ist Null, siehe Abbildung oben rechts).
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.
Im Gegensatz zu den bisher dargestellten Lösungen betrachten die Wissenschaftler Systeme mit ungleichen Massenwerten, bei denen Stabilität im Allgemeinen einfacher zu finden ist (siehe Periodic orbits for Newtonian planar three-body problem with unequal mass). Die Animation oben links zeigt eine Simulation mit den folgenden Massenwerten $m_1 = 2$, $m_2 = 0.9$ und $m_3 = 1$. Der Verlauf der Trajektorien, sowohl im Konfigurationsraum, als auch in den radialen Phasenräumen der Körper, sieht bei dieser Bewegung komplizierter aus. Auch die zeitliche Entwicklung der einzelnen Drehimpulskomponenten ist vielschichtiger als bei der 'Figur Acht' (siehe Abbildung oben rechts).
Alle bisher dargestellten Bewegungen waren aufgrund der speziellen Wahl der Anfangswerte von periodischer, regulärer Struktur, wobei einige (z. B. die dargestellte elliptische Lagrange-Lösung) nach einiger Zeit chaotisches Verhalten zeigten. Bewegungen von Dreikörper-Systemen mit gleichen Massenwerten sind in der Regel chaotisch und die Animation links unten zeigt ein typisches Beispiel einer chaotischen Bewegung ($t \in [0,45]$). 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, kehrt dann wieder zu den anderen beiden zurück und separiert sich dann jedoch erneut. Die Bewegung eines solchen Systems kann man nicht adäquat vorhersagen, da sich kleine, $\epsilon$-ähnliche Störungen nach einiger Zeit exponentiell vergrößern.
Die bisher betrachteten Bewegungen der drei Körper fanden alle auf der ($z=0$)-Ebene statt (planare Bewegungen). 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 (siehe Li, Xiaoming, and Shijun Liao. "Discovery of 10,059 new three-dimensional periodic orbits of general three-body problem."), von denen eine in der oberen rechten Animation dargestellt ist. Die Bewegung der drei Körper im dreidimensionalen ($x,y,z$)-Konfigurationsraum und im Phasenraum der Zentrumskoordinaten zeigt, dass sich die Trajektorien des roten und blauen Körpers überlagern. Die Autoren bezeichnen diese Bewegungsart als sogenannte "piano-trio orbits" (ein Trio für zwei Violinen und ein Klavier).
Numerische Lösung des Dreikörperproblems mit C++
Im Unterpunkt Abgeleitete Klassen, Vererbung von Klassenmerkmalen und Klassenhierarchien der Vorlesung hatten wir mittels der abstrakten Schnittstellen-Basisklasse 'dsolve' ein System von Differenzialgleichungen numerisch mittels des Runge-Kutta Ordnung vier Verfahrens gelöst. In dem folgenden C++ Header File N-BodyProblem.hpp ist die Subklasse 'nBody' von der abstrakten Klasse 'dsolve' abgeleitet. Mittels der Subklasse 'nBody' lassen sich Simulationen des allgemeinen n-Körper-Problems mit unterschiedlichen Massenwerten und Anfangsbedingungen erstellen. Die Subklasse überschreibt dabei die rein virtuelle Funktion 'dgls' und definiert somit die Bewegungsgleichung des n-Körper-Problems.
Mittels der zwei definierten Konstruktoren der 'nBody'-Klasse ist es nun möglich, eine Instanz eines 'nBody'-Objektes zu erstellen. Die erste Variante, Konstruktor 1, besitzt sechs Argumente ( nBody(Grav.Konstante G, Massenvektor der Körper, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte u_init) ) und mit ihm kann der Benutzer die Bewegung einer allgemeinen Anfangsbedingung des n-Körper-Problems simulieren. Die zweite Variante, Konstruktor 2, besitzt nur vier Argumente ( 'nBody(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Exzentrizität exz) ) und mit ihm kann der Benutzer die Bewegung einer elliptische Lagrangeschen homografischen Lösung mit Exzentrizität exz:=$e \in [0,1]$ des 3-Körper Problems simulieren. Zusätzlich werden in der 'nBody'-Klasse noch zwei öffentliche Methoden definiert, die den Ort und die Geschwindigkeit der Körper an einer gewünschten zeitlichen Stützstelle in kartesischen Koordinaten ausgeben.
/* Header-Datei * Das n-Körperproblem in drei Dimensionen * n gravitativ wirkende Körper bewegen sich durch ihre gegenseitige Newtonsche Anziehungskraft * Berechnung der Lösung eines Systems von 2*3*n Differentialgleichungen erster Ordnung * mittels Runge-Kutta Ordnung vier Methode * Verfahren zur Loesung der DGL ist in einer abstrakten Basisklasse ausgelagert * Berechnung der Lösung innerhalb einer Methode (Klassenfunktion) Runge-Kutta-Verfahren * Überschreibung der virtuellen Funktion der Bewegungsgleichungen in der abgeleiteten Klasse (Subklasse) 'nBody' * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b] * Unterklasse 'nBody' * Konstruktor 1: 'nBody'(Grav.Konstante G, Massenvektor der Körper, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte u_init) * Konstruktor 2: Elliptische Lagrangesche homografischen Lösungen (Exzentrizität exz in [0,1]) * 'nBody'(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Exzentrizität exz) * Ausgabe zum Plotten in Datei */ #include <cmath> // Bibliothek für mathematisches (e-Funktion, Betrag, ...) #include <vector> // Vector-Container der Standardbibliothek using namespace std; // Benutze den Namensraum std // Abstrakte Basisklasse 'dsolve' class dsolve { public: // Daten-Member der Klasse double a; // Untergrenze des Zeit-Intervalls double b; // Obergrenze des Intervalls int N; // Anzahl der Punkte vector<double> alpha; // Anfangswerte vector<vector<double>> ym; // Lösungsmatrix vector<double> t; // Zeit-Vektor // Virtuelle Funktion dgls, wird an anderer Stelle definiert virtual vector<double> dgls(double t, vector<double> u_vec) = 0; // Konstruktor mit vier Argumenten (Initialisierung der Daten-Member) dsolve(double a_, double b_, int N_, vector<double> alpha_) :a(a_), b(b_), N(N_), alpha(alpha_) { t.push_back(a_); // Zum Zeit-Vektor die Anfangszeit eintragen ym.push_back(alpha_); // Zur Loesungs-Matrix die Anfangswerte eintragen } // Methode der Runge-Kutta-Methode void solve() { double h = (b - a)/N; // Abstand dt zwischen den aequidistanten Punkten des t-Intervalls (h=dt) int n = alpha.size(); // Anzahl der Anfangswerte bzw. DGL-Gleichungen (hier n=2) vector<double> k1(n),k2(n),k3(n),k4(n); // Deklaration der vier Runge-Kutta Parameter fuer jede Loesung vector<double> y(n); // Temporaerer Hilfsvektor for (int i = 0; i < N; ++i) { // for-Schleife ueber die einzelnen Punkte des t-Intervalls for (int j = 0; j < n; ++j) k1[j] = h * dgls(t[i], ym[i])[j]; // for-Schleife ueber die n Runge-Kutta k1-Parameter for (int j = 0; j < n; ++j) y[j] = ym[i][j] + k1[j] / 2; // Temporaerer Hilfsvektor for (int j = 0; j < n; ++j) k2[j] = h * dgls(t[i] + h / 2, y)[j]; // for-Schleife ueber die n Runge-Kutta k2-Parameter for (int j = 0; j < n; ++j) y[j] = ym[i][j] + k2[j] / 2; // Temporaerer Hilfsvektor for (int j = 0; j < n; ++j) k3[j] = h * dgls(t[i] + h / 2, y)[j]; // for-Schleife ueber die n Runge-Kutta k3-Parameter for (int j = 0; j < n; ++j) y[j] = ym[i][j] + k3[j]; // Temporaerer Hilfsvektor for (int j = 0; j < n; ++j) k4[j] = h * dgls(t[i] + h, y)[j]; // for-Schleife ueber die n Runge-Kutta k4-Parameter for (int j = 0; j < n; ++j) y[j] = ym[i][j] + (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6; // Temporaerer Hilfsvektor ym.push_back(y); // Zum Loesungs-Vektor den neuen Wert eintragen t.push_back(a + (i + 1) * h); // Zum Zeit-Vektor die neue Zeit eintragen } // Ende for-Schleife ueber die einzelnen Punkte des t-Intervalls } // Ende des Methode 'solve()' }; // Ende der Klasse 'dsolve' // Sub-Klasse 'nBody' class nBody : public dsolve { public: double G; // Gravitationskonstante G vector<double> m; // Vektor der Massen der Körper unsigned n; // Anzahl der Körper // Überschreibung der virtuellen Funktion dgls // Definition der Bewegungsgleichung des n-Körper-Problems in drei Dimensionen vector<double> dgls(double t, vector<double> u_vec){ unsigned dim = 3 * n; // Die n-Körper bilden 3*n Dimensionen vector<double> du_dt(2 * dim); // Das System von 2*3*n DGLs erster Ordnung werden als Standard-Vektor definiert for(unsigned i=0; i < dim; ++i){ // 1. Block von 3*n-DGLs erster Ordnung du_dt[i] = u_vec[dim + i]; // drdt = v } // Ende 1. Block for(unsigned i=0; i < n; ++i){ // 2. Block von 3*n-DGLs erster Ordnung double res_x = 0.0, res_y = 0.0, res_z = 0.0; // x,y,z-Werte der DGL des i-ten Körpers double xi = u_vec[3 * i]; // x-Posiion von Körper i double yi = u_vec[3 * i + 1]; // y-Posiion von Körper i double zi = u_vec[3 * i + 2]; // z-Posiion von Körper i for(unsigned j=0; j < n; ++j){ // Schleife j über alle Körper if (i != j) { // j-ter Körper wirkt auf i-ten Körper double xj = u_vec[3 * j]; // x-Posiion von Körper j double yj = u_vec[3 * j + 1]; // y-Posiion von Körper j double zj = u_vec[3 * j + 2]; // z-Posiion von Körper j // alle j-Körper wirken gravitativ auf den i-ten Körper (in drei Dimensionen, xyz) double r2 = (xj-xi)*(xj-xi) + (yj-yi)*(yj-yi) + (zj-zi)*(zj-zi); double f = G * m[j] / (r2 * sqrt(r2)); res_x += f * (xj-xi); res_y += f * (yj-yi); res_z += f * (zj-zi); } // Ende if } // Ende Schleife j über alle Körper du_dt[dim + 3 * i] = res_x; // DGL des i-ten Körpers in x-Richtung du_dt[dim + 3 * i + 1] = res_y; // DGL des i-ten Körpers in y-Richtung du_dt[dim + 3 * i + 2] = res_z; // DGL des i-ten Körpers in z-Richtung } // Ende 2. Block return du_dt; // Rueckgabewert der DGLs als Standard-Vektor } // Ende: Definition der Bewegungsgleichung-Funktion // 1. Konstruktor mit sechs Argumenten: 'nBody'(Grav.Konstante G, Massenvektor der Körper, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte u_init) nBody(double G_ = 1.0, vector<double> m_ = {1,1,1}, double a_ = 0, double b_ = 45.0, int N_ = 10000, vector<double> alpha_ = {1.0,0.0,0.0,-0.5,0.8,0.0,-0.5,-0.8,0.0,0.0,0.5,0.0,0.2,-0.1,0.0,-0.2,-0.4,0.0}) : dsolve(a_, b_, N_, alpha_), G(G_), m(m_) {n = static_cast<unsigned>(m_.size());} // 2. Konstruktor mit vier Argumenten: 'nBody'(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Exzentrizität exz) nBody(double a_ = 0, double b_ = 45.0, int N_ = 10000, double exz_ = 1.0) : dsolve(a_, b_, N_, alpha_Lagr(exz_)), G(1.0), m({1,1,1}), n(3) {} // Methode der Ausgabe der Orte der Körper am Stützpunkt i vector<double> r(size_t i) { vector<double> w(3 * n, 0.0); if (i < t.size()) { for(unsigned j=0; j < 3*n; ++j){ w[j] = ym[i][j]; } } return w; } // Methode der Ausgabe der Geschwindigkeiten der Körper am Stützpunkt i vector<double> v(size_t i) { vector<double> w(3 * n, 0.0); if (i < t.size()) { for(unsigned j=0; j < 3*n; ++j){ w[j] = ym[i][j + 3*n]; } } return w; } private: // Elliptische Lagrangesche homografischen Lösungen (Exzentrizität exz) static vector<double> alpha_Lagr(double exz) { double c1 = 1.0 / sqrt(3.0); double c2 = sqrt(1.0 - exz); double c3 = sqrt(3.0*(1.0 - exz)) / 2.0; return {c1,0.0,0.0,-c1/2.0,0.5,0.0,-c1/2.0,-0.5,0.0,0.0,c2,0.0,-c3,-c2/2.0,0.0,c3,-c2/2.0,0.0}; } }; // Ende der Unterklasse 'nBody'
Am Ende des Header File wird schließlich die private Funktion 'alpha_Lagr' definiert, die für den zweiten Konstruktor die elliptische Lagrangesche Anfangsbedingungen erzeugt.
Im Hauptprogramm N-BodyProblem.cpp (siehe unteres C++ Programm) wird dann mit dem zweiten Konstruktor eine Instanz der 'nBody'-Klasse erstellt (nBody L = nBody(a,b,N,0.7);) und die vererbte Funktion 'solve()' aufgerufen. Anstatt der verwendeten Lagrangeschen Anfangsbedingungen hätte man auch mittels des ersten Konstruktors eine selbst erstellte Anfangsbedingung verwenden können (siehe auskommentierten Zeilen für die 'Figur 8'-Anfangsbedingungen).
/* Das n-Körperproblem in drei Dimensionen * n gravitativ wirkende Körper bewegen sich durch ihre gegenseitige Newtonsche Anziehungskraft * Berechnung der Lösung eines Systems von 2*3*n Differentialgleichungen erster Ordnung * mittels Runge-Kutta Ordnung vier Methode * Verfahren zur Loesung der DGL ist in einer abstrakten Basisklasse ausgelagert * Berechnung der Lösung innerhalb einer Methode (Klassenfunktion) Runge-Kutta-Verfahren * Überschreibung der virtuellen Funktion der Bewegungsgleichungen in der abgeleiteten Klasse (Subklasse) 'nBody' * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b] * Unterklasse 'nBody' * Konstruktor 1: 'nBody'(Grav.Konstante G, Massenvektor der Körper, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte u_init) * Konstruktor 2: Elliptische Lagrangesche homografischen Lösungen (Exzentrizität exz in [0,1]) * 'nBody'(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Exzentrizität exz) * Ausgabe zum Plotten in Datei */ #include "N-BodyProblem.hpp" // Header-Datei des n-Körperproblem in drei Dimensionen #include <iostream> // Standard Input- und Output Bibliothek in C, z.B. fprintf(...) #include <vector> // Vector-Container der Standardbibliothek using namespace std; // Benutze den Namensraum std // Hauptfunktion int main(){ double a = 0; // Untergrenze des Zeit-Intervalls double b = 12.4; // Obergrenze des Intervalls int N = 100000; // Anzahl der Punkte FILE *ausgabe; // Deklaration eines Files fuer die Ausgabedatei ausgabe = fopen("N-BodyProblem.dat", "w+"); // Ergebnisse werden in die Datei 'N-BodyProblem.dat' geschrieben // Anfangswerte für die 'Figur Acht' (Orte und Geschwindigkeiten der drei Körper in drei Dimensionen) // vector<double> u_init {0.97000436,-0.24308753,0.0,-0.97000436,0.24308753,0.0,0.0,0.0,0.0, // 0.46620368,0.43236573,0.0,0.46620368,0.43236573,0.0,-0.93240737,-0.86473146,0.0}; // nBody L = nBody(1,{1,1,1},a,b,N,u_init); // Konstruktor 1 bilded eine Instanz der Klasse nBody // Anfangswerte für die Elliptische Lagrangesche homografischen Lösungen (Exzentrizität exz=0.7) nBody L = nBody(a,b,N,0.7); // Konstruktor 2 bilded eine Instanz der Klasse nBody L.solve(); // Methode des Runge-Kutta-Verfahrens ausführen fprintf(ausgabe, "# 0: Index i \n# 1: t-Wert \n# 2:x_1 \n# 3:vx_1 \n# 4:y_1 \n# 5:vy_1 \n# 6:z_1 ... \n"); // Beschreibung der ausgegebenen Groessen for(size_t i=0; i < L.t.size(); i += 10){ // for-Schleife über die Zeitgitterpunkte (nur jeden 10-ten Wert ausgeben) fprintf(ausgabe, "%3ld %19.15f ",i,L.t[i]); // Ausgabe Index und Zeit for(unsigned j=0; j < 2*3*L.n; ++j){ // for-Schleife zur Ausgabe der Orte und Geschw. der n-Körper fprintf(ausgabe, "%19.15f %19.15f ",L.r(i)[j],L.v(i)[j]); // Ausgabe xyz-Werte und Geschwindigkeiten der Körper } // Ende der for-Schleife; n-Körper fprintf(ausgabe, "\n"); // Neue Zeile } // Ende der for-Schleife; Zeitgitterpunkte } // Ende der Hauptfunktion
Am Ende werden die berechneten Lösungswerte der Simulation in die Datei 'N-BodyProblem.dat' ausgegeben. Mittels des Python-Skriptes 3-BodyProblem.py kann man dann eine Animation der C++-Simulation erstellen.
Planetenbewegungen und das n-Körper-Problem
Wir möchten nun eine Simulation mit mehr als drei Körpern machen und als Beispiel die Planetenbewegungen unseres Sonnensystems betrachten. Die Bewegungsgleichungen des n-Körperproblems können mittels des C++ Header File N-BodyProblem.hpp simuliert werden, wobei man die unterschiedlichen Massenwerte und Anfangsbedingungen der Himmelskörper z. B. mittels der Python-Bibliothek Skyfield generiert. 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 verwenden im Folgenden die kombinierten Massen-Gravitationsparameter $\mu_j = G \cdot m_j$ und lesen die Positionen und Geschwindigkeiten der Himmelskörper an einem festen Datum (hier speziell 24. Dezember 2000 um 20:30 Uhr UTC) ein. Wir extrahierten die Daten in den Einheiten [AU] (Astronomische Einheit, 1 AU = 149597870700 m) und Tagen (d) und exportierten das gesamte Python-Dictionary in eine C++-map (siehe untere C++ Datei Planeten.cpp).
/* Planetenbewegungen in drei Dimensionen * basierend auf dem n-Körperproblem in drei Dimensionen * n gravitativ wirkende Körper bewegen sich durch ihre gegenseitige Newtonsche Anziehungskraft * Berechnung der Lösung eines Systems von 2*3*n Differentialgleichungen erster Ordnung * mittels Runge-Kutta Ordnung vier Methode * Verfahren zur Loesung der DGL ist in einer abstrakten Basisklasse ausgelagert * Berechnung der Lösung innerhalb einer Methode (Klassenfunktion) Runge-Kutta-Verfahren * Überschreibung der virtuellen Funktion der Bewegungsgleichungen in der abgeleiteten Klasse (Subklasse) 'nBody' * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b] * Unterklasse 'nBody', Konstruktor 1 verwendet * Konstruktor 1: 'nBody'(Grav.Konstante G, Massenvektor der Körper, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte u_init) * Konstruktor 2: Elliptische Lagrangesche homografischen Lösungen (Exzentrizität exz in [0,1]) * 'nBody'(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Exzentrizität exz) * Ausgabe zum Plotten in Datei "Planeten.dat" */ #include "N-BodyProblem.hpp" // Header-Datei des n-Körperproblem in drei Dimensionen #include <iostream> // Standard Input- und Output Bibliothek in C, z.B. fprintf(...) #include <vector> // Vector-Container der Standardbibliothek #include <string> // Standard-Header-Datei für Zeichenketten #include <map> // Standard-Header-Datei des assoziativer Containers 'map' für Key-Value Zuordnungen using namespace std; // Benutze den Namensraum std // Struktur zum georneten Speichern der Himmelskörper-Parameter struct Body { double gm; // gm=G*m: Gravitationskonstante * Masse des Körpers vector<double> r; // Anfangsort in drei Dimensionen vector<double> v; // Anfangsogeschwindigkeit in drei Dimensionen }; // Hauptfunktion int main(){ // Definition der gm-Parameter und Anfangswerte der Himmelskörper unseres Sonnensystems in einer C++-map // Ephemeriden-Datei 'de421.bsp': 24. Dezember 2000 um 20:30 Uhr UTC) // Länge in Astronomischen Einheitn (AU = 149.597.870,7 km), Zeit in Tagen [d] map<string, Body> bodies_data = { {"sun", {0.00029591220828559, {-0.00469206964217851, -0.00455284299129719, -0.00180351606161285}, {0.00000808182754769, -0.00000364991309399, -0.00000177739810148}}}, {"mercury", {0.00000000004912480, {0.00869081881722255, -0.41212824812482196, -0.22090270075603882}, {0.02248815298088397, 0.00281999190824476, -0.00082493351551472}}}, {"venus", {0.00000000072434525, {0.58492431879299430, 0.39106068887197776, 0.13886288700490920}, {-0.01176860984756512, 0.01468495634135895, 0.00735174095931183}}}, {"earth", {0.00000000088876924, {-0.06209037112534756, 0.89629552530222610, 0.38875814156074100}, {-0.01745436390436983, -0.00098293944800350, -0.00042570411094977}}}, {"moon", {0.00000000001093189, {-0.06238174828301801, 0.89381342029917210, 0.38778197364226596}, {-0.01689334413161172, -0.00103583077929171, -0.00050101150118707}}}, {"mars", {0.00000000009549535, {-1.65650069587999330, 0.01013256724249670, 0.04958858122190123}, {0.00024579442831315, -0.01163534295071834, -0.00534326234956417}}}, {"jupiter barycenter", {0.00000028253459095, {1.84307596674717920, 4.32528771663242900, 1.80908857439581870}, {-0.00711591051273230, 0.00280462876522796, 0.00137544423420692}}}, {"saturn barycenter", {0.00000008459706073, {4.71735462953881500, 7.27073386854747500, 2.80001882890842420}, {-0.00506407973083511, 0.00258153583418907, 0.00128420716954018}}}, {"uranus barycenter", {0.00000001292024826, {15.35097038912193100, -11.59975177639601200, -5.29752281748167600}, {0.00248407812132207, 0.00261666373062273, 0.00111089476005838}}}, {"neptune barycenter", {0.00000001524357348, {17.71921742890633700, -22.37141413515812300, -9.59788796335357700}, {0.00251753419769610, 0.00174985139276685, 0.00065354736058049}}}, {"pluto barycenter", {0.00000000000217844, {-8.78866167066838600, -28.36872537395957700, -6.20501655731820000}, {0.00306982805604962, -0.00102661679740455, -0.00124530173577011}}} }; double G = 1.0; // Gravitationskonstante G=1, da Massenwerte in G*m double a = 0; // Untergrenze des Zeit-Intervalls double b = 365.2307715; // Obergrenze des Intervalls in Tagen (1 Jahr aprox. 365.2307715 Tage) int N = 10000; // Anzahl der Punkte der Stützstellen 1000000 FILE *ausgabe; // Deklaration eines Files fuer die Ausgabedatei ausgabe = fopen("Planeten.dat", "w+"); // Ergebnisse werden in die Datei 'Planeten.dat' geschrieben // String-Vektor der zu simulierenden Himmelskörper // const vector<string> set_bodies {"sun","mercury","venus","earth","moon","mars"}; const vector<string> set_bodies {"sun","mercury", "venus","earth","moon","mars","jupiter barycenter", "saturn barycenter", "uranus barycenter", "neptune barycenter", "pluto barycenter"}; const size_t n = set_bodies.size(); // Anzahl der zu simulierenden Himmelskörper vector<double> m(n); // Vektor der Massen der Körper (gm-Werte, gm=G*m) vector<double> u_init(6*n); // Vektor der Anfangswerte für die Orte und Geschwindigkeiten der n Körper in drei Dimensionen for (size_t i = 0; i < n; ++i) { // for-Schleife über alle Körper const auto& data = bodies_data.at(set_bodies[i]); // Daten des i-ten Körpers als Referenz auf die C++-map m[i] = data.gm; // gm-Wert in den m-Vektor schreiben copy(data.r.begin(), data.r.end(), u_init.begin() + (i * 3)); // Orte an den Anfang von u_init schreiben (Position i * 3) copy(data.v.begin(), data.v.end(), u_init.begin() + (n + i) * 3); // Geschwindigkeiten in die zweite Hälfte von u_init schreiben (Position: (n + i) * 3) } // Ende for-Schleife über alle Körper nBody L = nBody(G,m,a,b,N,u_init); // Konstruktor 1 bilded eine Instanz der Klasse nBody L.solve(); // Methode des Runge-Kutta-Verfahrens ausführen fprintf(ausgabe, "# 0: Index i \n# 1: t-Wert \n# 2:x_1 \n# 3:vx_1 \n# 4:y_1 \n# 5:vy_1 \n# 6:z_1 ... \n"); // Beschreibung der ausgegebenen Groessen for(size_t i=0; i < L.t.size(); i += 1){ // for-Schleife über die Zeitgitterpunkte (i += 10 -> nur jeden 10-ten Wert ausgeben) fprintf(ausgabe, "%3ld %19.15f ",i,L.t[i]); // Ausgabe Index und Zeit for(unsigned j=0; j < 2*3*L.n; ++j){ // for-Schleife zur Ausgabe der Orte und Geschw. der n-Körper fprintf(ausgabe, "%19.15f %19.15f ",L.r(i)[j],L.v(i)[j]); // Ausgabe xyz-Werte und Geschwindigkeiten der Körper } // Ende der for-Schleife; n-Körper fprintf(ausgabe, "\n"); // Neue Zeile } // Ende der for-Schleife; Zeitgitterpunkte } // Ende der Hauptfunktion
Mittels des Vektors 'const vector<string> set_bodies {...}' wird nun die Liste der zu simulierenden Himmelskörper definiert, dann werden die Massenwerte und Anfangsbedingungen dieser von der C++-map kopiert und in den Konstruktor als Argumente eingefügt. Der Konstruktor erzeugt dann eine Instanz der 'nBody'-Klasse (nBody L = nBody(G,m,a,b,N,u_init);) und dann wird die vererbte Funktion 'solve()' aufgerufen. Am Ende werden die berechneten Lösungswerte der Simulation in die Datei 'Planeten.dat' ausgegeben.
Am Ende des Jupyter Notebooks 3-Body-Problem.ipynb simulieren wir die Planetenbewegung und das n-Körperproblem auch mittels Python. Zunächst führen wir dort eine Simulation durch, bei der wir die Bewegung der inneren Planeten für ca. ein Jahr simulieren. Zusätzlich simulieren wir auch die Bewegung des Mondes um die kommenden Sonnenfinsternisse vorherzusagen. 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 = \arccos \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)$). Die Animation unten links zeigt die Ergebnisse der Simulation.
Die Animation stellt einerseits die Bewegung der inneren Planeten um die Sonne dar (linke Abbildung, oben links) und andererseits zeigt sie auch die Bewegung des Mondes um die Erde (linke Abbildung, oben rechts), wobei die gelbe Linie den Verbindungsvektor Erde-Sonne und die blaue Linie den Verbindungsvektor Erde-Mond skizziert. Die zeitliche Entwicklung des Winkels $\alpha(t)$ wird in der linken Abbildung im unteren Bereich gezeigt. Die simulierten Berechnungen des Winkels $\alpha$ stimmen dabei sehr gut mit den Daten der zukünftigen Sonnenfinsternisse überein (12. August 2026 in Nordspanien (u. a. Mallorca und Ibiza) und 6. Februar 2027 in Argentinien und Brasilien).
Die Animation auf der rechten Seite zeigt eine Simulation aller Planeten über 26 Jahre, beginnend vom 24. Dezember 2000. Diese Simulation wurde sowohl mit Python, als auch mit dem C++-Programm mit möglichst hoher Genauigkeit durchgeführt (Details siehe im Jupyter Notebooks 3-Body-Problem.ipynb). Die Abbildung unten links vergleicht beide Simulationen miteinander. Man erkennt einen sehr kleinen Fehler von $\Delta x \approx 10^{-7}$. Betrachtet man Merkur und Erde nicht, dann reduziert sich der Fehler sogar auf $\Delta x \approx 10^{-10}$.
Vergleicht man die Simulationen mit den Daten der Ephemeriden-Datei, so stimmen, sowohl die Python als auch die C++ Simulation, 'nur' auf einen Fehler von $\Delta x \approx 10^{-5}$ überein (siehe Abbildung unten rechts). Der hauptsächliche Grund, für die Abweichungen liegt hierbei nicht in der Numerik, sondern es handelt sich 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.
