# Calculating distance between a set of random points in hyperbolic space

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?

• Unfortunately the math is well over my head here. But I wonder: would it be possible for you to simplify the problem to its simplest form and give a concrete example, preferably in Mathematica code? Jun 1, 2022 at 12:00
• Yes, that can be done, will update the post.
– apg
Jun 1, 2022 at 12:01
• The distance is defined along the geodesic.. I think it is not possible to get the distance without knowing the geodesics. Jun 1, 2022 at 14:25
• But, if you knew the geodesics, I suppose one could then use numerical integration to get their length?
– apg
Jun 1, 2022 at 14:42
• That is right. Try it out. Jun 1, 2022 at 19:05

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] == #2 &, {indices, pt0}],
MapThread[curve[#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[{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?

• Yes thank you, this is exactly what I need. I am reading about solving the geodesic equations now. So, once you find the geodesics of the metric (are these circles for all $t$? Or only for the case $t=1$? I suppose the code will reveal this.), you can then integrate the arc length along them. For the Poisson points, I can just select $n=\text{Poisson}(\lambda)$ points from some region of the upper half plane?
– apg
Jun 2, 2022 at 10:06