# 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 # u
y1 = y # u'
y2 = y # v
y3 = y # v'
dy = np.zeros_like(y)
dy = y1
dy = y3

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

dy = -C[0,0]*dy**2 -\
2*C[0,0]*dy*dy -\
C[0,1]*dy**2
dy = -C[1,0]*dy**2 -\
2*C[1,0]*dy*dy -\
C[1,1]*dy**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 '15 at 12:29
• No deviation at all occurs along the equator or from pole to pole. Will append that to the question. Oct 24 '15 at 14:15
• Would you mind posting a small complete test case that can be run directly? Oct 25 '15 at 0:11
• Sure, but it will take me some time to minimize the code for this example only. Oct 25 '15 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 '15 at 7:26

The reason the resulting geodesic curve was deviating was because the calculated Christoffel symbol of second kind was incorrect. Using the correct Christoffel symbol :

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


results in the proper output (as displayed in the reference) : Now, I suppose I have to figure out why the calculated Christoffel symbol was incorrect. But that is another question.

• Thanks for your well crafted question and answer! My gut feeling for this bug was a coordinate system transformation error. Bug report to the Sympy folks? Oct 28 '15 at 7:31
• I checked the calculated Christoffel symbol for other coordinate system, and they were calculated correctly. The issue was only with the spherical transformation apparently. I remember now that I had hand coded the symbol. So I introduced the bug. The reason I had to hand code the symbol was because I could not produce it with my code alone. However, I will ask the Sympy google group about why I could not produce the Christoffel symbol for the sphere. Oct 28 '15 at 14:05
• I finally found the courage to ask the question as to why I get the wrong Christoffel symbol : stackoverflow.com/questions/33453941/… I will ask at the SymPy google group as well. Oct 31 '15 at 16:59