# Numerical solution of Geodesic differential equations with Python

I have made a solver based on the SymPy.diffgeom library, where I use Scipy.Integrate to solve the following system of second order differential equations :

\begin{align} u'' &+ \Gamma^0_{00}(u')^2 + 2\Gamma^0_{01}u'v' + \Gamma^0_{11}(v')^2 = 0,\\ v'' &+ \Gamma^1_{00}(u')^2 + 2\Gamma^1_{01}u'v' + \Gamma^1_{11}(v')^2 = 0, \end{align} where $\Gamma^m_{ij}(u(s),v(s))$ is the Christoffel symbol of second kind. The geodesic solution $u = u(s),v = v(s)$ is a curve defined for the interval $s\in[s_0,s_1]$. These equations can be described as a system of first order differential equations by setting $p = u'$, and $q = v'$ : \begin{align} p' &+ \Gamma^0_{00}p^2 + 2\Gamma^0_{01}pq + \Gamma^0_{11}q^2 = 0\\ q' &+ \Gamma^1_{00}p^2 + 2\Gamma^1_{01}pq + \Gamma^1_{11}q^2 = 0 \end{align} with the initial condition given as $u(s_0) = u_0, u'(s_0) = du_0, v(s_0) = v_0, v'(s_0) = dv_0$.

Here is a code snippet of my implementation

def f(y,s,C,u,v):
y0 = y[0] # u
y1 = y[1] # u'
y2 = y[2] # v
y3 = y[3] # v'
dy = np.zeros_like(y)
dy[0] = y1
dy[2] = y3

C = C.subs({u:y0,v:y2}) # Evaluate C for u,v = (u0,v0)

dy[1] = -C[0,0][0]*dy[0]**2 -\
2*C[0,0][1]*dy[0]*dy[2] -\
C[0,1][1]*dy[2]**2
dy[3] = -C[1,0][0]*dy[0]**2 -\
2*C[1,0][1]*dy[0]*dy[2] -\
C[1,1][1]*dy[2]**2
return dy

def solve(C,u0,s0,s1,ds):
s = np.arange(s0,s1+ds,ds)
# The Christoffel symbol of 2nd kind, C, is a function of (u,v)
from sympy.abc import u,v
return sc.odeint(f,u0,s,args=(C,u,v))  # integration method : LSODA


I have implemented several generic test cases : torus, sphere, egg carton, and catenoid. However, there seem to be some issue with the solver. On a sphere, for example, the geodesic curve is the great circle (see reference). When I try to find the geodesic curve and plot on a sphere (with same parmaters as reference provided), the curve starts to veer off. There seem to be some sort of numerical instability that is altering the course of the geodesic curve over the interval $s\in[s_0,s_1]$. Is there any way I can make my solver more stable ? I have tried to reduce the step-size of the solver, but that has not made things any better (visually at least...I could probably try to estimate the convergence rate).

Edit 1:

I checked the code for the geodesic curve along the equator, and from pole to pole, and I get no deviation :

Edit 2:

I have pasted the test case on the following link : Geodesic on a sphere.

I forgot to mention this, but I am using the following versions

• SymPy : 0.7.7.dev
• SciPy : 0.16.0

It should now be possible for anyone to reproduce the same results.

• Do you still see this divergence given a geodesic along the equator or one going form pole to pole? Oct 24, 2015 at 12:29
• No deviation at all occurs along the equator or from pole to pole. Will append that to the question. Oct 24, 2015 at 14:15
• Would you mind posting a small complete test case that can be run directly? Oct 25, 2015 at 0:11
• Sure, but it will take me some time to minimize the code for this example only. Oct 25, 2015 at 8:39
• If it is true that start point anywhere leads back to the originating point then seems like you have a bug. I asked you to check out equatorial and polar trajectories because in each of these one of the angular dimensions ($theta$, $phi$) is always zero. Maybe this hints at a coordinate system error somewhere? Hard to say exactly. Oct 26, 2015 at 7:26

C =  Matrix([[(0, -tan(v)), (0,0)],[(sin(v)*cos(v),0),(0, 0)]])