2
$\begingroup$

Here I have a system of equations which takes $B$ matrix as input, and output the solutions. The code works very well for 4-rows matrices $B$ that I tried, however it is hard to give solutions in finite time when $B$ is equal to the following matrix.

Why this happens? Is it because the dimensions of the corresponding equations are too high to be dealt with by Solve (Reduce)? which level of complexity of equations that Mathematica Solve or Reduce functions cannot deal with? Is it possible to write a scalable code for this task?

B={{1,-1,0,0,1,0,0,0,0,-1},{0,1,-1,0,0,0,1,-1,0,0},{-1,0,0,0,0,1,-1,0,-1,0},{0,0,1,-1,0,0,0,0,1,1},{0,0,0,1,-1,-1,0,1,0,0}};
vec = {a, b, c, d, 0};
trigs = B . Sin[Transpose[B] . vec];
polys = TrigExpand /@ trigs /. {Sin[a] -> sa, Sin[b] -> sb, 
   Sin[c] -> sc, Sin[d] -> sd, Cos[a] -> ca, Cos[b] -> cb, 
   Cos[c] -> cc, Cos[d] -> cd};
allfuns = 
 Join[polys, {sa^2 + ca^2 - 1, sb^2 + cb^2 - 1, sc^2 + cc^2 - 1, 
   sd^2 + cd^2 - 1}];
gb = GroebnerBasis[allfuns, {sa, sb, sc, sd}, {ca, cb, cc, cd}];
sol = Reduce[gb == 0] /. {sa -> Sin[a], sb -> Sin[b], sc -> Sin[c], 
   sd -> Sin[d], ca -> Cos[a], cb -> Cos[b], cc -> Cos[c], 
   cd -> Cos[d]};
oldeqns = And @@ Thread[trigs == 0];
neweqns = FullSimplify[oldeqns && sol];
Reduce[neweqns && 0 <= a < 2 \[Pi] && 0 <= b < 2 \[Pi] && 
   0 <= c < 2 Pi && 0 <= d < 2 Pi, {a, b, c, d}, 
  Reals] // FullSimplify

Any suggestions? Thanks a lot!

Update

Here I give a try for the answer.

To look in detail of the Grobner basis produced in the code, I paste part of it due to its long length

