How to use Mathematica's GroebnerBasis to automate solving system of polynomials

Question

Currently, Only way I know to use GroebnerBasis to solve system of polynomials require lots of manual steps. I'd like to know how to automate this process and if Mathematica has commands to help with this.

For an example, given these 4 polynomials in x,y,z,r

vars = {x,y,z,r}
eq1  = 3 x^2+2 y z-2 x r==0;
eq2  = 2 x z - 2 y r==0;
eq3  = 2 x y-2 z-2 z r==0;
eq4  = x^2+y^2+z^2-1==0;
opt={MonomialOrder->Lexicographic};

One can solve these using Solve

Solve[{eq1,eq2,eq3,eq4},vars]

And get these solutions

 {{x->-1,y->0,z->0,r->-(3/2)},
  {x->-(2/3),y->-(1/3),z->-(2/3),r->-(4/3)},
   etc....
 }

To obtain the same solutions using GroebnerBasis, I start by finding basis

 g=GroebnerBasis[{eq1,eq2,eq3,eq4},vars,opt]

which gives

36 r-225 r^2-493 r^3-116 r^4+212 r^5+96 r^6
-4 z+25 r z+53 r^2 z+24 r^3 z
-189 r-4338 r^2-2076 r^3+1928 r^4+960 r^5+1105 z^2
5 r y+5 r^2 y+2 z+r z-r^2 z
4347 r+6291 r^2+12 r^3-2796 r^4-864 r^5+2210 y z+2210 r y z
-9945+10782 r+55171 r^2+25232 r^3-22556 r^4-13344 r^5+9945 y^2
13747 r+23859 r^2+4788 r^3-10604 r^4-5280 r^5+3315 x+3315 y z

There are 7 basis polynomials. Here comes the manual steps.

I look to see which polynomial above has only one variable in it. In this case, one is lucky. There is one. The first one depends only on r. So use that to solve for r

 Solve[g[[1]]==0,r]

Now I look to see which other polynomial in the basis has r and only one other variable to solve for.

This will be the second one in the list. So use it to solve for z using the first r value found above. (this process is very much like Gaussian elimination, the back substitution phase)

 Solve[(g[[2]]/.r->-3/2)==0,z]

Now we know r and z. So now will look for a basis polynomial which has r and z in it and only one more unknown to solve for. This will be the 4th one. Hence

 Solve[(g[[4]]/.{r->-3/2,z->0})==0,y]

Now we know r,z,y. Then look for polynomial which has these and one more unknown, which is x. This will be the last basis. Hence

 Solve[(g[[7]]/.{r->-3/2,z->0,y->0})==0,x]

So the one of the solution from above is

 {x->-1,y->0,z->0,r->-(3/2)}

This is one of the solutions found by Solve. To find the others, I repeat the above process, but now starting with the second r solution found in the first step of the process.

But this manual process is time consuming.

Does Mathematica have build in functions to automate this process once the GroebnerBasis are found?

Answer 1

As noted in comments, Solve does this. An efficient "by hand" approach is to solve for the first equation, plug the solutions into the next, solve, join corresponding solutions. This can be done in two steps. First the example set-up.

vars = {x, y, z, r};
eq1 = 3 x^2 + 2 y z - 2 x r == 0;
eq2 = 2 x z - 2 y r == 0;
eq3 = 2 x y - 2 z - 2 z r == 0;
eq4 = x^2 + y^2 + z^2 - 1 == 0;
opt = {MonomialOrder -> Lexicographic};
g = GroebnerBasis[{eq1, eq2, eq3, eq4}, vars, opt]

(* Out[108]= {36 r - 225 r^2 - 493 r^3 - 116 r^4 + 212 r^5 + 
  96 r^6, -4 z + 25 r z + 53 r^2 z + 24 r^3 z, -189 r - 4338 r^2 - 
  2076 r^3 + 1928 r^4 + 960 r^5 + 1105 z^2, 
 5 r y + 5 r^2 y + 2 z + r z - r^2 z, 
 4347 r + 6291 r^2 + 12 r^3 - 2796 r^4 - 864 r^5 + 2210 y z + 
  2210 r y z, -9945 + 10782 r + 55171 r^2 + 25232 r^3 - 22556 r^4 - 
  13344 r^5 + 9945 y^2, 
 13747 r + 23859 r^2 + 4788 r^3 - 10604 r^4 - 5280 r^5 + 3315 x + 
  3315 y z} *)

Here is a helper function. Given a partial solution as a list of rules, and a polynomial, plug values into the equation and solve it. I am skipping details about whether/when one should expect the new polynomial to be univariate: a sufficient condition is that the system have a finite solution set (this is a consequence of the Gianni/Kalkbrenner theorems from I believe 1987).

solvenext[lastsol_, poly_] := 
 Map[Join[lastsol, #] &, Solve[(poly /. lastsol) == 0]]

Now we deploy it, iterating from first polynomial in the Groebner basis to last. Fold is handy for this (okay, I first did it with a Do loop, then decided not to act my age and went with a functional programming construct instead).

Fold[
 Apply[Join, Table[solvenext[prevsol, #2], {prevsol, #1}]] &, {{}}, g]

(* Out[149]= {{r -> -(3/2), z -> 0, y -> 0, x -> -1}, {r -> -(4/3), 
  z -> -(2/3), y -> -(1/3), x -> -(2/3)}, {r -> -(4/3), z -> 2/3, 
  y -> 1/3, x -> -(2/3)}, {r -> -1, z -> -1, y -> 0, 
  x -> 0}, {r -> -1, z -> -1, y -> 0, x -> 0}, {r -> -1, z -> 1, 
  y -> 0, x -> 0}, {r -> -1, z -> 1, y -> 0, x -> 0}, {r -> 0, z -> 0,
   y -> -1, x -> 0}, {r -> 0, z -> 0, y -> 1, x -> 0}, {r -> 1/8, 
  z -> -(Sqrt[(11/2)]/8), y -> (3 Sqrt[11/2])/8, 
  x -> -(3/8)}, {r -> 1/8, z -> Sqrt[11/2]/8, 
  y -> -((3 Sqrt[11/2])/8), x -> -(3/8)}, {r -> 3/2, z -> 0, y -> 0, 
  x -> 1}} *)