Calculating distance between a set of random points in hyperbolic space

Question

Given the metric of the Poincaré upper half-plane model

$$(ds)^2 = \frac{(dx)^2 + (dy)^2}{y^2}$$

and two known points $(x_1, y_1)$ and $(x_2, y_2)$ in the corresponding hyperbolic space $\mathbb{H} = \{ (x ,y) \mid y > 0; x, y \in \mathbb{R} \}$, can Mathematica give me the distance between those points without resorting to the distance formula

$$d((x_1, y_1), (x_2, y_2)) = \operatorname{arcosh} \left( 1 + \frac{ {(x_2 - x_1)}^2 + {(y_2 - y_1)}^2 }{ 2 y_1 y_2 } \right)$$

for example, via numerical integration?

The problem I have is, for some $t \in [0,1]$, I need a distance formula for a hyperbolic metric

$$\frac{e^{2ty}(dx)^2+(dy)^2}{y^2}$$

but this is difficult to obtain, since I need to solve the geodesic equations. All I really need is a set of distances between some Poisson points drawn from $\mathbb{H}$ with this new metric, in order to form a random graph.

Is there a way of doing this using e.g. a built in function(s) in Mathematica?

Answer 1

One can both solve the geodesic equations and integrate to find the length of the geodesic numerically. There are packages for working with tensors (notably xAct), but here's a quick-and-dirty hacked-together workflow.

It's far, far from good code—notably, I assume you keep the variable t free throughout the whole thing (as I use it to parameterize time), and use t0 for the constant in the exponential. As such, I strongly recommend starting with a fresh Global context by evaluating the ClearAll["Global`*"] at the beginning. Nevertheless, it works, and thought I'd share!

Also, I do something a bit weird: I use x and y as inert symbols for the indices themselves. I refer to the coordinate values (what might be typically denoted $x^1$ and $x^2$, or simply $x$ and $y$) as X[x] and X[y]. My definitions for the metric g should make it clear what I mean.

Also, this is predicated on whether I've transcribed the definition of the Christoffel symbol correctly! You'll want to check over my math for sure, and maybe run some checks to be sure the distance function looks like you think it should in special cases.

ClearAll["Global`*"]

(* Basic tensor/index management infrastructure: *)

indices = {x, y};

ipatt = Alternatives @@ indices;

toVect[expr_, i_] := expr /. {i -> #} & /@ indices

contract[t0_, t1_, i_] := toVect[t0, i] . toVect[t1, i]

(* The metric tensor: *)

g[x, x] = E^(2 t0 X[y])/X[y]^2;
g[x, y] = 0;
g[y, x] = 0;
g[y, y] = 1/X[y]^2;

(* Manually define the inverse: *)

invg[x, x] = 1/g[x, x];
invg[y, y] = 1/g[y, y];
invg[x, y] = 0;
invg[y, x] = 0;

(* Christoffel symbols: *)

\[CapitalGamma][g_][c : ipatt, a : ipatt, b : ipatt] := 
 With[{c0 = abstractIndex[Unique[]]}, contract[invg[c, c0],
   (1/2) (partial[g[c0, a], b] + partial[g[c0, b], a] - 
      partial[g[a, b], c0]), c0]]

(* Coordinate derivative that only evaluates when no indices
   are explicitly abstract: *)

partial[expr_?(FreeQ[_abstractIndex]), i : Except[_abstractIndex]] := 
 D[expr, X[i]]

(* Geodesic equation in abstract index c; \[Gamma][i][t] is to be
   the i'th component of our geodesic curve at time t: *)

eq = \[Gamma][c]''[
   t] == -contract[
    contract[\[CapitalGamma][g][c, a, b], \[Gamma][a]'[t], 
     a], \[Gamma][b]'[t], b];

(* All equations, and replace the coordinate values with their values
   on the curve; in a more careful world, g[a,b] would have been a function
   that takes in points on the manifold, like any good tensor field! *)

eqs = toVect[eq, c] /. X[s_] :> \[Gamma][s][t];

(* Construct the Dirichlet boundary conditions: *)

dirichletconditions[curve_, var_, pt0_, pt1_] := 
 Join[MapThread[curve[#1][0] == #2 &, {indices, pt0}], 
  MapThread[curve[#1][1] == #2 &, {indices, pt1}]]

(* Compute the geodesic, given a numeric parameter value for t0; 
   geodesic[t0] is of the form {fx[t], fy[t]} (using the actual
   symbol t) *)

geodesic[t00_?NumericQ][pt0_List, pt1_List] := 
 Replace[{a_} :> a]@
  NDSolveValue[
   Join[eqs // Simplify, 
     dirichletconditions[\[Gamma], t, pt0, pt1]] /. t0 -> t00, 
   Through[(\[Gamma] /@ indices)[t]], {t, 0, 1}]

(* Compute the distance between two points: *)

distance[t00_?NumericQ][pt0_, pt1_] := 
 With[{f = geodesic[t00][pt0, pt1]}, 
  With[{df = D[f, t]}, 
   NIntegrate[
    Evaluate[
     Sqrt[
      df . (toVect[toVect[g[a, b], b], a] /. 
         Join[{t0 -> t00}, MapThread[X[#1] :> #2 &, {indices, f}]]) . 
       df]], {t, 0, 1}]]]

(* Examples: *)

distance[0.5][{1, 2}, {2, 3}]

(* 1.452699322 *)

distance[0.5][{1, 1}, {1, 1}]

(* 3.33904718*10^-16 (nearly 0, good) *)

distance[0][{1, 1}, {3, 5}]

hdist[{x1_, y1_}, {x2_, y2_}] := ArcCosh[1 + ((x2 - x1)^2 + (y2 - y1)^2)/(2 y1 y2)]

hdist[{1, 1}, {3, 5}] // N

(* 1.762747161 and 1.762747174; good, we agree (up to numerical errors)
   with the hyperbolic metric at t0 == 0 *)

I'd like to clean up this code and make it more general instead of being so special-case. Maybe I'll come back to it!

Also, watch out for this: some combinations of points and parameter values will take an inordinate amount of time for apparently no reason. Consider the following three expressions:

distance[1.69][{1, 1}, {1.2, 1.2}]

distance[1.71][{1, 1}, {1.2, 1.2}]

distance[1.7][{1, 1}, {1.2, 1.2}]

The first two evaluate just fine after a couple seconds, but the last hangs. Weird! The perils of numerical methods.

In the meantime, is this sufficient? You mentioned you wanted to draw Poisson points from the upper half-plane; is that already done, and all you need to do now is compute distances between them? Or is that part of the challenge too?