Allgemeine Relativitätstheorie mit dem Computer

General Theory of Relativity on the Computer

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

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

Frankfurt am Main 04.04.2021

Erster Vorlesungsteil: Allgemeine Relativitätstheorie mit Python

Eine kleine Einführung in Jupyter Notebooks

Einführung

In Jupyter Notebooks kann man die Programmiersprache Python in einer anwendungsfreundlichen Umgebung nutzen und die berechneten Ergebnisse auch gleich visualisieren. Sie sind somit in ihrem Erscheinungsbild den kommerziellen Software-Produkten Mathematica oder Maple sehr ähnlich.

Die in dieser Vorlesung bereitgestellten Jupyter Notebooks benutzen als Kernel Python 3. In diesem Notebook werden einige Python Module vorgestellt, die wir im Laufe der Vorlesung verwenden werden. Zunächst wird das Python Modul "sympy" eingebunden, das ein Computer-Algebra-System für Python bereitstellt und symbolische Berechnungen und im speziellen Matrix/Tensor-Berechnungen relativ einfach möglich macht. Falls Sie das "sympy" Modul das erste Mal verwenden, müssen Sie es zunächst in Ihrer Python 3 Umgebung installieren (z.B. in einem Linux Terminal mit "pip3 install sympy").

In [1]:
from sympy import *
init_printing()

Wir definieren uns eine (2x3)-Matrix $\hat{\bf {A}}$ und eine (3x3)-Matrix $\hat{\bf {B}}$:

In [2]:
MA = Matrix([[-1,3,9],[2,-5,6]])
MA
Out[2]:
$$\left[\begin{matrix}-1 & 3 & 9\\2 & -5 & 6\end{matrix}\right]$$
In [3]:
MB = Matrix([[-7,-1,5],[-9,-3,4],[1,8,-2]])
MB
Out[3]:
$$\left[\begin{matrix}-7 & -1 & 5\\-9 & -3 & 4\\1 & 8 & -2\end{matrix}\right]$$

Das Produkt der Matrizen $\hat{\bf {A}}$ und $\hat{\bf {B}}$ berechnet sich zu

In [4]:
MA*MB
Out[4]:
$$\left[\begin{matrix}-11 & 64 & -11\\37 & 61 & -22\end{matrix}\right]$$

, die transponierte Matrix $(\hat{\bf {A}})^T$

In [5]:
transpose(MA)
Out[5]:
$$\left[\begin{matrix}-1 & 2\\3 & -5\\9 & 6\end{matrix}\right]$$

und die Determinante der Matrix $\hat{\bf {B}}$ berechnet sich zu:

In [6]:
MB.det()
Out[6]:
$$-149$$

Definition von diagonalen Matrizen:

In [7]:
MC=diag(1, -1, -1, -1)
MC
Out[7]:
$$\left[\begin{matrix}1 & 0 & 0 & 0\\0 & -1 & 0 & 0\\0 & 0 & -1 & 0\\0 & 0 & 0 & -1\end{matrix}\right]$$
In [8]:
MC.det()
Out[8]:
$$-1$$

Weitere Matrix Operationen finden Sie unter http://docs.sympy.org/latest/tutorial/matrices.html.

Das Modul sympy ist auch dafür gemacht symbolische Berechnungen durchzuführen. Betrachten wir z.B. die folgende Gleichung $0=f(x)=x^2-3x-9$. Funktionen (hier f(x)) und Variablen (hier x) müssen in sympy zunächst als Symbole definiert werden.

In [9]:
x = symbols('x')
f = Function('f')(x)
f = x**2-3*x-9
Gl = Eq(f, 0)
Gl
Out[9]:
$$x^{2} - 3 x - 9 = 0$$

Ein Funktionswert der Funktion f kann man wie folgt ausgeben (z.B. $f(2)$ ):

In [10]:
f.subs(x, 2)
Out[10]:
$$-11$$

Die Lösungen der Gleichung $f(x)=0$ lauten:

In [11]:
LoesGL=solve(Gl)
LoesGL
Out[11]:
$$\left [ \frac{3}{2} + \frac{3 \sqrt{5}}{2}, \quad - \frac{3 \sqrt{5}}{2} + \frac{3}{2}\right ]$$

bzw. als Zahlenwert:

