### Jorge Cardona

I am a Mathematician and Engineer, eager to find the best intersection between both worlds. I am more than that really ...

# Numerical solution of the geodesic equation of the Schwarchild metric

sympy is a great tool for symbolic computations developed in Python. It has a module for differential geometry that allows you to compute some of the important objects as the Christoffel symbols and curvature objects. I wanted to give it a try and I created a jupyter notebook with a simple test using the Schwarzchild metric and numerically computing some geodesics.

In [1]:
%matplotlib inline

# General imports
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, Math

# Basic imports and functions
from sympy import latex, symbols, sin, cos, pi, simplify
from scipy.integrate import solve_ivp

from sympy.diffgeom import (
Manifold,
Patch,
CoordSystem,
metric_to_Christoffel_2nd,
TensorProduct as TP
)

def lprint(v):
display(Math(latex(v)))


Create now the coordinate system to define the metric as: $$c^2 d\tau^2 = \left(1 - \frac{r_s}{r}\right) \, c^2 \, d t^2 - \left(1 - \frac{r_s}{r}\right)^{-1}\, dr^2 - r^2 \left(d \theta^2 + \sin^2(\theta) d\phi^2 \right)$$

In [2]:
# Create a manifold.
M = Manifold('M', 4)

# Create a patch.
patch = Patch('P', M)

# Basic symbols
c, r_s = symbols('c r_s')

# Coordinate system
schwarzchild_coord = CoordSystem('schwarzchild', patch, ['t', 'r', 'theta', 'phi'])

# Get the coordinate functions
t, r, theta, phi = schwarzchild_coord.coord_functions()

# Get the base one forms.
dt, dr, dtheta, dphi = schwarzchild_coord.base_oneforms()

# Auxiliar terms for the metric.
dt_2 = TP(dt, dt)
dr_2 = TP(dr, dr)
dtheta_2 = TP(dtheta, dtheta)
dphi_2 = TP(dphi, dphi)
factor = (1 - r_s / r)

# Build the metric
metric = factor * c ** 2 * dt_2 - 1 / factor * dr_2 - r ** 2 * (dtheta_2 + sin(theta)**2 * dphi_2)
metric = factor * c ** 2 * dt_2 - 1 / factor * dr_2 - r ** 2 * (dtheta_2 + sin(theta)**2 * dphi_2)
metric = metric / c ** 2

In [3]:
# Get the Christoffel symbols of the second kind.
christoffel = metric_to_Christoffel_2nd(metric)

In [4]:
# Let's print this in an elegant way ;)
for i in range(4):
for j in range(4):
for k in range(4):
if christoffel[i, j, k] != 0:
display(Math('\Gamma^{}_{{{},{}}} = '.format(i, j, k) + latex(christoffel[i, j, k])))

