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 NDSolve
Reinitialize...
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
Attachments: