NSolve difficulty to find midpoint hyperbolic geodesic

Question

I want to find the midpoint of a geodesic in hyperbolic coordinates. I have already understood that, for general coordinates {r1,a2}, {r2,a2}, to find the midpoint {r3,a3} is difficult. So, to make things easier, I assume r1=1, a1=0, r2=2 and a2=1, so the system of equations becomes:

NSolve[Cosh[1]*Cosh[r3] - Sinh[1]*Sinh[r3]*Cos[-a3] == 
   0.5*(Cosh[1]*Cosh[2] - Sinh[1]*Sinh[2]*Cos[-1] ) && 
  0.5*(Cosh[1]*Cosh[2] - Sinh[1]*Sinh[2]*Cos[-1] ) == 
   Cosh[2]*Cosh[r3] - Sinh[2]*Sinh[r3]*Cos[1 - a3] , {r3, a3}, Reals]

I get

"NSolve[Cosh[1] Cosh[r3] - Cos[a3] Sinh[1] Sinh[r3] == 1.75122 && 1.75122 == Cosh[2] Cosh[r3] - Cos[1 - a3] Sinh[2] Sinh[r3], {r3, a3}, Reals]"

in return, instead of a solution. I think it is surprising that this problem is so hard to solve, even numerically. Or am I doing something wrong?

Answer 1

FindRoot can find specific roots given proper starting values:

eq = {Cosh[1]*Cosh[r3] - Sinh[1]*Sinh[r3]*Cos[-a3] == 
    0.5*(Cosh[1]*Cosh[2] - Sinh[1]*Sinh[2]*Cos[-1]), 
   0.5*(Cosh[1]*Cosh[2] - Sinh[1]*Sinh[2]*Cos[-1]) == 
    Cosh[2]*Cosh[r3] - Sinh[2]*Sinh[r3]*Cos[1 - a3]};

ContourPlot[Evaluate[eq], {r3, -Pi, Pi}, {a3, -Pi, Pi}, 
 FrameLabel -> {"r3", "a3"}]

FindRoot[eq, {r3, -1.8}, {a3, -2.5}]

{r3 -> -1.65282, a3 -> -2.53314}

FindRoot[eq, {r3, -1}, {a3, -2}]

{r3 -> -0.887494, a3 -> -1.95154}

Answer 2

Solve the first equation for r3 and plug into the second equation:

eqn = Cosh[1]*Cosh[r3] - Sinh[1]*Sinh[r3]*Cos[-a3] == 
    1/2*(Cosh[1]*Cosh[2] - Sinh[1]*Sinh[2]*Cos[-1]) && 
   1/2*(Cosh[1]*Cosh[2] - Sinh[1]*Sinh[2]*Cos[-1]) == 
    Cosh[2]*Cosh[r3] - Sinh[2]*Sinh[r3]*Cos[1 - a3];

solr3 = Reduce[First@eqn, {r3}, Reals];

$ineqPat = Less | Greater | LessEqual | GreaterEqual | Inequality;
tmp = solr3 /. e : ($ineqPat[___] | Cos[a3] == _) :> 
   Simplify[e, {a3, r3} ∈ Reals && -1 <= Cos[a3] <= 1];
sol1 = Solve[tmp, r3];

tmp = Join[#, {"a3" -> a3}] /. NSolve[Last@eqn && -Pi < a3 < Pi /. #, a3] & /@ sol1;
sol = Flatten[tmp /. "a3" -> a3, 1]
(*
  {{r3 -> 1.65282, a3 -> 0.608453}, {r3 -> 0.887494, a3 -> 1.19005},
   {r3 -> -1.65282, a3 -> -2.53314}, {r3 -> -0.887494, a3 -> -1.95154}}
*)

The general solution may be then represented by

(r3, a3 + 2 Pi C[1]} /. sol

where C[1] ∈ Integers.

Check:

ContourPlot[List @@ eqn // Evaluate, {a3, -Pi, Pi}, {r3, -2, 2},
 Epilog -> {Red, PointSize@Medium, Point[{a3, r3} /. sol]}]

Answer 3

Might make sense to solve for the trig and hyperbolic values because that can be done as a polynomial system. Then use those and perhaps myltiple inverses to recover the actual values.

To do this, expand the hyper/trig stuff and substitute with new variables. Also add polynomials to account for the basic identities.

eqns = {Cosh[1]*Cosh[r3] - Sinh[1]*Sinh[r3]*Cos[-a3] == 
    0.5*(Cosh[1]*Cosh[2] - Sinh[1]*Sinh[2]*Cos[-1]), 
   0.5*(Cosh[1]*Cosh[2] - Sinh[1]*Sinh[2]*Cos[-1]) == 
    Cosh[2]*Cosh[r3] - Sinh[2]*Sinh[r3]*Cos[1 - a3]};
exprs = TrigExpand[Apply[Subtract, eqns, {1}]];
subs = {Cosh[a_] :> ch[a], Cos[a_] :> c[a], Sinh[a_] :> sh[a], 
   Sin[a_] :> s[a]};
newpolys = {ch[r3]^2 - sh[r3]^2 - 1, c[a3]^2 + s[a3]^2 - 1};
exprs2 = Join[N[exprs] /. subs, newpolys]

(* Out[21]= {-1.75122291618 + 1.54308063482 ch[r3] - 
  1.17520119364 c[a3] sh[r3], 
 1.75122291618 - 3.76219569108 ch[r3] + 1.95960104142 c[a3] sh[r3] + 
  3.05189779915 s[a3] sh[r3], -1 + ch[r3]^2 - sh[r3]^2, -1 + c[a3]^2 +
   s[a3]^2} *)

solns = NSolve[exprs2]

(* Out[52]= {{c[a3] -> 0.820533531224, ch[r3] -> 2.70660617606, 
  s[a3] -> 0.571598394099, 
  sh[r3] -> 2.51509780969}, {c[a3] -> -0.820533531224, 
  ch[r3] -> 2.70660617606, s[a3] -> -0.571598394099, 
  sh[r3] -> -2.51509780969}, {c[a3] -> -0.371612463763, 
  ch[r3] -> 1.42036028817, s[a3] -> -0.928387945191, 
  sh[r3] -> -1.00867405449}, {c[a3] -> 0.371612463763, 
  ch[r3] -> 1.42036028817, s[a3] -> 0.928387945191, 
  sh[r3] -> 1.00867405449}} *)

You can pick up some solutions as below.

vars = 
  Cases[Variables[exprs2], _c | _ch] /. {c[a_] :> ArcCos[c[a]], 
    s[a_] :> ArcSin[s[a]], ch[a_] :> ArcCosh[ch[a]], 
    sh[a_] :> ArcSinh[sh[a]]};
realsols = vars /. solns

(* Out[54]= {{0.60845253015, 1.65282378268}, {2.53314012344, 
  1.65282378268}, {1.95154158767, 0.887493788519}, {1.19005106592, 
  0.887493788519}} *)

To get some others consider using different branches of the inverse trigs and hyperbolics.