Solve system of trigonometric polynomials

Question

I am trying to solve a real system of trigonometric polynomials, which contain 7 symbols {b[1],b[2],b[3],b[4],a[2],a[3],a[4]} between 7 equations:

equs = {
 eIX == 2 (b[1] b[2] Cos[a[2]] + b[3] b[4] Cos[a[3] - a[4]]),
 eXI == 2 (b[1] b[3] Cos[a[3]] + b[2] b[4] Cos[a[2] - a[4]]),
 eXX == 2 (b[2] b[3] Cos[a[2] - a[3]] + b[1] b[4] Cos[a[4]]),
 eIZ == b[1]^2 - b[2]^2 + b[3]^2 - b[4]^2,
 eZI == b[1]^2 + b[2]^2 - b[3]^2 - b[4]^2,
 eZZ == b[1]^2 - b[2]^2 - b[3]^2 + b[4]^2,
 1 == b[1]^2 + b[2]^2 + b[3]^2 + b[4]^2
}

All symbols are real and finite, and furthermore bounded:

assumps = {
 0 <= b[1] <= 1,
 0 <= b[2] <= 1,
 0 <= b[3] <= 1,
 0 <= b[4] <= 1,
 0 <= a[2] <= 2 Pi,
 0 <= a[3] <= 2 Pi,
 0 <= a[4] <= 2 Pi,
 -1 <= eIX <= 1,
 -1 <= eXI <= 1,
 -1 <= eXX <= 1,
 -1 <= eIZ <= 1,
 -1 <= eZI <= 1,
 -1 <= eZZ <= 1
}

Mathematica 11.2 on my MacBook Pro seems unable to efficiently solve this system for the 7 aforementioned symbols, via commands like:

Solve[
  Join[equs, assumps], 
  {b[1], b[2], b[3], b[4], a[2], a[3], a[4]}, 
  Reals]

Reduce[
  Join[equs, assumps], 
  {b[1], b[2], b[3], b[4], a[2], a[3], a[4]}, 
  Reals]

Worrying about the transcendental functions, its easy to recast this system into a strictly polynomial one of 10 variables {b[1],b[2],b[3],b[4],c[2],c[3],c[4],s[2],s[3],s[4]} (the new c and s replacing cos and sin) and 10 equations:

poly = Join[
  equs  /. Cos[a_ - b_] :> Cos[a] Cos[b] + Sin[a] Sin[b]
        /. {Cos[a[n_]] -> c[n], Sin[a[n_]] -> s[n]},
  Table[1 == c[n]^2 + s[n]^2, {n, 2, 4}]]

polyassumps = Join[
  DeleteCases[assumps, _ <= a[_] <= _], 
  Flatten @ Table[{-1 <= c[n] <= 1, -1 <= s[n] <= 1}, {n, 2, 4}]]

The system is now 10 degree-4 polynomial equations:

poly = {
 eIX == 2 (b[1] b[2] c[2] + b[3] b[4] (c[3] c[4] + s[3] s[4])),
 eXI == 2 (b[1] b[3] c[3] + b[2] b[4] (c[2] c[4] + s[2] s[4])),
 eXX == 2 (b[1] b[4] c[4] + b[2] b[3] (c[2] c[3] + s[2] s[3])),
 eIZ == b[1]^2 - b[2]^2 + b[3]^2 - b[4]^2,
 eZI == b[1]^2 + b[2]^2 - b[3]^2 - b[4]^2,
 eZZ == b[1]^2 - b[2]^2 - b[3]^2 + b[4]^2,
 1 == b[1]^2 + b[2]^2 + b[3]^2 + b[4]^2,
 1 == c[2]^2 + s[2]^2,
 1 == c[3]^2 + s[3]^2,
 1 == c[4]^2 + s[4]^2
}

polyassumps = {
 0 <= b[1] <= 1, 
 0 <= b[2] <= 1, 
 0 <= b[3] <= 1, 
 0 <= b[4] <= 1, 
 -1 <= eIX <= 1, 
 -1 <= eXI <= 1, 
 -1 <= eXX <= 1, 
 -1 <= eIZ <= 1, 
 -1 <= eZI <= 1, 
 -1 <= eZZ <= 1, 
 -1 <= c[2] <= 1, 
 -1 <= s[2] <= 1, 
 -1 <= c[3] <= 1, 
 -1 <= s[3] <= 1, 
 -1 <= c[4] <= 1, 
 -1 <= s[4] <= 1
}

Still, Solve and Reduce struggle:

Solve[
  Join[poly,polyassumps],
  {b[1], b[2], b[3], b[4], c[2], c[3], c[4], s[2], s[3], s[4]},
  Reals]

Is there a more efficient way to solve simultaneous equations of these kind in Mathematica?

Answer 1

The question requests that equs be solved symbolically to obtain {b[1], b[2], b[3], b[4], a[2], a[3], a[4]} in terms of {eIX, eXI, eXX, eIZ, eZI, eZZ}, subject to the constraints assumps. Obtaining solutions for {b[1], b[2], b[3], b[4]} can be obtained using Solve, producing sixteen sets of b[_], only one of which is entirely nonnegative.

Solve[equs[[4 ;;]], {b[1], b[2], b[3], b[4]}] // Last
(* {b[1] -> 1/2 Sqrt[1 + eIZ + eZI + eZZ], 
    b[2] -> 1/2 Sqrt[1 - eIZ + eZI - eZZ], 
    b[3] -> 1/2 Sqrt[1 + eIZ - eZI - eZZ], 
    b[4] -> 1/2 Sqrt[1 - eIZ - eZI + eZZ]} *)

Visibly, b[_] are real, only if

1 + eIZ + eZI + eZZ >= 0 && 1 - eIZ + eZI - eZZ >= 0 && 
    1 + eIZ - eZI - eZZ >= 0 && 1 - eIZ - eZI + eZZ >= 0;

