How to understand this Macaulay2 code?

224 Views Asked by At

I am trying to understand the m2 code (source) shown in the following:

restart
needsPackage "Visualize"
needsPackage "MinimalPrimes"
installMinprimes()
needsPackage "NumericalAlgebraicGeometry"
needsPackage "Bertini"

G = graph{{0,1}, {1,2}, {2,3}, {3,4}, {4,5}, {5,1}, {0,3}}

    openPort "8080"
    visualize G -- important: click on End session in browser window before continuing
    closePort()

makeIdeal = (kk) -> (
    S = kk[x_1..x_5,y_1..y_5];
    I = ideal(y_1+y_3,
        -x_2*y_1-x_5*y_1+x_1*y_2+x_1*y_5-y_1,
        x_2*y_1-x_1*y_2-x_3*y_2+x_2*y_3,
        x_3*y_2-x_2*y_3-x_4*y_3+x_3*y_4-y_3,
        x_4*y_3-x_3*y_4-x_5*y_4+x_4*y_5,
        x_5*y_1+x_5*y_4-x_1*y_5-x_4*y_5,
        x_1^2+y_1^2-1,
        x_2^2+y_2^2-1,
        x_3^2+y_3^2-1,
        x_4^2+y_4^2-1,
        x_5^2+y_5^2-1
        )
    )
I = makeIdeal QQ
S = ring I

netList I_*
gens gb I;
dim I -- so the dimension of the zero set in CC^14 is one.
degree I
C = decompose I -- there are 26 components defined over QQ.
C = C/trim;
netList C

C/dim//tally
C/degree

C1 = select(C, i -> dim i == 1)
C0 = select(C, i -> dim i == 0)
#C0
C0/degree
#C1
C1/degree
C1/degree//sum
for i in C1 list for j in C1 list (
    if i+j == 1 then " " else "+"
    )
netList(oo, Boxes=>false)

for i in C0 list ideal gens gb (i + minors(10, (jacobian I)%i))

  -- three of the points are singular (fat points).

-- The last curve is smooth.  The others are all singular, along 2 points of intersection of the pair of curves.
for i in C1 list ideal gens gb (i + minors(9, (jacobian I)%i))  
oo/decompose
flatten oo
oo/dim

C1_8
KK = frac(QQ[y_5])
SK = KK[x_1, y_1, x_2, y_2, x_3, y_3, x_4, y_4, x_5, MonomialOrder=>Lex]
P8 = sub(C1_8, SK)
G = flatten entries groebnerBasis P8
netList G

-- End of example....!

-- Numerically?

IC = makeIdeal CC 
IC = ideal drop(IC_*, {5,5})
SC = ring IC

-- These next two take too long right now...
elapsedTime numericalIrreducibleDecomposition IC
elapsedTime bertiniPosDimSolve IC_*

It is in general finding all isolated solutions to the polynomials $$y_1+y_3=0\\ -x_2*y_1-x_5*y_1+x_1*y_2+x_1*y_5-y_1=0\\ x_2*y_1-x_1*y_2-x_3*y_2+x_2*y_3=0\\ x_3*y_2-x_2*y_3-x_4*y_3+x_3*y_4-y_3=0\\ x_4*y_3-x_3*y_4-x_5*y_4+x_4*y_5=0\\ x_5*y_1+x_5*y_4-x_1*y_5-x_4*y_5=0\\ x_1^2+y_1^2-1=0\\ x_2^2+y_2^2-1=0\\ x_3^2+y_3^2-1=0\\ x_4^2+y_4^2-1=0\\ x_5^2+y_5^2-1=0$$

But I was not able to find which lines of code tell us the solutions. I am not familiar with algebraic geometry, but really want to have a quick understanding of how to solve polynomial equations by Macaulay2, under the help of the above example code. Could anyone help to briefly go through it?

Many thanks in advance to any help!

1

There are 1 best solutions below

1
On

I'm not sure where you got this code, but I think you're thinking too hard. Here's some (much shorter) code which also (numerically) computes the solutions to your family of polynomials:

-- this package exports the solveSystem function
needsPackage "NumericalAlgebraicGeometry"

-- build the polynomial ring
R = QQ[x_1..x_5,y_1..y_5]

-- this is a list of polynomials in R
system = { y_1+y_3,
           -x_2*y_1-x_5*y_1+x_1*y_2+x_1*y_5-y_1,
           x_2*y_1-x_1*y_2-x_3*y_2+x_2*y_3,
           x_3*y_2-x_2*y_3-x_4*y_3+x_3*y_4-y_3,
           x_4*y_3-x_3*y_4-x_5*y_4+x_4*y_5,
           x_5*y_1+x_5*y_4-x_1*y_5-x_4*y_5,
           x_1^2+y_1^2-1,
           x_2^2+y_2^2-1,
           x_3^2+y_3^2-1,
           x_4^2+y_4^2-1,
           x_5^2+y_5^2-1
         }

-- this does what it says on the tin
solveSystem system

This code runs in just under a minute for me, and outputs a list of (142, complex) solutions. For instance the first solution is:

{-.933179, 1, -.933179, .390388, .390388, -.359411, -2.49514e-14+3.63731e-14*ii, .359411, -.92065, .92065}

and as a quick sanity check, we can see that $y_1 = -.359411$ and $y_3 = .359411$ really do satisfy $y_1 + y_3 = 0$.

You can find documentation for this function here.


I hope this helps ^_^