In [12]:
(LoesGL[0].evalf() , LoesGL[1].evalf() )
Out[12]:
$$\left ( 4.85410196624968, \quad -1.85410196624968\right )$$

Überprüfen durch Einsetzen der Lösung (hier die erste, positive Lösung) in die Gleichung:

In [13]:
Gl.subs({(x,LoesGL[0])})
Out[13]:
$$\mathrm{True}$$

Darstellen der Funktion in einem xy-Diagramm:

In [15]:
plot(f,(x, -5, 8),xlabel="x",ylabel="f(x)");

Ableiten der Funktion nach x $\frac{df}{dx}$ :

In [16]:
f.diff(x)
Out[16]:
$$2 x - 3$$

Unbestimmtes Integral der Funktion $\int f(x)\,dx$

In [17]:
f.integrate(x)
Out[17]:
$$\frac{x^{3}}{3} - \frac{3 x^{2}}{2} - 9 x$$

Bestimmtes Integral der Funktion, z.B. $\int_{-4}^{6} f(x)\,dx$

In [18]:
f.integrate((x,-4,6))
Out[18]:
$$- \frac{80}{3}$$

Zusätzlich ist es in sympy auch möglich Differentialgleichungen (DGLs) zu lösen (falls eine analytische Lösung existiert). Betrachten wir z.B. die folgende DGL zweiter Ordnung: $\frac{d^2y(x)}{dx^2}=-3y(x)$

In [19]:
y = Function('y')(x)
DGL = Eq(y.diff(x,x), -3*y)
DGL
Out[19]:
$$\frac{d^{2}}{d x^{2}} y{\left (x \right )} = - 3 y{\left (x \right )}$$

Die allgemeine Lösung lautet:

In [20]:
dsolve(DGL)
Out[20]:
$$y{\left (x \right )} = C_{1} \sin{\left (\sqrt{3} x \right )} + C_{2} \cos{\left (\sqrt{3} x \right )}$$

Die allgemeine Lösung hängt von den zwei Parametern $C_1$ und $C_2$ ab, die durch eine Wahl der Anfangs/Randbedingungen festgelegt werden können. Wir wählen z.B.

In [21]:
Eq(y.subs(x,0),3)
Out[21]:
$$y{\left (0 \right )} = 3$$
In [22]:
Eq(y.diff(x).subs(x, 0),2)
Out[22]:
$$\left. \frac{d}{d x} y{\left (x \right )} \right|_{\substack{ x=0 }} = 2$$
In [23]:
Loes = dsolve(DGL,ics={y.subs(x,0):3,y.diff(x).subs(x, 0): 2})
Loes
Out[23]:
$$y{\left (x \right )} = \frac{2 \sqrt{3} \sin{\left (\sqrt{3} x \right )}}{3} + 3 \cos{\left (\sqrt{3} x \right )}$$

Darstellen der Lösung in einem xy-Diagramm:

In [24]:
plot(Loes.rhs,(x, 0 , 8),xlabel="x",ylabel="y(x)");

Es ist manchmal nötig die sympy-Ausdrücke (z.B. hier die Lösung der Differentialgleichung) in eine nummerische Funktion umzuwandeln. Diese Umwandlung wird mittels der sympy-Funktion "lambdify" gemacht. Die Lösung kann dann z.B. mit dem plot-Modul matplotlib visualisiert werden:

In [25]:
import matplotlib
import matplotlib.pyplot as plt 
import numpy as np

xval = np.linspace(0, 10, 1000)
func = lambdify(x, Loes.rhs)
plt.ylabel('y(x)')
plt.xlabel('x')
plt.plot(xval,func(xval));

Oft jedoch existiert eine analytische Lösung einer Differentialgleichung nicht und man muss die Lösung der Differentialgleichung nummerisch bestimmen. Wir zeigen im Folgenden wie man eine solche nummerische Lösung erhält und illustrieren dies am vorigen Beispiel: $\frac{d^2y(x)}{dx^2}=-3y(x)$

Um diese Differentialgleichung zweiter Ordnung mit Python nummerisch lösen zu können, müssen wir sie zunächst in ein System von zwei Differentialgleichungen erster Ordnung umschreiben. Wir definieren dazu formal: $y_1=\frac{dy}{dx}$, $y_2=\frac{dy_1}{dx}=\frac{d^2y}{dx^2}$.

