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?

"Just tried" - how, exactly? You didn't show the code you ran. — J. M.'s missing motivation, Jan 24 '17 at 12:41
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] && 0 < r3 < 2 && 0 < a3 < 1, {r3, a3}, Reals ] — DM037, Jan 24 '17 at 12:43
Then yes, you've hit the limit of what NSolve[] can do; it was primarily designed for algebraic equations, and has very limited support for transcendental equations. corey's answer shows a different way to proceed. — J. M.'s missing motivation, Jan 24 '17 at 12:46
It's a pretty solution but not really feasible for me; I would like to implement this function in another code to find midpoints by the dozen. Would there be a way to do this without having to estimate solutions yourself, e.g by taking gradients? Thanks! — DM037, Jan 24 '17 at 12:50
In general, solving transcendental equations without having a general understanding of where the roots could be (that is, localization) is unproductive. If you can determine bounds (and the geometric source of the equations seems to suggest a route), then FindRoot[] can polish those bounds to the actual roots themselves. — J. M.'s missing motivation, Jan 24 '17 at 12:54

score 6 · Answer 1 · answered Jan 24 '17 at 12:26

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}

People interested in the ContourPlot[]-based approach will want to see this. — J. M.'s missing motivation, Jan 24 '17 at 12:29

Michael E2 · Answer 2 · 2017-01-25T02:34:03.270

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]}]

Mathematica graphics

score 3 · Answer 3 · answered Jan 24 '17 at 21:18

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.

NSolve difficulty to find midpoint hyperbolic geodesic

3 Answers3

Linked