# Numerical solution of the geodesic equation of the Schwarzschild 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 Schwarzschild metric and numerically computing some geodesics.

Create now the coordinate system to define the metric as:

$c^2 \diff{\tau^2} = \left(1 - \frac{r_s}{r}\right) c^2 \diff{t^2} - \left(1 - \frac{r_s}{r}\right)^{-1} \diff{r^2} - r^2 \left(\diff{\theta^2} + \sin^2(\theta) \diff{\phi^2} \right)$

$\Gamma^0_{0,1} = \frac{r_{s}}{2 \mathbf{r}^{2} \left(1 - \frac{r_{s}}{\mathbf{r}}\right)}$

$\Gamma^0_{1,0} = \frac{r_{s}}{2 \mathbf{r}^{2} \left(1 - \frac{r_{s}}{\mathbf{r}}\right)}$

$\Gamma^1_{0,0} = - \frac{r_{s} \left(- c^{2} + \frac{c^{2} r_{s}}{\mathbf{r}}\right)}{2 \mathbf{r}^{2}}$

$\Gamma^1_{1,1} = \frac{c^{2} r_{s} \left(- c^{2} + \frac{c^{2} r_{s}}{\mathbf{r}}\right)}{2 \mathbf{r}^{2} \left(c^{2} - \frac{c^{2} r_{s}}{\mathbf{r}}\right)^{2}}$

$\Gamma^1_{2,2} = \frac{\mathbf{r} \left(- c^{2} + \frac{c^{2} r_{s}}{\mathbf{r}}\right)}{c^{2}}$

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

$\Gamma^2_{1,2} = \frac{1}{\mathbf{r}}$

$\Gamma^2_{2,1} = \frac{1}{\mathbf{r}}$

$\Gamma^2_{3,3} = - \sin{\left(\mathbf{\theta} \right)} \cos{\left(\mathbf{\theta} \right)}$

$\Gamma^3_{1,3} = \frac{1}{\mathbf{r}}$

$\Gamma^3_{2,3} = \frac{\cos{\left(\mathbf{\theta} \right)}}{\sin{\left(\mathbf{\theta} \right)}}$

$\Gamma^3_{3,1} = \frac{1}{\mathbf{r}}$

$\Gamma^3_{3,2} = \frac{\cos{\left(\mathbf{\theta} \right)}}{\sin{\left(\mathbf{\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),.$

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