$$\Gamma^0_{0,1} = \frac{r_{s}}{2 \boldsymbol{\mathrm{r}}^{2} \left(1 - \frac{r_{s}}{\boldsymbol{\mathrm{r}}}\right)}$$
$$\Gamma^0_{1,0} = \frac{r_{s}}{2 \boldsymbol{\mathrm{r}}^{2} \left(1 - \frac{r_{s}}{\boldsymbol{\mathrm{r}}}\right)}$$
$$\Gamma^1_{0,0} = - \frac{r_{s} \left(- c^{2} + \frac{c^{2} r_{s}}{\boldsymbol{\mathrm{r}}}\right)}{2 \boldsymbol{\mathrm{r}}^{2}}$$
$$\Gamma^1_{1,1} = \frac{c^{2} r_{s} \left(- c^{2} + \frac{c^{2} r_{s}}{\boldsymbol{\mathrm{r}}}\right)}{2 \boldsymbol{\mathrm{r}}^{2} \left(c^{2} - \frac{c^{2} r_{s}}{\boldsymbol{\mathrm{r}}}\right)^{2}}$$
$$\Gamma^1_{2,2} = \frac{\boldsymbol{\mathrm{r}} \left(- c^{2} + \frac{c^{2} r_{s}}{\boldsymbol{\mathrm{r}}}\right)}{c^{2}}$$
$$\Gamma^1_{3,3} = \frac{\boldsymbol{\mathrm{r}} \left(- c^{2} + \frac{c^{2} r_{s}}{\boldsymbol{\mathrm{r}}}\right) \sin^{2}{\left (\boldsymbol{\mathrm{\theta}} \right )}}{c^{2}}$$
$$\Gamma^2_{1,2} = \frac{1}{\boldsymbol{\mathrm{r}}}$$
$$\Gamma^2_{2,1} = \frac{1}{\boldsymbol{\mathrm{r}}}$$
$$\Gamma^2_{3,3} = - \sin{\left (\boldsymbol{\mathrm{\theta}} \right )} \cos{\left (\boldsymbol{\mathrm{\theta}} \right )}$$
$$\Gamma^3_{1,3} = \frac{1}{\boldsymbol{\mathrm{r}}}$$
$$\Gamma^3_{2,3} = \frac{\cos{\left (\boldsymbol{\mathrm{\theta}} \right )}}{\sin{\left (\boldsymbol{\mathrm{\theta}} \right )}}$$
$$\Gamma^3_{3,1} = \frac{1}{\boldsymbol{\mathrm{r}}}$$
$$\Gamma^3_{3,2} = \frac{\cos{\left (\boldsymbol{\mathrm{\theta}} \right )}}{\sin{\left (\boldsymbol{\mathrm{\theta}} \right )}}$$

# Solving the geodesic equation.¶

The geodesic equation is $${d^{2}x^{\mu } \over ds^{2}}+\Gamma ^{\mu }{}_{\alpha \beta }{dx^{\alpha } \over ds}{dx^{\beta } \over ds}=0$$

The term $x^\mu$ is a path depending on $s$.

Let $u(s) = x(s)$ and $v(s) = \dot{u}(s)$, then we can translate the geodesic equation into the equations $$\dot{u}^\mu(s) = F^\mu_1(u(s),v(s))$$ and $$\dot{v}^\mu(s) = F^\mu_2(u(v), v(s))$$ with $$F^\mu_1(u(s),v(s)) = v^\mu(s)$$ and $$F^\mu_2(u(s), v(s)) = - \Gamma^\mu_{\alpha \beta} v^\alpha(s) v^\beta(s)$$

In [5]:
## Specify c and r_s
christoffel_ = christoffel.subs({c: 2, r_s: 1})

def F(t, y):
u = y[0:4]
v = y[4:8]

du = v
dv = [0, 0, 0, 0]
for i in range(4):
for j in range(4):
for k in range(4):
dv[i] -= christoffel_.subs({r: u[1]})[i,j,k] * v[j] * v[k]

return np.concatenate((du, dv))


Use the method solve_ivp of sympy to compute numerically the geodesic starting at the event $(t,r,\theta, \phi) = (0, 10, 0, 0)$ in the direction of the tangent vector $(1, 0, 0, 0)$

In [54]:
T = 50
sol = solve_ivp(F, [0, T], [0, 10, 0, 0, 1, 0, 0, 0], t_eval=np.linspace(0, T, T * 123 + 1))

In [55]:
plt.figure(figsize=(14, 6),)
plt.plot(sol_1.y[0], sol_1.y[1])
ax = plt.gca()
ax.axhline(1, color="red", ls='--', lw=1)
plt.grid()
plt.ylim((0, 11))
plt.xlabel('t')
plt.ylabel('r')

Out[55]:
Text(0,0.5,'r')
In [67]:
T = 200
sol = solve_ivp(F, [0, T], [0, 10, 0, 0, 1, 0, 0.04, 0], t_eval=np.linspace(0, T, T * 123 + 1))

plt.figure(figsize=(14, 14),)
ax = plt.subplot(111, projection='polar')
ax.plot(sol.y[2], sol.y[1])
ax.grid(True)