Solving the Euler-Lagrange equation with NDSolve`ProcessEquations

Question

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 much time by using NDSolveProcessEquations coupled with NDSolveReinitialize.

Here are the Christoffel symbols:

Subscript[Γ, ϕ, ϕϕ][ϕ_, μ_] := 
  -((μ*(μ + 2*ϕ) + ϕ^2*(μ + ϕ) 
     ((-(ϕ/(μ + ϕ))^ϕ)*(μ + (μ + ϕ)*Log[ϕ/(μ + ϕ)])*PolyGamma[1, ϕ] - 
       (μ + ϕ)*(-1 + (ϕ/(μ + ϕ))^ϕ)*PolyGamma[2, ϕ]))/ 
          (2*ϕ*(μ + ϕ)*(μ + ϕ*(μ + ϕ)*(-1 + (ϕ/(μ + ϕ))^ϕ)*PolyGamma[1, ϕ]))); 
Subscript[Γ, ϕ, ϕμ][ϕ_, μ_] := 
  -((ϕ*(-1 + ϕ*(ϕ/(μ + ϕ))^ϕ*(μ + ϕ)*PolyGamma[1, ϕ])) /
    (2*(μ + ϕ)*(μ + ϕ*(μ + ϕ)*(-1 + (ϕ/(μ + ϕ))^ϕ)*PolyGamma[1, ϕ]))); 
Subscript[Γ, ϕ, μμ][ϕ_, μ_] := 
  ϕ/(2*(μ + ϕ)*(μ + ϕ*(μ + ϕ)*(-1 + (ϕ/(μ + ϕ))^ϕ)*PolyGamma[1, ϕ])); 
Subscript[Γ, μ, ϕϕ][ϕ_, μ_] := 
  (1/2)*μ*(1/(μ*ϕ + ϕ^2) - (ϕ/(μ + ϕ))^ϕ*PolyGamma[1, ϕ]); 
Subscript[Γ, μ, ϕμ][ϕ_, μ_] := μ/(2*ϕ*(μ + ϕ)); 
Subscript[Γ, μ, μμ][ϕ_, μ_] := -((2*μ + ϕ)/(2*μ*(μ + ϕ)));

Here is the Euler-Lagrange equation, depending on the Christoffel symbols (Chif is the Precision):

GeodE[{u_, v_}, t_, a_, b_] := 
  SetPrecision[
    {Derivative[2][u][t] == 
       -(Subscript[Γ, ϕ, ϕϕ][u[t], v[t]]*Derivative[1][u][t]^2 + 
          2*Subscript[Γ, ϕ, ϕμ][u[t], v[t]]*Derivative[1][u][t]*
            Derivative[1][v][t] + 
          Subscript[Γ, ϕ, μμ][u[t], v[t]]*Derivative[1][v][t]^2), 
     Derivative[2][v][t] == 
       -(Subscript[Γ, μ, ϕϕ][u[t], v[t]]*Derivative[1][u][t]^2 + 
         2*Subscript[Γ, μ, ϕμ][u[t], v[t]]*Derivative[1][u][t]*
           Derivative[1][v][t] + 
         Subscript[Γ, μ, μμ][u[t], v[t]]*Derivative[1][v][t]^2), 
     {u[0], v[0]} == a, {u[1], v[1]} == b}, 
    Chif];

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}};

The geodesic can be easily found from

{A, B} = {Easy4a, Easy4b};
solGlob =
  TimeConstrained[
    NDSolve[GeodE[{ϕ, μ}, t, A, B], {ϕ, μ}, {t, 0,1}, WorkingPrecision -> Chif], 
    240, "Pas fini"];

To determine geodesics between other pair of points, I defined alternatively (as in the documentation of "Components and Data Structures")

 EuLa = 
   First[
     NDSolve`ProcessEquations[
       {Derivative[2][u][t] == 
          -(Subscript[Γ, ϕ, ϕϕ][u[t], v[t]]*Derivative[1][u][t]^2 + 
          2*Subscript[Γ, ϕ, ϕμ][u[t], v[t]]*Derivative[1][u][t]*Derivative[1][v][t] + 
          Subscript[Γ, ϕ, μμ][u[t], v[t]]*Derivative[1][v][t]^2), 
        Derivative[2][v][t] == 
          -(Subscript[Γ, μ, ϕϕ][u[t], v[t]]*Derivative[1][u][t]^2 +  
          2*Subscript[Γ, μ, ϕμ][u[t], v[t]]*Derivative[1][u][t]*Derivative[1][v][t] + 
          Subscript[Γ, μ, μμ][u[t], v[t]]*Derivative[1][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:

Centre = {0.7767`, 11.207780999999999`};
TestC = {0.7767`, 87.26795000000001`};
B = TestC;
A = Centre;
solBis = NDSolve`Reinitialize[EuLa, {u[0], v[0]} == A, {u[1], v[1]} == B];

I corrected "Eula", but there is still an error:

NDSolve`Reinitialize::ndcinit: Initial conditions should be specified at a single point.

I'm using NDSolve (and the "components and data structure" possibilities) for the first time, and I can't see where the problem is. Is the formuation of EuLa ill-written?

Is it impossible?

Yes, it seems. I have to solve some "à la Dirichlet" problem:

{A, B} = {Easy4a, Easy4b};
EuLa = First[NDSolve`ProcessEquations[{Derivative[2][u][t] == -(Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Phi]][u[t], v[t]]*Derivative[1][u][t]^2 + 
    2*Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Mu]][u[t], v[t]]*Derivative[1][u][t]*Derivative[1][v][t] + Subscript[\[CapitalGamma], \[Phi], \[Mu]\[Mu]][u[t], v[t]]*Derivative[1][v][t]^2), 
 Derivative[2][v][t] == -(Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Phi]][u[t], v[t]]*Derivative[1][u][t]^2 + 2*Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Mu]][u[t], v[t]]*Derivative[1][u][t]*
     Derivative[1][v][t] + Subscript[\[CapitalGamma], \[Mu], \[Mu]\[Mu]][u[t], v[t]]*Derivative[1][v][t]^2), {u[0], v[0]} == {A[[1]], A[[2]]}, 
 {u[1], v[1]} == {B[[1]], B[[2]]}}, {u, v}, t]]
NDSolve`Iterate[EuLa, 0]
NDSolve`Iterate[EuLa, 1]
initial = NDSolve`ProcessSolutions[EuLa];
Show[Graphics[{PointSize[0.03], Red, Point[A]}], Graphics[{PointSize[0.07], GrayLevel[0.6], Point[B]}],
ParametricPlot[Evaluate[{u[t], v[t]} /. initial], {t, 0, 1}, PlotRange -> All], Frame -> True, AspectRatio -> 1]

slope = {Derivative[1][u][0], Derivative[1][v][0]} /. initial

Thanks to "slope", we obtain exactly the same solution by solving with NDSolve`Reinitialize the equivalent à la Cauchy problem:

new1 = First[NDSolve`Reinitialize[EuLa, {{u[0], v[0]} == A, {Derivative[1][u][0], Derivative[1][v][0]} == slope}]]
NDSolve`Iterate[new1, 0]
NDSolve`Iterate[new1, 1]
sol1 = NDSolve`ProcessSolutions[new1]
Show[Graphics[{PointSize[0.03], Red, Point[A]}], Graphics[{PointSize[0.07], GrayLevel[0.6], Point[B]}],` ParametricPlot[Evaluate[{u[t], v[t]} /. sol1], {t, 0, 1}, PlotRange -> All], Frame -> True, AspectRatio -> 1]`

The Cauchy boundary condition can be changed, giving rise to another valid solution:

{A, B} = {TestC, Centre}
new2 = First[NDSolve`Reinitialize[EuLa, {{u[0], v[0]} == A, {Derivative[1][u][0], Derivative[1][v][0]} == B}]]
NDSolve`Iterate[new2, 0]
NDSolve`Iterate[new2, 1]
sol2 = NDSolve`ProcessSolutions[new2]
Show[Graphics[{PointSize[0.03], Red, Point[A]}],
ParametricPlot[Evaluate[{u[t], v[t]} /. sol2], {t, 0, 1},
PlotRange -> All], Frame -> True, AspectRatio -> 1]

It could be interesting to compute the exponential map, but it seems impossible to solve my Dirichlet problem under these conditions:

new3 = First[NDSolve`Reinitialize[EuLa, {{u[0], v[0]} == A, {u[1], v[1]} == B}]]
NDSolve`Iterate[new3, 0]
NDSolve`Iterate[new3, 1]
sol3 = NDSolve`ProcessSolutions[new3]
NDSolve`Reinitialize::ndcinit: Initial conditions should be specified at a single point.

Am I right?

Claude