Could Macaulay2 can solve polynomial equations?

123 views
Skip to first unread message

Arrigo Coen

unread,
Oct 27, 2020, 3:19:32 PM10/27/20
to Macaulay2
Hello,

Could Macaulay2 can solve polynomial equations equivalent to Mathematica's `Solve` function?

Something like
Solve[a x + b y + c z == d && 3 x + 4 y + 5 z == 6 && 7 x + 7 y + 8 z == 9, {x, y, z}]
to get
{{x -> (3 (b - 2 c + d))/(3 a - 11 b + 7 c), 
  y -> -((3 a - 15 c + 11 d)/(3 a - 11 b + 7 c)), 
  z -> -((-6 a + 15 b - 7 d)/(3 a - 11 b + 7 c))}}

Thanks!

jche...@gmail.com

unread,
Oct 27, 2020, 5:38:21 PM10/27/20
to Macaulay2
This is certainly possible, but (understandably) requires slightly different setup for M2 syntax. One way to do this (for a linear system, as in your example) is as follows:

K = frac(QQ[a,b,c,d])
R = K[x,y,z]
m = matrix{{a*x + b*y + c*z, 3*x + 4*y + 5*z, 7*x + 7*y + 8*z}}
B = matrix{{d}, {6}, {9}}
A = transpose(m // vars R)
(inverse A)*B

(Of course, solving nonlinear systems requires different methods, but M2 can certainly handle many non-linear systems as well.)

Justin

Arrigo Coen

unread,
Oct 27, 2020, 6:26:30 PM10/27/20
to Macaulay2

Thank Justin for your comment. In fact, I have a nonlinear system (nonlinear in the sense that there are products between the variables, e.g., z0*z11). An example of my equations is:

```
Solve[{
1-(2/3)*z11  ==  (1-gamR)*(1-(2/3)*zR0*zR10)+gamR*(1-(2/3)*zR10) , 
(1/3)*z11  ==  (1-gamR)*(1/3)*zR0*zR10+gamR*(1/3)*zR10 , 
(1/3)*z11  ==  (1-gamR)*(1/3)*zR0*zR10+gamR*(1/3)*zR10 , 
(1-gam)*(1-(2/3)*z11)+gam*(1-(2/3)*z0*z11)  ==  (1-gamR)*(1-(2/3)*zR1*zR0*zR10)+gamR*(1-(2/3)*zR10) , 
(1-gam)*(1/3)*z11+gam*(1/3)*z0*z11  ==  (1-gamR)*(1/3)*zR1*zR0*zR10+gamR*(1/3)*zR10 , 
(1-gam)*(1/3)*z11+gam*(1/3)*z0*z11  ==  (1-gamR)*(1/3)*zR1*zR0*zR10+gamR*(1/3)*zR10 , 
(1-gam)*(1-(2/3)*z1)+gam*(1/3)*z0  ==  (1-gamR)*(1-(2/3)*zR1)+gamR*(1/3)*zR0 ,
(1-gam)*(1/3)*z1+gam*(1-(2/3)*z0)  ==  (1-gamR)*(1/3)*zR1+gamR*(1-(2/3)*zR0) ,
(1-gam)*(1/3)*z1+gam*(1/3)*z0  ==  (1-gamR)*(1/3)*zR1+gamR*(1/3)*zR0 , 
(1-gam)*(1-(2/3)*z1*z11)+gam*(1-(2/3)*z11)  ==  1-(2/3)*zR0*zR10 , 
(1-gam)*(1/3)*z1*z11+gam*(1/3)*z11  ==  (1/3)*zR0*zR10 , 
(1-gam)*(1/3)*z1*z11+gam*(1/3)*z11  ==  (1/3)*zR0*zR10 }, 
{ z0, z1, z11, zR0, zR1, zR10, gam, gamR }]
```


And what I am looking for are the regions like this
```
{z1 -> zR1, z11 -> 0, zR0 -> 0, gam -> 0, gamR -> 0}
```
where all the equations are fulfilled. These regions tell me the restrictions that fulfill the equations. 

Do you know how can I found these regions? 

Best regards, 
A

jche...@gmail.com

unread,
Oct 27, 2020, 7:59:55 PM10/27/20
to Macaulay2
Yes, M2 can find these "regions" for you (and indeed this is one of the core functions of Macaulay2). The solution set to a polynomial system is called a(n algebraic) variety, and the "regions" you're looking for are called the components of the variety. Finding these regions goes by the name of primary decomposition in commutative algebra:

R = QQ[z0,z1,z11,zR0,zR1,zR10,gam,gamR]
I = ideal(1-(2/3)*z11-(1-gamR)*(1-(2/3)*zR0*zR10)+gamR*(1-(2/3)*zR10),(1/3)*z11-(1-gamR)*(1/3)*zR0*zR10+gamR*(1/3)*zR10,(1/3)*z11-(1-gamR)*(1/3)*zR0*zR10+gamR*(1/3)*zR10,(1-gam)*(1-(2/3)*z11)+gam*(1-(2/3)*z0*z11)-(1-gamR)*(1-(2/3)*zR1*zR0*zR10)+gamR*(1-(2/3)*zR10),(1-gam)*(1/3)*z11+gam*(1/3)*z0*z11-(1-gamR)*(1/3)*zR1*zR0*zR10+gamR*(1/3)*zR10,(1-gam)*(1/3)*z11+gam*(1/3)*z0*z11-(1-gamR)*(1/3)*zR1*zR0*zR10+gamR*(1/3)*zR10,(1-gam)*(1-(2/3)*z1)+gam*(1/3)*z0-(1-gamR)*(1-(2/3)*zR1)+gamR*(1/3)*zR0,(1-gam)*(1/3)*z1+gam*(1-(2/3)*z0)-(1-gamR)*(1/3)*zR1+gamR*(1-(2/3)*zR0),(1-gam)*(1/3)*z1+gam*(1/3)*z0-(1-gamR)*(1/3)*zR1+gamR*(1/3)*zR0,(1-gam)*(1-(2/3)*z1*z11)+gam*(1-(2/3)*z11)-1-(2/3)*zR0*zR10,(1-gam)*(1/3)*z1*z11+gam*(1/3)*z11-(1/3)*zR0*zR10,(1-gam)*(1/3)*z1*z11+gam*(1/3)*z11-(1/3)*zR0*zR10)
primaryDecomposition I

The output of the last command above shows you that there are 4 components, which are in fact linear spaces of dimension 3 (including the one you already mentioned). 

Justin

jche...@gmail.com

unread,
Oct 27, 2020, 8:10:58 PM10/27/20
to Macaulay2
I just noticed that I was too hasty in claiming linearity of the components above. There are in fact 2 linear (= degree 1) components and 2 quadric (= degree 2) components. The full output of primaryDecomposition I is:

{ideal(gamR,gam,zR10,z11,z1-zR1), ideal(gamR,gam,zR0,z11,z1-zR1),
     ideal(gamR,zR10,z11,z0-1,z1*gam-z1+zR1-gam),
     ideal(gamR,zR0,z11,z0-1,z1*gam-z1+zR1-gam)}

Justin

Arrigo Coen

unread,
Oct 28, 2020, 5:51:45 PM10/28/20
to Macaulay2
Thank you a lot Justin!!!
I just run the Macaulay2 code that you post, but besides the "primaryDescomposition of I", I get also another output

```
              2                2           2            2             1                1           1            1     1                1           1            1       2                    2               2             2          2            2             1                    1               1             1          1            1     1                    1               1             1          1            1     1         2         1           2           2     2                    2         1         2           1           1     1                  1         1         1           1           1     1     2             2         2           2           1             1         1           1           1             1         1           1
o2 = ideal (- -zR0*zR10*gamR + -zR0*zR10 - -zR10*gamR - -z11 + 2gamR, -zR0*zR10*gamR - -zR0*zR10 + -zR10*gamR + -z11, -zR0*zR10*gamR - -zR0*zR10 + -zR10*gamR + -z11, - -zR0*zR1*zR10*gamR + -zR0*zR1*zR10 - -z0*z11*gam + -z11*gam - -zR10*gamR - -z11 + 2gamR, -zR0*zR1*zR10*gamR - -zR0*zR1*zR10 + -z0*z11*gam - -z11*gam + -zR10*gamR + -z11, -zR0*zR1*zR10*gamR - -zR0*zR1*zR10 + -z0*z11*gam - -z11*gam + -zR10*gamR + -z11, -z0*gam + -z1*gam + -zR0*gamR - -zR1*gamR - -z1 + -zR1 - gam + gamR, - -z0*gam - -z1*gam - -zR0*gamR + -zR1*gamR + -z1 - -zR1 + gam + gamR, -z0*gam - -z1*gam + -zR0*gamR + -zR1*gamR + -z1 - -zR1, -z1*z11*gam - -z1*z11 - -zR0*zR10 - -z11*gam, - -z1*z11*gam + -z1*z11 - -zR0*zR10 + -z11*gam, - -z1*z11*gam + -z1*z11 - -zR0*zR10 + -z11*gam)
              3                3           3            3             3                3           3            3     3                3           3            3       3                    3               3             3          3            3             3                    3               3             3          3            3     3                    3               3             3          3            3     3         3         3           3           3     3                    3         3         3           3           3     3                  3         3         3           3           3     3     3             3         3           3           3             3         3           3           3             3         3           3

o2 : Ideal of R

i3 : primaryDecomposition I

o3 = {ideal (gamR, zR10, z11, z0 - 1, z1*gam - z1 + zR1 - gam), ideal (gamR, gam, zR10, z11, z1 - zR1), ideal (gamR, zR0, z11, z0 - 1, z1*gam - z1 + zR1 - gam), ideal (gamR, gam, zR0, z11, z1 - zR1)}
```

What is the meaning of o2? 

When I run the equivalent in Mathematica I get those four components of the variety, but also others. My guess is that o2 are these others components, am I right?

jche...@gmail.com

unread,
Oct 28, 2020, 7:58:26 PM10/28/20
to Macaulay2
o2 just refers to "output 2", i.e. the output of the second line. This is the ideal defined by the polynomial system, whose zero set is the variety you want. The output of the primary decomposition (the 3rd line) shows that this ideal is an intersection of 4 prime ideals, so there should only be 4 components (if Mathematica gives other components, I'd be very interested to see what they are.) For instance, the component you gave

{z1 -> zR1, z11 -> 0, zR0 -> 0, gam -> 0, gamR -> 0} 

is the last element in the list o3, corresponding to the (vanishing set of the) ideal (gamR, gam, zR0, z11, z1 - zR1).

Justin

Arrigo Coen

unread,
Oct 29, 2020, 2:43:30 PM10/29/20
to Macaulay2
Mathematica gave me this output:
```
{{z0 -> (gam - gamR + gamR*zR0)/gam, 
  z1 -> (gam - gamR - zR1 + gamR*zR1)/(-1 + gam), 
     z11 -> 0, zR10 -> 0}, {z1 -> zR1, z11 -> 0, zR0 -> 0, gam -> 0,   gamR -> 0}, 
   {z1 -> 1, z11 -> zR0*zR10, zR1 -> 1, gam -> 0, gamR -> 0}, 
   {z1 -> zR1, z11 -> 0, zR10 -> 0, gam -> 0, gamR -> 0}, 
   {z1 -> 1, z11 -> zR10, zR0 -> 1, gam -> 0, gamR -> 1}, 
   {z1 -> 1, z11 -> zR10, zR0 -> 1, zR1 -> 1, gam -> 0}, 
   {z1 -> gamR + zR1 - gamR*zR1, z11 -> 0, zR0 -> 1, zR10 -> 0,   gam -> 0}, 
   {z0 -> 1, z11 -> zR0*zR10, zR1 -> 1, gam -> 1, gamR -> 0}, 
   {z0 -> 1, z11 -> zR10, zR0 -> 1, gam -> 1, gamR -> 1}, 
   {z0 -> zR0, z11 -> 0, zR10 -> 0, gam -> 1, gamR -> 1}, 
   {z0 -> 1, z11 -> zR10, zR0 -> 1, zR1 -> 1, gam -> 1}, 
   {z0 -> 1 - gamR + gamR*zR0, z11 -> 0, zR1 -> 1, zR10 -> 0,   gam -> 1}, 
   {z0 -> 1, z1 -> (gam - zR1)/(-1 + gam), z11 -> 0, zR0 -> 0,   gamR -> 0}, 
   {z0 -> 1, z1 -> 1, z11 -> zR0*zR10, zR1 -> 1, gamR -> 0}, 
   {z0 -> 1, z1 -> (gam - zR1)/(-1 + gam), z11 -> 0, zR10 -> 0,  gamR -> 0}, 
   {z0 -> 1, z1 -> 1, z11 -> zR10, zR0 -> 1, gamR -> 1}, 
   {z0 -> 1, z1 -> 1, z11 -> zR10, zR0 -> 1, zR1 -> 1}, 
   {z0 -> 1, z1 -> (gam - gamR - zR1 + gamR*zR1)/(-1 + gam),   z11 -> 0, zR0 -> 1,      zR10 ->    0}, 
   {z0 -> (gam - gamR - 2*gam*gamR + 2*gamR^2 + gam*gamR^2)/(gam*(-1 +gamR)^2),z1 -> (-gam + gamR + gam*gamR)/((-1 + gam)*(-1 + gamR)),  z11 -> -((gamR*zR10)/(-1 + gamR)), zR0 -> gamR^2/(-1 + gamR)^2,   zR1 -> gamR^2/(-1 + gamR)^2}, 
   {z0 -> (gam - gamR - 2*gam*gamR + 2*gamR^2 +  gam*gamR^2)/(gam*(-1 + gamR)^2),  z1 -> (gam - gamR - zR1 + gamR*zR1)/(-1 + gam), z11 -> 0, zR0 -> gamR^2/(-1 + gamR)^2, zR10 -> 0}, 
   {z0 -> 1, z11 -> 0,   zR0 -> 0, zR1 -> 1,  gam -> 1, gamR -> 0}, {z0 -> 1, z11 -> 0, zR1 -> 1, zR10 -> 0, gam -> 1, gamR -> 0}, 
   {z0 -> 1, z11 -> 0, zR0 -> 1, zR1 -> 1, zR10 -> 0, gam -> 1}, 
   {z0 -> (1 - 3*gamR + 3*gamR^2)/(-1 + gamR)^2, z11 -> 0,  zR0 -> gamR^2/(-1 + gamR)^2, zR1 -> 1, zR10 -> 0, gam -> 1}}
```
If I understand you correctly, this Mathematica output is giving intersections of the prime ideals, is that right?

Arrigo

Arrigo Coen

unread,
Oct 29, 2020, 2:45:11 PM10/29/20
to Macaulay2
Also,  I forgot to mention that, when I run the Mathematica code I get the message
```

```

jche...@gmail.com

unread,
Oct 29, 2020, 5:29:32 PM10/29/20
to Macaulay2
Thanks for the reply - it's helpful to see what Mathematica gives as output (although it seems the Mathematica message didn't actually get posted). Unfortunately, although I don't know how Mathematica is producing this output, I can say that some of it is simply wrong. E.g. the 3rd component that Mathematica gives, namely

{z1 -> 1, z11 -> zR0*zR10, zR1 -> 1, gam -> 0, gamR -> 0} 

does not actually solve the original polynomial system, as can be verified by substituting these values back in the system: you could do this in M2 with 

sub(I, {z1 => 1, z11 => zR0*zR10, zR1 => 1, gam => 0, gamR => 0})

or an equivalent in Mathematica. Of the 24 components that Mathematica gives, only 6 of these (namely the 2nd, 4th, 13th, 15th, 21st, 22nd components) are actually contained in the true solution set: 4 of them are the components that M2 gives, and the other 2 are strictly smaller subsets of one of the 4 true components (and thus are not actually components).

Arrigo Coen

unread,
Nov 2, 2020, 3:33:51 PM11/2/20
to Macaulay2
Thank you, Justin. All your comments help me a lot!

I think I now understand enough to solve my equations with Macaulay2. 

(Pd. Sorry for the delay in this message. I don't know why it wasn't sent three days ago)

Reply all
Reply to author
Forward
0 new messages