4
$\begingroup$

I have a system of algebraic equations. For example: $$ \left\{ \begin{aligned} x^2 y + 2 x y + 2 &= 0,\\ y^3 + 2 x + 1 &= 0. \end{aligned} \right. $$ (In my problem the system contains 16 equations of degrees 1–3).

I need to know the number of complex roots (the system in question has a finite number of roots).

I believe, the number of roots can be determined using GroebnerBasis. But what is the most elegant way to do it?

$\endgroup$
  • $\begingroup$ Your question would get more attention if you provide Mathematica code instead of $\LaTeX$ :) $\endgroup$ – Öskå Jul 1 '14 at 10:50
  • $\begingroup$ {x^2 y + 2 x y + 2 == 0, y^3 + 2 x + 1 == 0} But the question is not about solving this particular system. $\endgroup$ – bcp Jul 1 '14 at 11:43
  • $\begingroup$ @bcp Getting help is sometimes about making it easy for potential answerers to verify their ideas on your illustrative example. $\endgroup$ – Michael E2 Jul 2 '14 at 0:49
5
$\begingroup$

The simplest method, that will usually be correct, is to just use NSolve.

polys = {x^2 y + 2 x y + 2, y^3 + 2 x + 1};
Length[NSolve[polys]]

(* Out[421]= 7 *)

Below is from code I have used for solution counting and related (it is related to code inside NSolve, to answer the question of "Why not just use NSolve?"). Some is also used in this MSE response. The underlying method is described in Gröbner Bases: A Computational Approach to Commutative Algebra, by T. Becker and V. Weispfenning (with H. Kredel).

normalSet[elist_List, maxlist_, nvars_Integer] /; 
  Length[elist] == nvars := 
 Module[{jj, indices, iterators}, indices = Array[jj, {nvars}];
  iterators = Table[{indices[[j]], 0, maxlist[[j]] - 1}, {j, nvars}];
  Flatten[Table[indices, Evaluate[Apply[Sequence, iterators]]], 
   nvars - 1]]

getRoofTuples[elist_, height_] := 
 Module[{len = Length[elist], expvec, expvecs, tt, rest}, 
  expvecs = tt[];
  Do[expvec = elist[[jj]];
   If[expvec[[height]] =!= 
      0 && (height === 1 || Max[Take[expvec, height - 1]] === 0), 
    expvecs = tt[expvec, expvecs]], {jj, len}];
  expvecs = Apply[List, Flatten[expvecs, Infinity, tt]];
  rest = Complement[elist, expvecs];
  expvecs = Map[Drop[#, height - 1] &, expvecs];
  {expvecs, rest}]

isUnder[t1_, t2_] := Apply[And, Thread[t1 <= t2]]

hitRoof[tup_, roofs_] := 
 Module[{jj, len = Length[roofs]}, 
  For[jj = 1, jj <= len, jj++, 
   If[isUnder[roofs[[jj]], tup], Return[True]]];
  False]

normalSet[elist_List, maxes_, nvars_Integer] := 
 Module[{tuples, ntups, tuple, ntup, mixedexpons, tt}, 
  mixedexpons = 
   Select[elist, (Max[Apply[Sequence, #]] != Apply[Plus, #]) &];
  tuples = Table[{j}, {j, 0, maxes[[nvars]] - 1}];
  Do[{roofs, mixedexpons} = getRoofTuples[mixedexpons, jj];
   ntups = tt[];
   Do[tuple = tuples[[kk]];
    ntups = tt[Prepend[tuple, 0], ntups];
    Do[ntup = Prepend[tuple, mm];
     If[hitRoof[ntup, roofs], Break[], ntups = tt[ntup, ntups]], {mm, 
      maxes[[jj]] - 1}], {kk, Length[tuples]}];
   tuples = Flatten[ntups, Infinity, tt], {jj, nvars - 1, 1, -1}];
  Sort[Apply[List, tuples]]]

solutionSetSize[polys_, vars_] := 
 Catch[Module[{lv = Length[vars], leadvecs, dtl, maxes, nset, 
    pureterms},
   gb = GroebnerBasis[polys, vars, 
     MonomialOrder -> DegreeReverseLexicographic];
   dtl = First[
     GroebnerBasis`DistributedTermsList[gb, vars, 
      MonomialOrder -> DegreeReverseLexicographic]];
   leadvecs = Map[#[[1, 1]] &, dtl];
   pureterms = Select[leadvecs, Count[#, a_ /; a =!= 0] === 1 &];
   If[Length[pureterms] < lv, Throw[Infinity]];
   maxes = Apply[Plus, pureterms];
   nset = normalSet[leadvecs, maxes, lv];
   Throw[Length[nset]]
   ]]

Back to the example:

solutionSetSize[polys, {x, y}]

(* Out[423]= 7 *)
$\endgroup$
1
$\begingroup$

The documentation is very instructive.

For your particular example:

gb = GroebnerBasis[{x^2 y + 2 x y + 2, y^3 + 2 x + 1}, {x, y}]

Finding roots:

sol = N@Solve[gb[[1]] == 0, y]

=>{{y -> -1.29819}, {y -> -0.871821 - 1.29622 I}, {y -> -0.871821 + 1.29622 I}, {y -> 0.274176 - 1.20102 I}, {y -> 0.274176 + 1.20102 I}, {y -> 1.24674 - 0.331078 I}, {y -> 1.24674 + 0.331078 I}}

sol2 = Solve[gb[[2]] == 0, x] /. sol

=> {{{x -> 0.593927}}, {{x -> -2.36593 + 0.388879 I}}, {{x -> -2.36593 - 0.388879 I}}, {{x -> 0.0829231 - 0.730783 I}}, {{x -> 0.0829231 + 0.730783 I}}, {{x -> -1.26396 + 0.753779 I}}, {{x -> -1.26396 - 0.753779 I}}}

Just for fun showing real root:

ContourPlot[
 Evaluate@Thread[{x^2 y + 2 x y + 2, y^3 + 2 x + 1} == 0], {x, -10, 
  10}, {y, -10, 10}, 
 Epilog -> {PointSize[0.02], Red, 
   Point[{{First[x /. First[sol2]], First[y /. sol]}}]}]

enter image description here

$\endgroup$
  • $\begingroup$ So the question is whether I can determine the number of roots without actually finding all of them. (e.g. if I had only one single-variable polynomial equation, I would take its degree as the answer) $\endgroup$ – bcp Jul 1 '14 at 13:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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