Wir definieren das System der gekoppelten zwei Differentialgleichungen:

In [26]:
def DGLsys(vy, x):
    y, y1 = vy
    dy = y1
    dy1 = -3*y
    return np.array([dy,dy1])

Das nummerische Lösen der Differentialgleichung machen wir mit der im Modul "scipy" definierten Funktion integrate.odeint.

In [27]:
import numpy as np
import matplotlib.pyplot as plt 
from scipy import integrate
In [28]:
xval = np.linspace(0, 10, 10001)

x0=0.0
y0=3.0
dy0=2.0
initialval = np.array([y0,dy0])

LoesNum = integrate.odeint(DGLsys, initialval, xval)

Die Lösung können wir uns wieder grafisch darstellen:

In [29]:
plt.ylabel('y(x)')
plt.xlabel('x')
plt.plot(xval,LoesNum[:, 0]);

Die Funktionen, bzw. Variablendeklaration in sympy kann man sich anzeigen lassen mit:

In [30]:
type(x)
Out[30]:
sympy.core.symbol.Symbol

und löschen mit dem folgenden Befehl:

In [31]:
del(x,y)

Die Visualisierung von Funktionen die von mehreren Variablen abhängen kann man auf unterschiedlichen Weisen in Jupyter Notebooks realisieren. Betrachten wir z.B. die Funktion $f(x,y)={\rm sin}(x)+{\rm cos}(y)$

In [32]:
def f(x,y):
  val=np.sin(x) + np.cos(y)
  return val

Die Funktion stellt eine Fläche im Raum dar und wir visualisieren diese zunächst mittels matplotlib:

In [33]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib
from matplotlib import cm
In [34]:
X = np.linspace(0, 12, 300)
Y = np.linspace(0, 12, 300)
X, Y = np.meshgrid(X, Y)
ZF = np.array(list(f(X, Y)))
In [35]:
params = {
    'figure.figsize'    : [11,7],
#    'text.usetex'       : True,
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
}
matplotlib.rcParams.update(params)
In [36]:
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, ZF, cmap=cm.Blues, linewidth=0, alpha=1)
ax.view_init(azim=35, elev=55)
fig.colorbar(surf, shrink=0.5, aspect=10)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("f");

Man kann in Jupyter Notebooks auch animierte Grafiken erzeugen.

In der folgenden Animation wurde die Funktion $f(x,y,t)={\rm sin}(\frac{x}{\frac{t}{25}+1})+{\rm cos}(\frac{y}{\frac{t}{25}+1})$ in einem xy-Diagramm dargestellt und der Parameter $t$ im Bereich $t \in [0,20]$ variiert.

In [37]:
import matplotlib.animation as animation
from IPython.display import HTML

def f(x,y,t):
  val=np.sin(x/(t/25+1)) + np.cos(y/(t/25+1))
  return val

fig = plt.figure()
ax = fig.gca(projection='3d')

def init():
    ZF = np.array(list(f(X, Y,0)))
    ax.plot_surface(X, Y, ZF, cmap=cm.Blues, linewidth=0, alpha=1)
    return fig,

def animate(t):
    ax.cla()
    ZF = np.array(list(f(X, Y,t)))
    ax.plot_surface(X, Y, ZF, cmap=cm.Blues, linewidth=0, alpha=1)
    return fig,

ani = animation.FuncAnimation(fig,animate,init_func=init,frames=20,interval=50)
fig.colorbar(surf, shrink=0.5, aspect=10)
ax.view_init(azim=35, elev=55)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("f")

plt.close(ani._fig)

HTML(ani.to_html5_video())
Out[37]:

Die Visualisierung kann man auch mittels der "Plotly Python Open Source Graphing Library" erstellen. In dieser Bibliothek ist es möglich interaktive Abbildungen in Python darzustellen ( siehe https://plotly.com/python/ ).

In [38]:
import plotly.graph_objects as go

fig = go.Figure(go.Surface(x=X, y=Y, z=ZF, colorscale='Blues', showscale=True))

fig.update_layout(autosize=True,
                  width=700, height=700,
                  margin=dict(l=15, r=50, b=65, t=90))
fig.update_layout(scene = dict(
                    xaxis_title='x',
                    yaxis_title='y',
                    zaxis_title='f'))