{-1215 sd + 3888 sb^2 sd - 2592 sb^4 sd - 2304 sb^6 sd + 
  2304 sb^8 sd + 3888 sb sc sd - 5184 sb^3 sc sd - 6912 sb^5 sc sd + 
  9216 sb^7 sc sd + 3888 sc^2 sd - 7776 sb^2 sc^2 sd - 
  11520 sb^4 sc^2 sd + 16896 sb^6 sc^2 sd + 4096 sb^8 sc^2 sd - 
  5184 sb sc^3 sd - 11520 sb^3 sc^3 sd + 18432 sb^5 sc^3 sd + 
  16384 sb^7 sc^3 sd - 2592 sc^4 sd - 11520 sb^2 sc^4 sd + 
  19200 sb^4 sc^4 sd + 32768 sb^6 sc^4 sd - 6912 sb sc^5 sd + 
  18432 sb^3 sc^5 sd + 40960 sb^5 sc^5 sd - 2304 sc^6 sd + 
  16896 sb^2 sc^6 sd + 32768 sb^4 sc^6 sd + 9216 sb sc^7 sd + 
  16384 sb^3 sc^7 sd + 2304 sc^8 sd + 4096 sb^2 sc^8 sd + 
  3888 sb sd^2 - 5184 sb^3 sd^2 - 6912 sb^5 sd^2 + 9216 sb^7 sd^2 + 
  3888 sc sd^2 - 864 sb^2 sc sd^2 - 38016 sb^4 sc sd^2 + 
  32256 sb^6 sc sd^2 + 14336 sb^8 sc sd^2 - 864 sb sc^2 sd^2 - 
  64512 sb^3 sc^2 sd^2 + 50688 sb^5 sc^2 sd^2 + 
  73728 sb^7 sc^2 sd^2 - 5184 sc^3 sd^2 - 64512 sb^2 sc^3 sd^2 + 
  53760 sb^4 sc^3 sd^2 + 174080 sb^6 sc^3 sd^2 + 
  16384 sb^8 sc^3 sd^2 - 38016 sb sc^4 sd^2 + 53760 sb^3 sc^4 sd^2 + 
  247808 sb^5 sc^4 sd^2 + 65536 sb^7 sc^4 sd^2 - 6912 sc^5 sd^2 + 
  50688 sb^2 sc^5 sd^2 + 247808 sb^4 sc^5 sd^2 + 
  114688 sb^6 sc^5 sd^2 + 32256 sb sc^6 sd^2 + 
  174080 sb^3 sc^6 sd^2 + 114688 sb^5 sc^6 sd^2 + 9216 sc^7 sd^2 + 
  73728 sb^2 sc^7 sd^2 + 65536 sb^4 sc^7 sd^2 + 14336 sb sc^8 sd^2 + 
  16384 sb^3 sc^8 sd^2 + 3888 sd^3 - 7776 sb^2 sd^3 - 
  11520 sb^4 sd^3 + 16896 sb^6 sd^3 + 4096 sb^8 sd^3 ...

Nonetheless, it is already enough to see that the Grobner basis involves 8-degree terms, maybe this is the reason that Mathematica cannot find the solution. The Galois theory asserts that

a general polynomial of degree at least five cannot be solved by radicals

That is to say, there is no general formula to compute the roots of polynomials of higher degree.

$\endgroup$
2
  • $\begingroup$ Rational roots of polynomials over the rationals of arbitrary degree are easily found. $\endgroup$
    – user64494
    Feb 18, 2022 at 8:09
  • $\begingroup$ @user64494 Does this has contradiction to Galois theory ? $\endgroup$
    – M.K
    Feb 18, 2022 at 8:20

3 Answers 3

2
$\begingroup$

How to get at least a few solutions.

Reduce first for Sin of variables and then Solve for variables.

(B = {{1, -1, 0, 0, 1, 0, 0, 0, 0, -1}, {0, 1, -1, 0, 0, 0, 1, -1, 0, 
 0}, {-1, 0, 0, 0, 0, 1, -1, 0, -1, 0}, {0, 0, 1, -1, 0, 0, 0, 0, 
 1, 1}, {0, 0, 0, 1, -1, -1, 0, 1, 0, 0}}) // MatrixForm

vec = {a, b, c, d, 0};
(trigs = B.Sin[Transpose[B].vec]) // TableForm

trigs0 = Thread[TrigExpand[trigs] == 0];

condsinlist = {-1 <= Sin[a] <= 1, -1 <= Sin[b] <= 1, -1 <= Sin[c] <= 
1, -1 <= Sin[d] <= 1};

varslist = {a, b, c, d};

condvars = -Pi <= a <= Pi && -Pi <= b <= Pi && -Pi <= c <= 
Pi && -Pi <= d <= Pi;

redall0 = Reduce[Join[trigs0, condsinlist], Reals];

le0b = LogicalExpand[#] & /@ redall0;

TraditionalForm[
  le0b //. Or -> 
Composition[(Column[#, Right, Background -> {{White, LightGray}}, 
   Frame -> All] &), List]]

(tabb = Table[{TimeConstrained[
  Solve[le0b[[i]] && condvars, varslist, Reals, Method -> Reduce],
   3], "i =", i}, {i, 1, Length[le0b]}]) // TableForm

enter image description here

trigs0 /. DeleteCases[Flatten[tabb[[All, 1]], 1], $Aborted]

(*   True   *)

Examine the $Aborted equations (aborted after 3 seconds) further to may be get more solutions.

$\endgroup$
4
  • $\begingroup$ I got the error Join::heads: Heads Equal and List at positions 1 and 2 are expected to be the same. Reduce::naqs: Join[B.Sin[Transpose[B].{a,b,c,d,0}]==0,{-1<=Sin[a]<=1,-1<=Sin[b]<=1,-1<=Sin[c]<=1,-1<=Sin[d]<=1}] is not a quantified system of equations and inequalities. $\endgroup$
    – M.K
    Feb 16, 2022 at 17:45
  • $\begingroup$ Seems you have not defined B. Start with a fresh Kenel (Quit ) and copy/paste my code. Works well in version 8.0 $\endgroup$
    – Akku14
    Feb 16, 2022 at 19:23
  • $\begingroup$ Thanks, it works. Does the code find all solutions or part of them? $\endgroup$
    – M.K
    Feb 17, 2022 at 12:17
  • $\begingroup$ As i wrote "Examine the $Aborted equations (aborted after 3 seconds) further to may be get more solutions." May be you can get some more analytical/exact solutions. But i think, some can only be found numericaly. $\endgroup$
    – Akku14
    Feb 17, 2022 at 12:42
1
$\begingroup$

If I correctly understand your "and what kind of complexity of equations that Mathematica cannot solve?", here is an example of simple trig equation on which Reduce, Solve, and FindInstance stuck:

FindInstance[Sin[x]*Sin[1755*x]*Sin[2011*x] == 1 && x > -Pi && x <= Pi, x, Reals]

The above equation is solved by hand as follows. The equation Sin[x]*Sin[1755*x]*Sin[2011*x] == 1 over the reals implies Sin[x]==1 or Sin[x]==-1. It remains to test x==Pi/2 and x==-Pi/2. Only the former satisfies the equation.

$\endgroup$
11
  • $\begingroup$ Can the down vote be explained? What is wrong in my answer? Deep regard. $\endgroup$
    – user64494
    Feb 14, 2022 at 19:01
  • $\begingroup$ BTW, the equation under consideration has complex roots too, e.g. $- 3.00144974730400- 0.000885659243711501 \,\mathrm{I}$. $\endgroup$
    – user64494
    Feb 14, 2022 at 19:08
  • $\begingroup$ Thanks a lot! Yes, this example gets stuck as well. Could you explain more about why this happens? (I just noticed that the down vote is canceled, maybe someone click it by mistake: ) ) $\endgroup$
    – M.K
    Feb 14, 2022 at 20:04
  • 2
    $\begingroup$ @user64494 , Reduce can do your equation Reduce[{s1 s2 s3 == 1, s1 == Sin[x], s2 == Sin[1755 x ], s3 == Sin[2011 x], -Pi < x <= Pi, -1 <= s1 <= 1, -1 <= s2 <= 1, -1 <= s3 <= 1}, {x, s1, s2, s3}, Reals] yields x == \[Pi]/2 && s1 == 1 && s2 == -1 && s3 == -1 (version 8.0) $\endgroup$
    – Akku14
    Feb 14, 2022 at 20:48
  • 1
    $\begingroup$ @user64494 c = Catalan; Solve[{s1 s2 s3 == 1, s1 == Sin[x^c], s2 == Sin[1755 x^c], s3 == Sin[2011 x^c], -Pi < x <= Pi, -1 <= s1 <= 1, -1 <= s2 <= 1, -1 <= s3 <= 1}, {x, s1, s2, s3}, Reals, Method -> Reduce] $\endgroup$
    – cvgmt
    Feb 15, 2022 at 3:11
1
$\begingroup$

A lot of problems cannot be solved using exact techniques and require numeric solution. Search this site for examples of FindRoot and NSolve. For the example provided by user64494,

Clear["Global`*"]

eqn = Sin[x]*Sin[1755*x]*Sin[2011*x] == 1;

The approximate numeric solution is

soln = NSolveValues[eqn && -Pi < x <= Pi, x, Reals]

(* {1.5708} *)

The exact value is

sol = RootApproximant[soln[[1]]/Pi]*Pi

(* π/2 *)

Verifying,

eqn /. x -> sol

(* True *)

EDIT: When NSolve fails, FindRoot with a good initial estimate may work.

eqn2 = 
  Sin[x^Catalan]*Sin[1755*x^Catalan]*Sin[2011*x^Catalan] == 1;

Plot[Evaluate[List @@ eqn2], {x, 0, Pi}]

enter image description here

arg2 = NArgMax[{eqn2[[1]], 1 < x < 2}, x]

(* 1.63724 *)

sol2 = FindRoot[eqn2, {x, arg2}, WorkingPrecision -> 15]

(* {x -> 1.63724137606215} *)

eqn2 /. sol2

(* True *)

Another example,

EDIT 3:

eqn3 = Sin[Surd[x, 3]]*Sin[1755*Surd[x, 3]]*Sin[2011*Surd[x, 3]] == 1;

eqn3a = eqn3 /. Surd[x, 3] :> y

(* Sin[y] Sin[1755 y] Sin[2011 y] == 1 *)

This is the same equation as eqn

sol3a = NSolve[{eqn3a, -10 <= y <= 10}, y]

(* {{y -> -4.71239}, {y -> 1.5708}, {y -> 7.85398}} *)

Then the solution to eqn3 on the interval -10 <= x <= 10 is

sol3 = sol3a /. (y -> z_) :> (x -> (RootApproximant[z/Pi]*Pi)^3)

(* {{x -> -((27 π^3)/8)}, {x -> π^3/8}, {x -> (125 π^3)/8}} *)

Verifying,

eqn3 /. sol3

(* {True, True, True} *)
$\endgroup$
11
  • $\begingroup$ Got it! So if I want to have a scalable code that works for different $B$ with different dimensions, I should choose numerical approaches (numerical algebraic geometry methods). Am I right? $\endgroup$
    – M.K
    Feb 14, 2022 at 21:18
  • $\begingroup$ Think of Sin[x^Catalan]*Sin[1755*x^Catalan]*Sin[2011*x^Catalan]==1;, where NSolveValues incorrectly returns {}. $\endgroup$
    – user64494
    Feb 14, 2022 at 21:18
  • $\begingroup$ Also think of Sin[Surd[x, 3]]*Sin[1755*Surd[x, 3]]*Sin[2011*Surd[x, 3]] == 1; where NSolveValues is running without any response for a long time. $\endgroup$
    – user64494
    Feb 14, 2022 at 23:32
  • $\begingroup$ @user64494 - When NSolve fails use FindRoot. See edit. $\endgroup$
    – Bob Hanlon
    Feb 15, 2022 at 0:13
  • $\begingroup$ Yes, FindRoot[eqn2, {x, arg2}, WorkingPrecision -> 15] numerically solves eqn2 = Sin[x^Catalan]*Sin[1755*x^Catalan]*Sin[2011*x^Catalan] == 1;. However, how to identify it as a symbolic result? $\endgroup$
    – user64494
    Feb 15, 2022 at 4:35

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.