I would like to compute the shortest distance between a number on pairs of points on a 2-dimensional Riemanniann statistical manifold (Negative Binomial distributions manifold, equipped with the Rao's metrics).

This is sometimes VERY long, because I suspect there is a cut point on the geodesic (at least in a numerical sense, say). So, I thought it was possible to save time by using NDSolveProcessEquations and NDSolveReinitialize...

Here are the Christoffel symbols, in one of the standard parametrizations:

Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Phi]][\[Phi]_, \[Mu]_] := \
-((\[Mu] (\[Mu] + 
          2 \[Phi]) + \[Phi]^2 (\[Mu] + \[Phi]) (-(\[Phi]/(\[Mu] + \
\[Phi]))^\[Phi] (\[Mu] + (\[Mu] + \[Phi]) Log[\[Phi]/(\[Mu] + \
\[Phi])]) PolyGamma[
            1, \[Phi]] - (\[Mu] + \[Phi]) (-1 + (\[Phi]/(\[Mu] + \
\[Phi]))^\[Phi]) PolyGamma[
            2, \[Phi]]))/(2 \[Phi] (\[Mu] + \[Phi]) (\[Mu] + \[Phi] (\
\[Mu] + \[Phi]) (-1 + (\[Phi]/(\[Mu] + \[Phi]))^\[Phi]) PolyGamma[
           1, \[Phi]])));
Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Mu]][\[Phi]_, \[Mu]_] := \
-((\[Phi] (-1 + \[Phi] (\[Phi]/(\[Mu] + \[Phi]))^\[Phi] (\[Mu] + \
\[Phi]) PolyGamma[
           1, \[Phi]]))/(2 (\[Mu] + \[Phi]) (\[Mu] + \[Phi] (\[Mu] + \
\[Phi]) (-1 + (\[Phi]/(\[Mu] + \[Phi]))^\[Phi]) PolyGamma[
           1, \[Phi]])));
Subscript[\[CapitalGamma], \[Phi], \[Mu]\[Mu]][\[Phi]_, \[Mu]_] := \
\[Phi]/(2 (\[Mu] + \[Phi]) (\[Mu] + \[Phi] (\[Mu] + \[Phi]) (-1 + (\
\[Phi]/(\[Mu] + \[Phi]))^\[Phi]) PolyGamma[1, \[Phi]]));
Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Phi]][\[Phi]_, \[Mu]_] := 
  1/2 \[Mu] (1/(\[Mu] \[Phi] + \[Phi]^2) - (\[Phi]/(\[Mu] + \[Phi]))^\
\[Phi] PolyGamma[1, \[Phi]]);
Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Mu]][\[Phi]_, \[Mu]_] := \
\[Mu]/(2 \[Phi] (\[Mu] + \[Phi]));
Subscript[\[CapitalGamma], \[Mu], \[Mu]\[Mu]][\[Phi]_, \[Mu]_] := -((
   2 \[Mu] + \[Phi])/(2 \[Mu] (\[Mu] + \[Phi])));

and the associated Euler-Lagrange equation (sorry, I tried to write something nice, see rather the attached notebook)

GeodE[{u_, v_}, t_, a_, b_] := 
  SetPrecision[{u''[
      t] == -(Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Phi]][u[t], 
          v[t]] (u'[t]^2) + 
        2 Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Mu]][u[t], 
          v[t]] u'[t] v'[t] + 
        Subscript[\[CapitalGamma], \[Phi], \[Mu]\[Mu]][u[t], 
          v[t]] (v'[t]^2)),
    v''[t] == -(Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Phi]][u[t], 
          v[t]] (u'[t]^2) + 
        2 Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Mu]][u[t], 
          v[t]] u'[t] v'[t] + 
        Subscript[\[CapitalGamma], \[Mu], \[Mu]\[Mu]][u[t], 
          v[t]] (v'[t]^2))
    , {u[0], v[0]} == a, {u[1], v[1]} == b}, Chif];

I consider a pair of points such that the equation above can be quickly solved by NDSolve.

{Easy4a, Easy4b} = {{0.0317527, 0.44814}, {1.46715, 9.50924}};

To determine geodesics between other pair of points, I defined (as in the documentation)

EuLa = First[
  NDSolve`ProcessEquations[{u''[
      t] == -(Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Phi]][u[t], 
          v[t]] (u'[t]^2) + 
        2 Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Mu]][u[t], 
          v[t]] u'[t] v'[t] + 
        Subscript[\[CapitalGamma], \[Phi], \[Mu]\[Mu]][u[t], 
          v[t]] (v'[t]^2)),
    v''[t] == -(Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Phi]][u[t], 
          v[t]] (u'[t]^2) + 
        2 Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Mu]][u[t], 
          v[t]] u'[t] v'[t] + 
        Subscript[\[CapitalGamma], \[Mu], \[Mu]\[Mu]][u[t], 
          v[t]] (v'[t]^2))
    , {u[0], v[0]} == {A[[1]], A[[2]]}, {u[1], v[1]} == {B[[1]], 
      B[[2]]}}, {u, v}, {t, 0, 1}]]

It seems that there is no problem here. Then I tried to find another geodesic, linking two other points:

B = TestC;
A = Centre; 
solBis = NDSolve`Reinitialize[
   Eula, {u[0], v[0]} == A, {u[1], v[1]} == B];

Total fiasco:

NDSolve::vdobj: Encountered Eula, which is not a valid NDSolve`StateData expression.

I'm using NDSolve for the first time, and I can't see where the problem is. Is the formuation of EuLa ill-written? Is it impossible?

Thanks, Claude