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] -\
    dy[3] = -C[1,0][0]*dy[0]**2 -\
           2*C[1,0][1]*dy[0]*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).

enter image description here

Edit 1:

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

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.

  • $\begingroup$ Do you still see this divergence given a geodesic along the equator or one going form pole to pole? $\endgroup$
    – boyfarrell
    Commented Oct 24, 2015 at 12:29
  • $\begingroup$ No deviation at all occurs along the equator or from pole to pole. Will append that to the question. $\endgroup$
    – imranal
    Commented Oct 24, 2015 at 14:15
  • $\begingroup$ Would you mind posting a small complete test case that can be run directly? $\endgroup$
    – Kirill
    Commented Oct 25, 2015 at 0:11
  • $\begingroup$ Sure, but it will take me some time to minimize the code for this example only. $\endgroup$
    – imranal
    Commented Oct 25, 2015 at 8:39
  • $\begingroup$ 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. $\endgroup$
    – boyfarrell
    Commented Oct 26, 2015 at 7:26

1 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) :

enter image description here

Now, I suppose I have to figure out why the calculated Christoffel symbol was incorrect. But that is another question.

  • $\begingroup$ 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? $\endgroup$
    – boyfarrell
    Commented Oct 28, 2015 at 7:31
  • $\begingroup$ 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. $\endgroup$
    – imranal
    Commented Oct 28, 2015 at 14:05
  • 1
    $\begingroup$ 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. $\endgroup$
    – imranal
    Commented Oct 31, 2015 at 16:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.