Im Folgenden wird die Geodätengleichung in vorgegebener Schwarzschild Raumzeit betrachtet. Die Geodätengleichung beschreibt wie sich ein Probekörper (Masse Probekörper $m << M$ Masse schwarzes Loch) im Raum bewegt und sagt voraus, dass diese Bewegung sich stets entlang der kürzesten Kurve, in der durch die Metrik beschriebenen gekrümmten Raumzeit, vollzieht. Zunächst wird das Python-Modul "OGRePy" eingebunden, welches auf dem Modul SymPy basiert und symbolische Berechnungen in der Allgemeinen Relativitätstheorie relativ einfach möglich macht. Eine ausführliche Einführung in das Modul "OGRePy" finden sie unter dem folgenden Link: OGRePy-Documentation.
Definition der Koordinaten und der kovarianten Raumzeit-Metrik der Schwarzschildmetrik: $$ g_{\mu\nu}=\left( \begin{array}{ccc} 1-\frac{2\,M}{r} & 0 & 0 & 0\\ 0& -\frac{1}{1-\frac{2\,M}{r}}& 0&0 \\ 0& 0& -r^2& 0\\ 0& 0& 0& -r^2 \hbox{sin}^2(\theta)\end{array} \right)$$
import OGRePy as gr
from sympy import *
init_printing()
By Barak Shoshany (baraksh@gmail.com) (baraksh.com)
v2.0.0 (2026-05-11)
GitHub repository: https://github.com/bshoshany/OGRePy
Documentation: .ipynb, .pdf, .html
Um symbolische Rechnungen in einer gekrümmten Raumzeit machen zu können, müssen wir zunächst wieder das verwendete Koordinatensystem und die gekrümmte Metrik der Raumzeit (hier die Schwarzschild-Metrik) definieren.
from OGRePy.abc import t, phi, theta
r = gr.sym("r", nonnegative=True)
M = gr.sym("M", nonnegative=True)
Spherical = gr.Coordinates(t, r, theta, phi)
g = gr.Metric(
coords = Spherical,
components = gr.diag((1-2*M/r), -1/(1-2*M/r), -r**2, -r**2*gr.s.sin(theta)**2),
)
g.show()
$$g{}_{\mu}{}_{\nu}\Bigg|_{\left(t, r, \theta, \phi\right)} = \left(\begin{matrix}- \frac{2 M}{r} + 1 & 0 & 0 & 0\\0 & - \frac{1}{- \frac{2 M}{r} + 1} & 0 & 0\\0 & 0 & - r^{2} & 0\\0 & 0 & 0 & - r^{2} \sin^{2}{\left(\theta \right)}\end{matrix}\right)$$
Die Chistoffel Symbole (zweiter Art): $$ \Gamma^{\sigma}_{\mu \nu} = \frac{1}{2}g^{\sigma \alpha} \left( g_{\nu \alpha| \mu} + g_{\mu \alpha| \nu} - g_{\mu \nu| \alpha} \right)$$
Hier speziell $$ \Gamma^{1}_{1 1} = \Gamma^{r}_{r r}$$
g.christoffel().list(indices=(-1,-1,1))
$$\begin{align*} \Gamma{}_{t}{}_{t}{}^{r} &= \frac{M \left(2 M - r\right)}{r^{3}} \\ \Gamma{}_{t}{}_{r}{}^{t} &= \frac{M}{r \left(- 2 M + r\right)} \\ \Gamma{}_{r}{}_{t}{}^{t} = \Gamma{}_{r}{}_{r}{}^{r} &= \frac{M}{r \left(2 M - r\right)} \\ \Gamma{}_{r}{}_{\theta}{}^{\theta} = \Gamma{}_{r}{}_{\phi}{}^{\phi} = -\Gamma{}_{\theta}{}_{r}{}^{\theta} = -\Gamma{}_{\phi}{}_{r}{}^{\phi} &= - \frac{1}{r} \\ \Gamma{}_{\theta}{}_{\theta}{}^{r} &= - 2 M + r \\ \Gamma{}_{\theta}{}_{\phi}{}^{\phi} = -\Gamma{}_{\phi}{}_{\theta}{}^{\phi} &= - \frac{1}{\tan{\left(\theta \right)}} \\ \Gamma{}_{\phi}{}_{\phi}{}^{r} &= \left(- 2 M + r\right) \sin^{2}{\left(\theta \right)} \\ \Gamma{}_{\phi}{}_{\phi}{}^{\theta} &= \frac{\sin{\left(2 \theta \right)}}{2} \end{align*}$$
g.christoffel().components(indices=(-1,-1,1),coords=Spherical)[1,1,1]
Der Riemann Tensor: $$ R^{\mu}{}\!_{\nu \rho \sigma} = \frac{\partial \Gamma^{\mu}_{\nu \sigma}}{\partial x^{\rho}} - \frac{\partial \Gamma^{\mu}_{\nu \rho}}{\partial x^{\sigma}} + \Gamma^{\alpha}_{\nu \sigma}\Gamma^{\mu}_{\rho \alpha} - \Gamma^{\alpha}_{\nu \rho}\Gamma^{\mu}_{\sigma \alpha} $$
g.riemann().list(exact=True)
$$\begin{align*} R{}^{t}{}_{r}{}_{t}{}_{r} = -R{}^{t}{}_{r}{}_{r}{}_{t} &= \frac{2 M}{r^{2} \left(- 2 M + r\right)} \\ R{}^{t}{}_{\theta}{}_{t}{}_{\theta} = -R{}^{t}{}_{\theta}{}_{\theta}{}_{t} = R{}^{r}{}_{\theta}{}_{r}{}_{\theta} = -R{}^{r}{}_{\theta}{}_{\theta}{}_{r} &= - \frac{M}{r} \\ R{}^{t}{}_{\phi}{}_{t}{}_{\phi} = -R{}^{t}{}_{\phi}{}_{\phi}{}_{t} = R{}^{r}{}_{\phi}{}_{r}{}_{\phi} = -R{}^{r}{}_{\phi}{}_{\phi}{}_{r} &= - \frac{M \sin^{2}{\left(\theta \right)}}{r} \\ R{}^{r}{}_{t}{}_{t}{}_{r} = -R{}^{r}{}_{t}{}_{r}{}_{t} &= \frac{2 M \left(- 2 M + r\right)}{r^{4}} \\ R{}^{\theta}{}_{t}{}_{t}{}_{\theta} = -R{}^{\theta}{}_{t}{}_{\theta}{}_{t} = R{}^{\phi}{}_{t}{}_{t}{}_{\phi} = -R{}^{\phi}{}_{t}{}_{\phi}{}_{t} &= \frac{M \left(2 M - r\right)}{r^{4}} \\ R{}^{\theta}{}_{r}{}_{r}{}_{\theta} = -R{}^{\theta}{}_{r}{}_{\theta}{}_{r} = R{}^{\phi}{}_{r}{}_{r}{}_{\phi} = -R{}^{\phi}{}_{r}{}_{\phi}{}_{r} &= \frac{M}{r^{2} \left(- 2 M + r\right)} \\ R{}^{\theta}{}_{\phi}{}_{\theta}{}_{\phi} = -R{}^{\theta}{}_{\phi}{}_{\phi}{}_{\theta} &= \frac{2 M \sin^{2}{\left(\theta \right)}}{r} \\ R{}^{\phi}{}_{\theta}{}_{\theta}{}_{\phi} = -R{}^{\phi}{}_{\theta}{}_{\phi}{}_{\theta} &= - \frac{2 M}{r} \end{align*}$$
Der Ricci Tensor:
$$ R_{\mu \nu} = R^{\sigma}{}\!_{\mu \sigma \nu} = \frac{\partial \Gamma^{\sigma}_{\mu \nu}}{\partial x^{\sigma}} - \frac{\partial \Gamma^{\sigma}_{\mu \sigma}}{\partial x^{\nu}} + \Gamma^{\rho}_{\mu \nu}\Gamma^{\sigma}_{\rho \sigma} - \Gamma^{\sigma}_{\mu \rho}\Gamma^{\rho}_{\nu \sigma} $$
g.ricci_tensor()
$$R{}_{\mu}{}_{\nu}\Bigg|_{\left(t, r, \theta, \phi\right)} = \left(\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right)$$
Wir berechnen nun das infinitesimale Weglängenelement $ds$. Es ist wie folgt definiert
$$ds^2=g_{\mu\nu} dx^\mu dx^\nu \quad ,$$
wobei $dx^\mu = \left(dx^0,dx^1,dx^2,dx^3 \right) = \left(dt,dr,d\theta,d\phi \right)$ den Abstand zweier infinitesimal benachbarter Ereignisse (eine infinitesimale Raumzeitvariation) beschreibt.
Das Quadrat des Linienelment $ds^2 = g_{\mu\nu} dx^\mu dx^\nu$ berechnet sich zu
ds2 = g.line_element()
ds2
Geodätische Bewegung eines Probekörpers¶
Die Geodätengleichung lässt sich durch folgendes Variationsprinzip herleiten
\begin{equation} \int_{x^\mu(A)}^{x^\mu(B)} ds= \int_{x^\mu(A)}^{x^\mu(B)} \sqrt{g_{\mu\nu} dx^\mu dx^\nu}=\int_A^B \sqrt{g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}} d\lambda \rightarrow \hbox{Extremal} \end{equation}
, wobei man, bei festgehaltenen Anfangs- und Endpunkten ($x^\mu(A)$ und $x^\mu(B)$), nach den extremalsten Kurven $x^\mu(\lambda)$ in der Raumzeit sucht. Den Parameter $\lambda$ bezeichnet man als affiner Parameter, welcher z.B. die Eigenzeit $\tau$ des Körpers sein kann. Die Geodätengleichung ergibt sich dann mittels der Euler-Lagrange Gleichungen, wobei man als Lagrangedichte den Ausdruck $L = \sqrt{g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}}$, bzw. alternativ $L = g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}$ verwendet.
\begin{equation} \frac{d}{d\lambda} \left( \frac{\partial L}{\partial \frac{d x^\alpha}{d \lambda}} \right) - \frac{\partial L}{\partial x^\alpha} = 0 \quad \xrightarrow{\cdot \frac{1}{2}g^{\mu \alpha }} \quad \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 \end{equation}
$\Gamma^\mu_{\nu \rho}$ sind die Christoffel Symbole zweiter Art und $\lambda$ ein affiner Parameter (z.B. die Eigenzeit $\tau$ des Probekörpers). Die Geodätengleichung ist ein System gekoppelter Differenzialgleichungen, das wir uns im Folgenden näher betrachten werden.
Das Modul "OGRePy" verwendet die Lagrangedichte $L = g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}$ zur Berechnung der Geodätengleichung:
g.lagrangian()
$$L\Bigg|_{\left(t, r, \theta, \phi\right)} = \frac{\dot{r}^{2} r^{2} - \dot{t}^{2} \left(2 M - r\right)^{2} - r^{3} \left(2 M - r\right) \left(\dot{\phi}^{2} \sin^{2}{\left(\theta \right)} + \dot{\theta}^{2}\right)}{r \left(2 M - r\right)}$$
Die Geodätengleichung mittels der Lagrangedichte $L$ berechnet lautet:
g.geodesic_from_lagrangian()
$$0{}_{\mu}\Bigg|_{\left(t, r, \theta, \phi\right)} = \left(\begin{matrix}\partial_{\lambda} \left(\frac{\dot{t} \left(- 2 M + r\right)}{r}\right)\\\frac{- 4 M^{3} \dot{t}^{2} + 4 M^{2} \dot{\phi}^{2} r^{3} \sin^{2}{\left(\theta \right)} + 4 M^{2} \dot{\theta}^{2} r^{3} + 4 M^{2} \dot{t}^{2} r + 4 M^{2} r^{2} \partial_{\lambda} \left(\frac{\dot{r} r}{2 M - r}\right) - 4 M \dot{\phi}^{2} r^{4} \sin^{2}{\left(\theta \right)} - 4 M \dot{\theta}^{2} r^{4} - M \dot{r}^{2} r^{2} - M \dot{t}^{2} r^{2} - 4 M r^{3} \partial_{\lambda} \left(\frac{\dot{r} r}{2 M - r}\right) + \dot{\phi}^{2} r^{5} \sin^{2}{\left(\theta \right)} + \dot{\theta}^{2} r^{5} + r^{4} \partial_{\lambda} \left(\frac{\dot{r} r}{2 M - r}\right)}{r^{2} \left(4 M^{2} - 4 M r + r^{2}\right)}\\\frac{\dot{\phi}^{2} r^{2} \sin{\left(2 \theta \right)}}{2} + \partial_{\lambda} \left(- \dot{\theta} r^{2}\right)\\\partial_{\lambda} \left(- \dot{\phi} r^{2} \sin^{2}{\left(\theta \right)}\right)\end{matrix}\right)$$
Wir benutzen im Folgenden jedoch die direkte Berechnung der Geodätengleichung mittels der Christoffel Symbole.
g.geodesic_from_christoffel()
$$0{}^{\mu}\Bigg|_{\left(t, r, \theta, \phi\right)} = \left(\begin{matrix}\frac{- 2 M \dot{r} \dot{t} + \ddot{t} r \left(2 M - r\right)}{r \left(2 M - r\right)}\\\frac{M \dot{r}^{2} r^{2} - M \dot{t}^{2} \left(2 M - r\right)^{2} + r^{3} \left(2 M - r\right) \left(\ddot{r} + \dot{\phi}^{2} \left(2 M - r\right) \sin^{2}{\left(\theta \right)} + \dot{\theta}^{2} \left(2 M - r\right)\right)}{r^{3} \left(2 M - r\right)}\\\ddot{\theta} - \frac{\dot{\phi}^{2} \sin{\left(2 \theta \right)}}{2} + \frac{2 \dot{\theta} \dot{r}}{r}\\\ddot{\phi} + \frac{2 \dot{\phi} \dot{\theta}}{\tan{\left(\theta \right)}} + \frac{2 \dot{\phi} \dot{r}}{r}\end{matrix}\right)$$
Die Geodätengleichung ist somit ein System von nichtlinearen gekoppelten Differenzialgleichungen zweiter Ordnung. Die entsprechenden vier Gleichungen lauten:
GeodEq = g.geodesic_from_christoffel().components(indices=(1,),coords=Spherical)
GeodEq0 = Eq(0,simplify(GeodEq[0]))
GeodEq0
bzw. umgestellt:
param = Spherical.of_param()
dot = Spherical.of_param_dot()
ddot = Spherical.of_param_ddot()
GeodEq0_a = Eq(ddot[0],solve(GeodEq0,ddot[0])[0])
GeodEq0_a
GeodEq1 = Eq(0,simplify(GeodEq[1]))
GeodEq1_a = Eq(ddot[1],solve(GeodEq1,ddot[1])[0])
GeodEq1_a
GeodEq2 = Eq(0,simplify(GeodEq[2]))
GeodEq2_a = Eq(ddot[2],solve(GeodEq2,ddot[2])[0])
GeodEq2_a
GeodEq3 = Eq(0,simplify(GeodEq[3]))
GeodEq3_a = Eq(ddot[3],solve(GeodEq3,ddot[3])[0])
GeodEq3_a
Radiales Einfallen eines Probekörpers in ein schwarzes Loch¶
Wir lassen im Folgenden nur ebene Bewegungen zu ($\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$) und beschränken uns zusätzlich auf ein radial in das schwarze Loch einfallenden Probekörper ($\frac{d \phi}{d\tau}=0, \frac{d^2 \phi}{d\tau^2}=0$). Unter diesen Annahmen sind nur die ersten beiden Geodätengleichungen relevant und diese lauten dann:
GeodEq0_a
OGRePy: You have the latest version of the package.
GeodEq1_a = GeodEq1_a.subs({(param[2],pi/2),(dot[2],0),(dot[3],0)})
GeodEq1_a
Es ist manchmal hilfreich auch eine Latex-Ausgabe der DGL zu erzeugen:
print(latex(GeodEq0_a))
print(latex(GeodEq1_a))
\frac{d^{2}}{d \lambda^{2}} t{\left(\lambda \right)} = \frac{2 M \frac{d}{d \lambda} r{\left(\lambda \right)} \frac{d}{d \lambda} t{\left(\lambda \right)}}{\left(2 M - r{\left(\lambda \right)}\right) r{\left(\lambda \right)}}
\frac{d^{2}}{d \lambda^{2}} r{\left(\lambda \right)} = \frac{M \left(2 M - r{\left(\lambda \right)}\right)^{2} \left(\frac{d}{d \lambda} t{\left(\lambda \right)}\right)^{2} - M r^{2}{\left(\lambda \right)} \left(\frac{d}{d \lambda} r{\left(\lambda \right)}\right)^{2}}{\left(2 M - r{\left(\lambda \right)}\right) r^{3}{\left(\lambda \right)}}
Um dieses System von zwei Differenzialgleichungen zweiter Ordnung mit Python numerisch lösen zu können, müssen wir es zunächst in ein System von vier 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) = y_1(\tau) \, , \quad r(\tau) = y_2(\tau) \, , \quad \frac{d t(\tau)}{d\tau} = y_3(\tau) \, , \quad \frac{d r(\tau)}{d\tau} = y_4(\tau) $$
und schreiben unser DGL-System wie folgt um:
$$ \begin{eqnarray} \frac{d^2 t(\tau)}{d\tau^2} &=& \frac{d y_3(\tau)}{d\tau} = f_3(\tau,y_1,y_2,y_3,y_4) = \frac{2 M \frac{d}{d \tau} r{\left(\tau \right)} \frac{d}{d \tau} t{\left(\tau \right)}}{r \left(2 M - r\right)} \\ \frac{d^2 r(\tau)}{d\tau^2} &=& \frac{d y_4(\tau)}{d\tau} = f_4(\tau,y_1,y_2,y_3,y_4) = \\ &=& \frac{M \left(- 4 M^{2} \left(\frac{d}{d \tau} t{\left(\tau \right)}\right)^{2} + 4 M r \left(\frac{d}{d \tau} t{\left(\tau \right)}\right)^{2} + r^{2} \left(\frac{d}{d \tau} r{\left(\tau \right)}\right)^{2} - r^{2} \left(\frac{d}{d \tau} t{\left(\tau \right)}\right)^{2}\right)}{r^{3} \left(- 2 M + r\right)} \end{eqnarray} $$
y1, y2, y3, y4 = symbols('y_1, y_2, y_3, y_4')
SysDGL_0 = Eq(ddot[0],GeodEq0_a.rhs.subs({(param[1],y2),(dot[0],y3),(dot[1],y4)}))
SysDGL_0
SysDGL_1 = Eq(ddot[1],GeodEq1_a.rhs.subs({(param[1],y2),(dot[0],y3),(dot[1],y4)}))
SysDGL_1
Diese beiden DGLs werden wir später zusätzlich mittels der Programmiersprache C++ lösen. Wir exportieren deshalb schon jetzt das zugrundeliegende System von DGLs als C++ Code:
from sympy.codegen.ast import Assignment
wert = symbols('wert')
print(ccode(Assignment(wert,SysDGL_0.rhs)))
print(ccode(Assignment(wert,SysDGL_1.rhs)))
wert = 2*M*y_3*y_4/(y_2*(2*M - y_2)); wert = (-M*pow(y_2, 2)*pow(y_4, 2) + M*pow(y_3, 2)*pow(2*M - y_2, 2))/(pow(y_2, 3)*(2*M - y_2));
Es ist manchmal hilfreich auch eine Latex-Ausgabe der DGL zu erzeugen:
print(latex(SysDGL_0))
print(latex(SysDGL_1))
\frac{d^{2}}{d \lambda^{2}} t{\left(\lambda \right)} = \frac{2 M y_{3} y_{4}}{y_{2} \left(2 M - y_{2}\right)}
\frac{d^{2}}{d \lambda^{2}} r{\left(\lambda \right)} = \frac{- M y_{2}^{2} y_{4}^{2} + M y_{3}^{2} \left(2 M - y_{2}\right)^{2}}{y_{2}^{3} \left(2 M - y_{2}\right)}
Nun lösen wir die DGLs mittels Python und setzen im Folgenden $M=1$:
f_3 = lambdify((y1,y2,y3,y4), SysDGL_0.rhs.subs(M,1))
f_4 = lambdify((y1,y2,y3,y4), SysDGL_1.rhs.subs(M,1))
Wir berechnen nun die numerische Lösung des oben beschriebenen Systems bestehend aus zwei gekoppelten DGLs.
Wir benutzen hierfür das Python-Modul SciPy, welches eine breite Kollektion von mathematischen Algorithmen und Funktionen bereitstellt. Im Speziellen werden wir in diesem Semester die Funktion 'solve_ivp(...)' verwenden, die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differenzialgleichungen bereitstellt. Zusätzlich berechnen wir die numerische Lösung auch mittels der älteren Methode "integrate.odeint()", die ebenfalls in dem Python-Modul "scipy" definiert ist.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
from scipy import integrate
Wir definieren das System der gekoppelten vier Differenzialgleichungen erster Ordnung:
def DGLsys(time, vy):
y1,y2,y3,y4 = vy
dy1 = y3
dy2 = y4
dy3 = f_3(y1,y2,y3,y4)
dy4 = f_4(y1,y2,y3,y4)
return np.array([dy1,dy2,dy3,dy4])
Zusätzlich auch für die ältere Methode "integrate.odeint()"
def DGLsys_odeint(vy, time):
y1,y2,y3,y4 = vy
dy1 = y3
dy2 = y4
dy3 = f_3(y1,y2,y3,y4)
dy4 = f_4(y1,y2,y3,y4)
return np.array([dy1,dy2,dy3,dy4])
Zum Vergleich lösen wir auch die Bewegungsgleichung nach Newton:
$$\frac{d^2r}{dt^2}=-\frac{1}{r^2}$$
def DGLsysNewton(time, vy):
r, y2 = vy
dr = y2
dy2 = -1/r**2
return np.array([dr,dy2])
Zusätzlich auch für die ältere Methode "integrate.odeint()"
def DGLsysNewton_odeint(vy, time):
r, y2 = vy
dr = y2
dy2 = -1/r**2
return np.array([dr,dy2])
In Abhängigkeit von den Anfangswerten können unterschiedliche Bahnen der Bewegung entstehen. Wir wählen zunächst als Beispiel die Anfangswerte eines Probekörpers, der sich zur Zeit $t=\tau=0$ bei einem Radius von $r_0=10=5 R_S$ befindet und eine Anfangsgeschwindigkeit in radialer Richtung von $\left. \frac{dr}{d\tau} \right|_{\tau=0}=0$ hat. Die weitere nötige Anfangsbedingung ( $\left. \frac{dt}{d\tau} \right|_{\tau=0}$ ) erhalten wir mittels der Bedingung $\left( \frac{ds}{d\tau} \right)^2= g_{\mu\nu} \frac{dx^\mu}{d\tau} \frac{dx^\nu}{d\tau} = g_{\mu\nu} u^\mu u^\nu = 1$, die für massive Probekörper immer erfüllt sein muss ($u^\mu$ ist die Vierergeschwindigkeit des Probekörpers). Wir beschreiben die Bewegung aus der Sichtweise eines im Unendlichen ruhenden Beobachters (Schwarzschild Koordinaten).
N = 10000
tval = np.linspace(0, 33.7, N+1)
tvalNewton = np.linspace(0, 35.12, N+1)
t0 = 0.0
r0 = 10.0
dr0 = 0.0
dt0 = 1/np.sqrt(1-2/r0)
initialval = np.array([t0,r0,dt0,dr0])
initialvalNewton = np.array([r0,dr0])
Loes = integrate.solve_ivp(DGLsys, [tval[0], tval[-1]], initialval, t_eval=tval)
LoesNewton = integrate.solve_ivp(DGLsysNewton, [tvalNewton[0], tvalNewton[-1]], initialvalNewton, t_eval=tvalNewton)
Loes_odeint = integrate.odeint(DGLsys_odeint, initialval, tval)
LoesNewton_odeint = integrate.odeint(DGLsysNewton_odeint, initialvalNewton, tvalNewton)
rcParams.update({
'figure.figsize' : [10,8],
'axes.titlesize' : 14,
'axes.labelsize' : 16,
'xtick.labelsize' : 14 ,
'ytick.labelsize' : 14
})
plt.cla()
plt.ylabel(r"$\rm Koordinaten \, Zeit \,\,t$")
plt.xlabel(r"$\rm Radius \,\,r$")
plt.xlim(0,10.1)
plt.ylim(0,57)
plt.plot(Loes.y[1],Loes.y[0],c="blue", linewidth=1.5, linestyle='-',label="Allgemeine Relativitätstheorie");
plt.plot(LoesNewton.y[0],tval,c="red", linewidth=1.5, linestyle='-',label="Newton");
plt.legend(loc='upper right', fontsize=18)
plt.axvline(x=2, linewidth=1.5, color="black", linestyle=':')
plt.text(2, 5, r"$R_S=2M$", size=20, rotation=0.,ha="center", va="center",
bbox=dict(boxstyle="round",ec="black",fc="white"));
Zusätzlich auch für die ältere Methode "integrate.odeint()"
plt.cla()
plt.ylabel(r"$\rm Koordinaten \, Zeit \,\,t$")
plt.xlabel(r"$\rm Radius \,\,r$")
plt.xlim(0,10.1)
plt.ylim(0,57)
plt.plot(Loes_odeint[:, 1],Loes_odeint[:, 0],c="blue", linewidth=1.5, linestyle='-',label="Allgemeine Relativitätstheorie");
plt.plot(LoesNewton_odeint[:, 0],tval,c="red", linewidth=1.5, linestyle='-',label="Newton");
plt.legend(loc='upper right', fontsize=18)
plt.axvline(x=2, linewidth=1.5, color="black", linestyle=':')
plt.text(2, 5, r"$R_S=2M$", size=20, rotation=0.,ha="center", va="center",
bbox=dict(boxstyle="round",ec="black",fc="white"));
Der in das schwarze Loch fallende Probekörpers bleibt somit, für einen äußeren Beobachter, anscheinend am Ereignishorizont bei $R_S=2M$ stehen (siehe die blaue Kurve im oberen Raumzeitdiagramm). Seine Bewegung in radialer Richtung wird nahe dem Ereignishorizont immer langsamer und das Bild, welches ein äußerer Beobachter wahrnimmt, friert am Ereignishorizont ein. Man kann zeigen, dass das wahrgenommene Bild gleichzeitig unendlich rotverschoben wird, sodass es nahe dem Ereignishorizont, langsam den wahrnehmbaren, sichtbaren Bereich verlässt.
Eine gute und zugleich lehrreiche Verbildlichung dieser und weiterer seltsamen Eigenschaften von schwarzen Löchern findet sich in dem Essay Black Holes and the German Reichstag (Schwarze Löcher und der Deutsche Reichstag) welcher 2005 im Magazin Physics World erschien (siehe Black Holes and the German Reichstag).
Die untere Abbildung veranschaulicht zusätzlich zur Lösung des radialen Falls des Probekörpers einige radial nach Außen gesendete Lichtstrahlen (gepunktete rote Kurven).
plt.cla()
plt.ylabel(r"$\rm Koordinaten \, Zeit \,\,t$")
plt.xlabel(r"$\rm Radius \,\,r$")
plt.xlim(0,10.1)
plt.ylim(0,80)
plt.plot(Loes.y[1],Loes.y[0],c="blue", linewidth=1.5, linestyle='-',label="Allgemeine Relativitätstheorie");
positions = [8, 6, 4, 3, 2.3, 2.1, 2.05, 2.01]
for pos in positions:
rval = np.linspace(pos, 10, 1000)
it = np.where(Loes.y[1] < pos)
PhotonOut = 2*np.log(np.abs(rval-2))+rval-(2*np.log(np.abs(pos-2))+pos)+Loes.y[0][it[0][0]]
plt.plot(rval, PhotonOut, c="red", linewidth=1.5, linestyle=":")
plt.legend(loc='lower center', fontsize=18)
plt.axvline(x=2, linewidth=1.5, color="black", linestyle=':')
plt.text(2, 15, r"$R_S=2M$", size=20, rotation=0.,ha="center", va="center",
bbox=dict(boxstyle="round",ec="black",fc="white"));
Zusätzlich auch für die ältere Methode "integrate.odeint()"
plt.cla()
plt.ylabel(r"$\rm Koordinaten \, Zeit \,\,t$")
plt.xlabel(r"$\rm Radius \,\,r$")
plt.xlim(0,10.1)
plt.ylim(0,80)
plt.plot(Loes_odeint[:, 1],Loes_odeint[:, 0],c="blue", linewidth=1.5, linestyle='-',label="Allgemeine Relativitätstheorie");
positions = [8, 6, 4, 3, 2.3, 2.1, 2.05, 2.01]
for pos in positions:
rval = np.linspace(pos, 10, 1000)
it = np.where(Loes_odeint[:, 1] < pos)
PhotonOut = 2*np.log(np.abs(rval-2))+rval-(2*np.log(np.abs(pos-2))+pos)+Loes_odeint[it[0][0],0]
plt.plot(rval, PhotonOut, c="red", linewidth=1.5, linestyle=":")
plt.legend(loc='lower center', fontsize=18)
plt.axvline(x=2, linewidth=1.5, color="black", linestyle=':')
plt.text(2, 15, r"$R_S=2M$", size=20, rotation=0.,ha="center", va="center",
bbox=dict(boxstyle="round",ec="black",fc="white"));
