Solving the Euler-Lagrange equation with NDSolveProcessEquations

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[u][t] ==
-(Subscript[Γ, ϕ, ϕϕ][u[t], v[t]]*Derivative[u][t]^2 +
2*Subscript[Γ, ϕ, ϕμ][u[t], v[t]]*Derivative[u][t]*
Derivative[v][t] +
Subscript[Γ, ϕ, μμ][u[t], v[t]]*Derivative[v][t]^2),
Derivative[v][t] ==
-(Subscript[Γ, μ, ϕϕ][u[t], v[t]]*Derivative[u][t]^2 +
2*Subscript[Γ, μ, ϕμ][u[t], v[t]]*Derivative[u][t]*
Derivative[v][t] +
Subscript[Γ, μ, μμ][u[t], v[t]]*Derivative[v][t]^2),
{u, v} == a, {u, v} == 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[
NDSolveProcessEquations[
{Derivative[u][t] ==
-(Subscript[Γ, ϕ, ϕϕ][u[t], v[t]]*Derivative[u][t]^2 +
2*Subscript[Γ, ϕ, ϕμ][u[t], v[t]]*Derivative[u][t]*Derivative[v][t] +
Subscript[Γ, ϕ, μμ][u[t], v[t]]*Derivative[v][t]^2),
Derivative[v][t] ==
-(Subscript[Γ, μ, ϕϕ][u[t], v[t]]*Derivative[u][t]^2 +
2*Subscript[Γ, μ, ϕμ][u[t], v[t]]*Derivative[u][t]*Derivative[v][t] +
Subscript[Γ, μ, μμ][u[t], v[t]]*Derivative[v][t]^2),
{u, v} == {A[], A[]},
{u, v} == {B[], B[]}},
{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 = NDSolveReinitialize[EuLa, {u, v} == A, {u, v} == B];

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

NDSolveReinitialize::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[NDSolveProcessEquations[{Derivative[u][t] == -(Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Phi]][u[t], v[t]]*Derivative[u][t]^2 +
2*Subscript[\[CapitalGamma], \[Phi], \[Phi]\[Mu]][u[t], v[t]]*Derivative[u][t]*Derivative[v][t] + Subscript[\[CapitalGamma], \[Phi], \[Mu]\[Mu]][u[t], v[t]]*Derivative[v][t]^2),
Derivative[v][t] == -(Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Phi]][u[t], v[t]]*Derivative[u][t]^2 + 2*Subscript[\[CapitalGamma], \[Mu], \[Phi]\[Mu]][u[t], v[t]]*Derivative[u][t]*
Derivative[v][t] + Subscript[\[CapitalGamma], \[Mu], \[Mu]\[Mu]][u[t], v[t]]*Derivative[v][t]^2), {u, v} == {A[], A[]},
{u, v} == {B[], B[]}}, {u, v}, t]]
NDSolveIterate[EuLa, 0]
NDSolveIterate[EuLa, 1]
initial = NDSolveProcessSolutions[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[u], Derivative[v]} /. initial

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

new1 = First[NDSolveReinitialize[EuLa, {{u, v} == A, {Derivative[u], Derivative[v]} == slope}]]
NDSolveIterate[new1, 0]
NDSolveIterate[new1, 1]
sol1 = NDSolveProcessSolutions[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[NDSolveReinitialize[EuLa, {{u, v} == A, {Derivative[u], Derivative[v]} == B}]]
NDSolveIterate[new2, 0]
NDSolveIterate[new2, 1]
sol2 = NDSolveProcessSolutions[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[NDSolveReinitialize[EuLa, {{u, v} == A, {u, v} == B}]]
NDSolveIterate[new3, 0]
NDSolveIterate[new3, 1]
sol3 = NDSolveProcessSolutions[new3]
NDSolveReinitialize::ndcinit: Initial conditions should be specified at a single point.

Am I right?

Claude

• Please show us the code text of the example, not pictures of part of it. If you don't know how to, check this post: meta.mathematica.stackexchange.com/q/1584/1871 – xzczd Feb 8 '17 at 10:57
• A first error that I can see is that you are assigning your equations to a variable named EuLa, but when you want to solve them you call to Eula. – dpravos Feb 8 '17 at 11:01
• Value of Chif is missing, and the part defining solGlob causes infy warning and fails, at least in v9.0.1. Please check the sample carefully before posting it here. Then, have you tried ParametricNDSolve? – xzczd Feb 9 '17 at 10:19
• Chif is the precision (MachinePrecision, or another reasonable number). I use version 11.0.1, and I can plot solGlob without any problem (what is infy?). I also tried ParametricNDSolve ; the result was the same (excepted for pairs of very distant points, where the geodesic is hard to obtain), but the computation was much lengthy. Thus, I chose NDSolve, even if ParametricNDSolve seemed to work better in few pathological situations. – C. Manté Feb 9 '17 at 10:52
• Just tested on Cloud, OK, seems that NDSolve is improved in v11. In v9.0.1 I only got Power::infy warning, and even if I manually set up "Shooting" method, I need to give a very good initial guess and manually set up "ImplicitSolver" as Method -> {"Shooting", "StartingInitialConditions" -> With[{b = 0}, {{\[Phi][b], \[Mu][b]} == Rationalize[A, 0], {Derivative[\[Phi]][b], Derivative[\[Mu]][b]} == Rationalize[{0.22659514943257, 4.5093449250158}, 0]}], "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}}. Then, back to your question… – xzczd Feb 9 '17 at 11:35