(A nearly identical constraint can be obtained from Solve[equs[[4 ;;]], {b[1], b[2], b[3], b[4]}, Reals] and then extracting the condition included in the solution, but doing so is more cumbersome than the approach above.) It is instructive to plot this constraint:

RegionPlot3D[%, {eIZ, -1, 1}, {eZI, -1, 1}, {eZZ, -1, 1}, PlotPoints -> 100, 
    ViewPoint -> {-2, -2, 1.4}, AxesLabel -> {eIZ, eZI, eZZ}, ImageSize -> Large, 
    LabelStyle -> {15, Bold, Black}, Mesh -> None]

Thus, the constraint derived here satisfies assumps[[11 ;; 13]] but is more restrictive. For some values of {eIZ, eZI, eZZ}, there are no nonnegative solutions. For instance,

NSolve[Join[equs[[4 ;;]] /. {eIZ -> -.9, eZI -> -.9, eZZ -> -.9}, assumps[[;; 4]]], 
    {b[1], b[2], b[3], b[4]}, Reals]
(* {} *)

More generally, only 1/3 of {eIZ, eZI, eZZ} chosen randomly from the cube in the plot result in solutions.

Volume@ImplicitRegion[1 + eIZ + eZI + eZZ >= 0 && 1 - eIZ + eZI - eZZ >= 0 && 
    1 + eIZ - eZI - eZZ >= 0 && 1 - eIZ - eZI + eZZ >= 0, {eIZ, eZI, eZZ}]
(* 8/3 *)

Next, consider equs[[ ;; 3]]. To take advantage of Mathematica's powerful capabilities for solving polynomials, employ the Weierstrass Substitution.

equw = Simplify[TrigExpand[equs] /. {Sin[a[i_]] -> 2 t[i]/(1 + t[i]^2), 
    Cos[a[i_]] -> (1 - t[i]^2)/(1 + t[i]^2)}];
equw[[ ;; 3]]
(* {eIX + (2 b[1] b[2] (-1 + t[2]^2))/(1 + t[2]^2) == (2 b[3] b[4] (1 + 4 t[3] t[4] - 
    t[4]^2 + t[3]^2 (-1 + t[4]^2)))/((1 + t[3]^2) (1 + t[4]^2)), 
    eXI + (2 b[1] b[3] (-1 + t[3]^2))/(1 + t[3]^2) == (2 b[2] b[4] (1 + 4 t[2] t[4] - 
    t[4]^2 + t[2]^2 (-1 + t[4]^2)))/((1 + t[2]^2) (1 + t[4]^2)), 
    eXX + (2 b[1] b[4] (-1 + t[4]^2))/(1 + t[4]^2) == (2 b[2] b[3] (1 + 4 t[2] t[3] - 
    t[3]^2 + t[2]^2 (-1 + t[3]^2)))/((1 + t[2]^2) (1 + t[3]^2)) *)

(Of course, equw[[4 ;; ]] = equs[[4 ;; ]].) The transformed equations are very useful for obtaining numerical solutions. For instance,

{eIX -> -0.15, eXI -> -0.72, eXX -> 0.35, eIZ -> -0.55, eZI -> -0.21, eZZ -> -0.23};
NSolve[Join[equw /. %, assumps[[;; 4]]], {b[1], b[2], b[3], b[4], t[2], t[3], t[4]}, 
    Reals] /. Rule[t[i_], n_] -> Rule[a[i], 2 ArcTan[n]]
(* {{b[1] -> 0.05, b[2] -> 0.626498, b[3] -> 0.471699, b[4] -> 0.618466, 
       a[2] -> -1.06941, a[3] -> -1.89781, a[4] -> 2.50076}, 
    {b[1] -> 0.05, b[2] -> 0.626498, b[3] -> 0.471699, b[4] -> 0.618466, 
       a[2] -> 1.06941, a[3] -> 1.89781, a[4] -> -2.50076}, 
    {b[1] -> 0.05, b[2] -> 0.626498, b[3] -> 0.471699, b[4] -> 0.618466, 
       a[2] -> 2.39568, a[3] -> 1.34392, a[4] -> -0.406055}, 
    {b[1] -> 0.05, b[2] -> 0.626498, b[3] -> 0.471699, b[4] -> 0.618466, 
       a[2] -> -2.39568, a[3] -> -1.34392, a[4] -> 0.406055}} *)

I have chosen to display a[_] in the range {-Pi, Pi} rather than the equivalent {0, 2 Pi} in order to highlight the symmetries between the first and second, and between the third and fourth solutions. I have solved equw for a few other sets of parameters, each yielding four solutions with the same symmetries. As above, it is instructive to plot the region in {eIX, eXI, eXX} space for which solutions exist for the values of b[_] just computed. (Warning: Plot requires several minutes to complete.)

equs[[;; 3, 2]] /. %[[1, ;; 4]];
Region[ParametricRegion[%, {{a[2], -Pi, Pi}, {a[3], -Pi, Pi}, {a[4], -Pi, Pi}}], 
    Boxed -> True, Axes -> True, ImageSize -> Large], AxesLabel -> {eIX, eXI, eXX}, 
    PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}, LabelStyle -> {15, Bold, Black},
    ViewPoint -> {-2, -2, 1.4}]

Not only is the region shown small but it is hollow. For example,

Element[{0, 0, 0}, ParametricRegion[%%, {{a[2], -Pi, Pi}, {a[3], -Pi, Pi}, 
    {a[4], -Pi, Pi}}]]
(* False *)

I have attempted to obtain a symbolic solution for equw[[ ;; 3]] but so far without success. Reduce, the obvious choice, ran for 12 hours without yielding an answer.