Numerical solution of Geodesic differential equations with Python

Question

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.

Answer 1

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.