NSolve on bounded domain

Question

I don't expect NSolve[] to be able to obtain all solutions in generalfor example. In principle there should be an algorithm that finds all solution to $f(x)=0$ for $x \in D$ in a bounded (compact) domain provided appropriate bounds are available that allow to rule out solutions in a sufficiently large ball around an evaluated point.

Can NSolve[] do this? Is such an algorithm implemented in Mathematica or another mathematical software package?

The problem at hand is the following

 NSolve[{1/2 Sqrt[3/2] y Cos[2 c] Sin[b]^2 + 
   x (1/2 Cos[2 c] Sin[2 b] Sin[a] + Cos[a] Sin[b] Sin[2 c]) == 0 &&
  0 <= c <= 2 π && 0 <= a <= 2 π && 
  0 <= b <= 2 π && -2 <= x <= 2 && -2 <= y <= 
  2, -Sqrt[(3/2)] y Cos[c] Sin[b]^2 Sin[c] + 
   x Sin[b] (Cos[a] Cos[2 c] - 2 Cos[b] Cos[c] Sin[a] Sin[c]) == 0 &&
  0 <= c <= 2 π && 0 <= a <= 2 π && 
  0 <= b <= 2 π && -2 <= x <= 2 && -2 <= y <= 2, 
  Sqrt[3/2] y Cos[b] Cos[c] Sin[b] + 
   x (Cos[2 b] Cos[c] Sin[a] + Cos[b] Cos[a] Sin[c]) == 0 && 
  0 <= c <= 2 π && 0 <= a <= 2 π && 
  0 <= b <= 2 π && -2 <= x <= 2 && -2 <= y <= 
  2, -Sqrt[(3/2)] y Cos[b] Sin[b] Sin[c] + 
   x (Cos[b] Cos[a] Cos[c] - Cos[2 b] Sin[a] Sin[c]) == 1 && 
  0 <= c <= 2 π && 0 <= a <= 2 π && 
  0 <= b <= 2 π && -2 <= x <= 2 && -2 <= y <= 2, 
  1/4 y (1 + 3 Cos[2 b]) - Sqrt[6] x Cos[b] Sin[b] Sin[a] == 0 && 
  0 <= c <= 2 π && 0 <= a <= 2 π && 
  0 <= b <= 2 π && -2 <= x <= 2 && -2 <= y <= 2}, {x, y, c, a, 
  b}, Reals]

Michael E2 · Accepted Answer · 2015-02-05T23:47:03.403

[Update notice: It seemed convenient to give a parametrization of the solutions. See mysols below.]

Probably NSolve should be able to do this, given enough time and memory, because the equations are reducible to a polynomial system. NSolve gives up after a while, though, and says it can't solve it. So we can manually convert it to a polynomial system by replacing sine and cosine by variables and adding a constraint that their squares sum to 1. There turns out to another issue in the corresponding large system with 8 variables, namely, precision. See below.

Clear[a, b, c, d, e, f, x, y];
eqns = {
 1/2 Sqrt[3/2] y Cos[2c] Sin[b]^2 + x (1/2 Cos[2c] Sin[2b] Sin[a] + Cos[a] Sin[b] Sin[2c]) == 0,
 -Sqrt[3/2] y Cos[c] Sin[b]^2 Sin[c] + x Sin[b] (Cos[a] Cos[2c] - 2 Cos[b] Cos[c] Sin[a] Sin[c]) == 0,
 Sqrt[3/2] y Cos[b] Cos[c] Sin[b] + x (Cos[2b] Cos[c] Sin[a] + Cos[b] Cos[a] Sin[c]) == 0,
 -Sqrt[3/2] y Cos[b] Sin[b] Sin[c] + x (Cos[b] Cos[a] Cos[c] - Cos[2 b] Sin[a] Sin[c]) == 1, 
 1/4 y (1 + 3 Cos[2 b]) - Sqrt[6] x Cos[b] Sin[b] Sin[a] == 0};
polys = eqns /.
 {Cos[2 v_] :> 2 v[1]^2 - 1, Cos[v_] :> v[1], Sin[2 v_] :> 2 v[1] v[2], Sin[v_] :> v[2]} /. 
 Thread[{a[1], b[1], c[1], a[2], b[2], c[2]} :> {a, b, c, d, e, f}];
cons = {
 a^2 + d^2 == 1, b^2 + e^2 == 1, c^2 + f^2 == 1,
 -1 <= c <= 1 && -1 <= a <= 1 && -1 <= b <= 1 && -1 <= d <= 1 && -1 <= e <= 1 && -1 <= f <= 1 &&
 -2 <= x <= 2 && -2 <= y <= 2};

With MachinePrecision no errors are reported but we get solutions that seemed suspicious.

NSolve[Join[polys, cons], {x, y, a, b, c, d, e, f}, Reals]
(*
  {{x -> 1., y -> 0, a -> 0, b -> 0, c -> 0, d -> -1., e -> -1., f -> -1.},
   ... (* 6 solutions with {d, e, f} -> ±1 omitted *)
   {x -> 1., y -> 0, a -> 0, b -> 0, c -> 0, d -> 1., e -> -1., f -> 1.}}
*)

Checking the equations with these solutions showed that there was some slight error, but enough to have some equations evaluate to False. This suggested checking NSolve again but with a higher WorkingPrecision. We get a different and interesting result:

sols = NSolve[Join[polys, cons], {x, y, a, b, c, d, e, f}, Reals,
  WorkingPrecision -> $MachinePrecision]

NSolve::infsolns: Infinite solution set has dimension at least 1. Returning intersection of solutions with (169443 a)/120544-(152085 b)/120544-(178587 c)/120544+(132357 d)/120544-(34779 e)/30136-(132147 f)/120544+(157165 x)/120544-(138533 y)/120544 == 1. >>

(*
  {{x -> 1.00000000000000, y -> 0, a -> 0.88420577536926, 
    b -> 1.00000000000000, c -> 0.88420577536926, 
    d -> 0.467097577390045, e -> 0, f -> -0.467097577390045},
    ...}  (* three solutions omitted *)
*)

Instead of a random hyperplane, might try to parametrize the solutions. The option Method -> {"UseSlicingHyperplanes" -> False} to NSolve will try to generate parametrized solution sets, but it ran for longer than I cared to wait. One could settle for a numeric parametrization:

mysols[a0_?NumericQ] := NSolve[Join[polys, cons, {a == a0}],
  {x, y, a, b, c, d, e, f}, Reals,
  WorkingPrecision -> 20]

Just a quick check of the solutions above as well as to show how to convert back to angles:

sols = Thread[{x, y, a, b, c} -> ({x, y, ArcTan[a, d], ArcTan[b, e], ArcTan[c, f]} /. #)] & /@ sols0;
eqns /. sols
(*
  {{True, True, True, True, True}, {True, True, True, True, True},
   {True, True, True, True, True}, {True, True, True, True, True},
   {True, True, True, True, True}, {True, True, True, True, True},
   {True, True, True, True, True}, {True, True, True, True, True}}
*)

Precision is still an issue. If we evaluate the equations at MachinePrecision using sols, some equations have enough of an error to evaluate to False:

eqns /. N@sols
(*
  {{False, False, False, True, False}, {True, True, False, True, True},
   {True, True, True, True, True}, {True, True, True, True, True},
   {True, True, False, True, True}, {False, False, False, True, False},
   {False, False, True, True, False}, {False, False, True, True, False}}
*)

NSolve on bounded domain

1 Answers1