3
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ Did you try to evaluate the parameters eIX,e1X,... and after that to use NSolve? $\endgroup$
    – user64494
    Commented Mar 16, 2021 at 15:54
  • $\begingroup$ No luck either with NSolve (yet). The expressions to the RHS of eIX... are already the outputs of FullSimplify given the assumptions assumps (if that's what you meant by evaluate) $\endgroup$
    – Anti Earth
    Commented Mar 16, 2021 at 21:00
  • $\begingroup$ You might do better in NSolve by removing the inequality restrictions and the domain restriction. This will require that the parameters e.g. eIX be given specific numeric values though. $\endgroup$ Commented Mar 16, 2021 at 22:57

1 Answer 1

1
$\begingroup$

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]

enter image description here

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

enter image description here

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.