In diesem Unterpunkt werden wir zeigen, wie man die Geodätengleichung numerisch mittels eines C++ Programms löst. Die Geodätengleichung $\frac{d^2 x^\mu}{d\lambda^2} + \Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} ~=~ 0$ ist ein System von vier gekoppelten nichtlinearen Differenzialgleichungen zweiter Ordnung \[ \begin{eqnarray} && \frac{d^2 x^0}{d\lambda^2} = \frac{d^2 t}{d\lambda^2} = - \Gamma^0_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} \\ && \frac{d^2 x^1}{d\lambda^2} = \frac{d^2 r}{d\lambda^2} = - \Gamma^1_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda}\\ && \frac{d^2 x^2}{d\lambda^2} = \frac{d^2 \theta}{d\lambda^2} = - \Gamma^2_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda}\\ && \frac{d^2 x^3}{d\lambda^2} = \frac{d^2 \phi}{d\lambda^2} = - \Gamma^3_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} \end{eqnarray} \] , wobei wir wieder ein sphärisches Koordinatensystem benutzen ($x^\mu=\left(t,r,\theta,\phi \right)$) und im Folgenden den affinen Parameter $\lambda$ gleich der Eigenzeit $\tau$ des Probekörpers setzen. In dem Jupyter Notebook Klassifizierung unterschiedlicher Bahnbewegungen ('GeodSchwarzschildII_OGRePy.ipynb' , View Notebook, Download Notebook) betrachten wir die Geodätengleichung in vorgegebener Schwarzschild-Raumzeit und beschränkten uns dabei auf ebene Bahnbewegungen des Testteilchens ($\theta=\frac{\pi}{2}, {\rm cos}(\theta)=0, {\rm sin}(\theta)=1, \frac{d \theta}{d\tau}=0, \frac{d^2 \theta}{d\tau^2}=0$). Unter diesen Annahmen sind nur drei der vier Geodätengleichungen relevant und diese lauten: \[ \begin{eqnarray} \frac{d^2 t}{d\tau^2} &=& \frac{2 M }{\left(2 M - r\right) r} \cdot \frac{d r}{d \tau} \cdot \frac{d t}{d \tau} \\ \frac{d^2 r}{d\tau^2} &=& \frac{M \left(2 M - r\right)}{r^{3}} \cdot \left(\frac{d t}{d \tau} \right)^{2} - \frac{M}{\left(2 M - r\right) r} \cdot \left(\frac{d r}{d \tau} \right)^{2} + \left(- 2 M + r\right) \cdot \left(\frac{d \phi}{d \tau}\right)^{2}\\ \frac{d^2 \phi}{d\tau^2} &=& - \frac{2}{r} \cdot \frac{d \phi}{d \tau} \cdot \frac{d r}{d \tau} \end{eqnarray} \]
Um dieses System von drei Differenzialgleichungen zweiter Ordnung numerisch lösen zu können, müssen wir es zunächst in ein System von sechs Differenzialgleichungen erster Ordnung umschreiben. Details zum numerischen Lösung von Differenzialgleichungen finden Sie in der Vorlesung Einführung in die Programmierung für Studierende der Physik in den Unterpunkten Differenzialgleichungen: Numerische Lösung von Anfangswertproblemen und Systeme von gekoppelten Differenzialgleichungen und Differenzialgleichungen zweiter Ordnung. Wir definieren dazu formal \[ t(\tau) = u_1(\tau) \, , \quad r(\tau) = u_2(\tau) \, , \quad \phi(\tau) = u_3(\tau) \, , \quad \frac{d t(\tau)}{d\tau} = u_4(\tau) \, , \quad \frac{d r(\tau)}{d\tau} = u_5(\tau) \, , \quad \frac{d \phi(\tau)}{d\tau} = u_6(\tau) \] und schreiben unser DGL-System wie folgt um: \[ \begin{eqnarray} \frac{d t(\tau)}{d\tau} &=& \frac{d u_1(\tau)}{d\tau} = u_4(\tau) \\ \frac{d r(\tau)}{d\tau} &=& \frac{d u_2(\tau)}{d\tau} = u_5(\tau) \\ \frac{d \phi(\tau)}{d\tau} &=& \frac{d u_3(\tau)}{d\tau} = u_6(\tau) \\ \frac{d^2 t(\tau)}{d\tau^2} &=& \frac{d u_4(\tau)}{d\tau} = \frac{2 M u_{4} u_{5}}{u_{2} \left(2 M - u_{2}\right)} \\ \frac{d^2 r(\tau)}{d\tau^2} &=& \frac{d u_5(\tau)}{d\tau} = - \frac{M u_{5}^{2}}{u_{2} \left(2 M - u_{2}\right)} + \frac{M u_{4}^{2} \left(2 M - u_{2}\right)}{u_{2}^{3}} + u_{6}^{2} \left(- 2 M + u_{2}\right) \\ \frac{d^2 \phi(\tau)}{d\tau^2} &=& \frac{d u_6(\tau)}{d\tau} = - \frac{2 u_{5} u_{6}}{u_{2}} \end{eqnarray} \]
Im Unterpunkt Abgeleitete Klassen, Vererbung von Klassenmerkmalen und Klassenhierarchien der Vorlesung hatten wir mittels der abstrakten Schnittstellen-Basisklasse 'dsolve' ein allgemeines System von Differenzialgleichungen numerisch mittels des Runge-Kutta Ordnung vier Verfahrens gelöst. Diese Basisklasse verwendeten wir dann z.B. zum numerischen Lösen der Planetenbewegungen (siehe Planetenbewegungen und das n-Körper-Problem) und das dort beschriebene C++ Programm verwenden wir nun auch, als Grundgerüst, zum Lösen der Geodätengleichung.
In dem folgenden C++ Header File GeodesicSchwarzschild.hpp ist die Subklasse 'geodesic' von der abstrakten Klasse 'dsolve' abgeleitet. Mittels der Subklasse 'geodesic' lassen sich Simulationen der Geodätengleichung mit unterschiedlichen Massenwerten und Anfangsbedingungen erstellen. Die Subklasse überschreibt dabei die rein virtuelle Funktion 'dgls' und definiert somit die Bewegungsgleichung der Geodätengleichung.
Mittels des Konstruktors der 'geodesic'-Klasse ist es nun möglich, eine Instanz eines 'geodesic'-Objektes zu erstellen. Der Konstruktor besitzt dabei fünf Argumente ( 'geodesic'(Masse des schwarzen Lochs M, Anfangswert der Eigenzeit a, Endwert der Eigenzeit b, Anzahl der Punkte N, Anfangswerte des Probekörpers u_init) ) und mit ihm kann der Benutzer die Bewegung eines Probekörpers in der Schwarzschild-Raumzeit simulieren.
/* Header-Datei * Siehe VARTC 2026/27 (siehe http://itp.uni-frankfurt.de/~hanauske/VARTC/index.html , Vorlesung 6) * basierend auf VPROG 2026 (siehe https://itp.uni-frankfurt.de/~hanauske/VPROG/index.html , Vorlesung 8,9,10) * C++ Programm zum Loesen der Geodaetengleichung * Berechnung der Lösung eines Systems von sechs Differenzialgleichungen erster Ordnung * in vorgegebener Schwarzschild-Raumzeit (nicht-rotierendes, ungeladenes Schwarzes Loch) * Die Loesung benutzt die Runge-Kutta Ordnung vier Methode (RK4) * Koordinaten (Koordinatenzeit: t=y_1 , Radius: r=y_2, Winkel: theta=pi/2=const, Winkel: phi=y_3) * Die zeitliche Entwicklung erfolgt in Eigenzeit tau von y_1(tau), y_2(tau), y_3(tau) * 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) 'geodesic' * Zeitentwicklung fuer unterschiedliche Eigenzeit-Werte tau in [a,b] * Unterklasse 'geodesic' * Konstruktor mit fuenf Argumenten: 'geodesic'(Masse des schwarzen Lochs M, Anfangswert der Eigenzeit a, Endwert der Eigenzeit b, Anzahl der Punkte N, Anfangswerte u_init) * u_init = alpha = {t,r,phi,d_t,d_r,d_phi} */ #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 (Eigenzeit) double b; // Obergrenze des Intervalls (Eigenzeit) int N; // Anzahl der Punkte vector<double> alpha; // Anfangswerte vector<vector<double>> ym; // Lösungsmatrix vector<double> tau; // Eigenzeit-Vektor // Virtuelle Funktion dgls, wird an anderer Stelle definiert virtual vector<double> dgls(double tau, 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_) { tau.push_back(a_); // Zum Eigenzeit-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 dtau zwischen den aequidistanten Punkten des tau-Intervalls (h=dtau) 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 tau-Intervalls for (int j = 0; j < n; ++j) k1[j] = h * dgls(tau[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(tau[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(tau[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(tau[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 tau.push_back(a + (i + 1) * h); // Zum Eigenzeit-Vektor die neuen Wert eintragen } // Ende for-Schleife ueber die einzelnen Punkte des tau-Intervalls } // Ende des Methode 'solve()' }; // Ende der Klasse 'dsolve' // Sub-Klasse 'geodesic' class geodesic : public dsolve { public: double M; // Masse M des Schwarzen Lochs // Überschreibung der virtuellen Funktion dgls, Definition der Geodätengleichung vector<double> dgls(double tau, vector<double> u_vec){ vector<double> du_dt(6); // Die sechs DGLs erster Ordnung werden als Standard-Vektor definiert du_dt[0] = u_vec[3]; // Sechs DGLs erster Ordnung du_dt[1] = u_vec[4]; // ... du_dt[2] = u_vec[5]; du_dt[3] = 2*M*u_vec[3]*u_vec[4]/(u_vec[1]*(2*M - u_vec[1])); du_dt[4] = -M*pow(u_vec[4], 2)/(u_vec[1]*(2*M - u_vec[1])) + M*pow(u_vec[3], 2)*(2*M - u_vec[1])/pow(u_vec[1], 3) + pow(u_vec[5], 2)*(-2*M + u_vec[1]); du_dt[5] = -2*u_vec[4]*u_vec[5]/u_vec[1]; return du_dt; // Rueckgabewert der DGLs als Standard-Vektor } // Ende: Definition der Bewegungsgleichung-Funktion // Konstruktor mit fuenf Argumenten: 'geodesic'(Masse M, Anfangswert der Eigenzeit a, Endwert der Eigenzeit b, Anzahl der Punkte N, Anfangswerte u_init) // u_init = alpha = {t,r,phi,d_t,d_r,d_phi} geodesic(double M_ = 1.0, double a_ = 0, double b_ = 800.0, int N_ = 10000, vector<double> alpha_ = {0.0,10.0,0.0,1.208356321620407,0.0,0.041}) : dsolve(a_, b_, N_, alpha_), M(M_) {} }; // Ende der Unterklasse 'geodesic'
Im Hauptprogramm GeodesicSchwarzschild.cpp (siehe unteres C++ Programm) wird dann mit dem Konstruktor eine Instanz der 'geodesic'-Klasse erstellt und die vererbte Funktion 'solve()' aufgerufen.
/* Siehe VARTC 2026/27 (siehe http://itp.uni-frankfurt.de/~hanauske/VARTC/index.html , Vorlesung 6) * basierend auf VPROG 2026 (siehe https://itp.uni-frankfurt.de/~hanauske/VPROG/index.html , Vorlesung 8,9,10) * C++ Programm zum Loesen der Geodaetengleichung * in vorgegebener Schwarzschild-Raumzeit (nicht-rotierendes, ungeladenes Schwarzes Loch) * Die Loesung benutzt die Runge-Kutta Ordnung vier Methode (RK4) * Koordinaten (Koordinatenzeit: t=y_1 , Radius: r=y_2, Winkel: theta=pi/2=const, Winkel: phi=y_3) * Die zeitliche Entwicklung erfolgt in Eigenzeit tau von y_1(tau), y_2(tau), y_3(tau) * Vewendet die Header-Datei GeodesicSchwarzschild.hpp * Konstruktor mit fuenf Argumenten: 'geodesic'(Masse M, Anfangswert der Eigenzeit a, Endwert der Eigenzeit b, Anzahl der Punkte N, Anfangswerte u_init) * u_init = alpha = {t,r,phi,d_t,d_r,d_phi} * Ausgabe zum Plotten mittels Python Jupyter Notebook GeodSchwarzschildII_OGRePy.ipynb * bzw. mittels des Python-Skrips GeodSchwarzschild.py */ #include "GeodesicSchwarzschild.hpp" // Header-Datei zum Loesen der Geodaetengleichung #include <iostream> // Standard Input- und Output Bibliothek #include <vector> // Vector-Container der Standardbibliothek using namespace std; // Benutze den Namensraum std // Hauptfunktion int main(){ double a = 0; // Untergrenze des Eigenzeit-Intervalls double b = 800.0; // Obergrenze des Eigenzeit-Intervalls int N = 100000; // Anzahl der Punkte FILE *ausgabe; // Deklaration eines Files fuer die Ausgabedatei ausgabe = fopen("geodesic.dat", "w+"); // Ergebnisse werden in die Datei 'geodesic.dat' geschrieben geodesic L = geodesic(1,a,b,N,{0.0,10.0,0.0,1.208356321620407,0.0,0.041}); // Konstruktor erzeugt eine Instanz der Klasse geodesic L.solve(); // Methode des Runge-Kutta-Verfahrens ausführen fprintf(ausgabe, "# 0: Index i \n# 1: tau \n# 2: t \n# 3:r \n# 4:phi \n# 5:d_t \n# 6:d_r \n# 7:d_phi \n"); // Beschreibung der ausgegebenen Groessen for(size_t i = 0; i < L.tau.size(); i += 10){ // for-Schleife über die Zeitgitterpunkte (nur jeden 10-ten Wert ausgeben) fprintf(ausgabe, "%3zu %19.15f ", i, L.tau[i]); // Ausgabe Index und Zeit for (const double value : L.ym[i]) { // Bereichsbasierte for-Schleife zur fprintf(ausgabe, "%19.15f ", value); // Ausgabe von Ort und Geschw. des Probekoerpers } // 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 'geodesic.dat' ausgegeben. Mittels des Python-Skriptes GeodesicSchwarzschild.py kann man dann die folgende Animation der C++-Simulation erstellen.
Die oben abgebildete Animation zeigt das effektive Potenzial als Funktion des Radius und die Trajektorie des Probekörpers in einem eingebetteten Raumzeit-Diagramm der Schwarzschild-Metrik.