Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

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

2 Answers 2

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 Gro"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 *)
share|improve this answer

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

share|improve this answer
    
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) –  bcp Jul 1 at 13:48

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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