-- always include this first
restart;
-- Macaulay2 does include a package but it needs to be loaded
loadPackage "NumericalAlgebraicGeometry";
-- univariate complex polynomials
R = CC[z];
-- example 1: solving z^4-1=0
F1 = {z^4 - 1} -- the equation to solve
zs1 = solveSystem F1
xsol = sort(apply(zs1, zi -> (coordinates zi)#0))
<< "list of all solutions: " << xsol << endl
zr1 = realPoints zs1
xreal = sort(apply(zr1, zi -> (coordinates zi)#0))
<< "list of all real solutions: " << xreal << endl;
peek zs1#0
-- example 2: simple polynomial
F2 = {(z-5)^3*(z+2)^2}
zs2 = solveSystem F2
<< "solutions: " << zs2 << endl;
-- example 2: print the coordinates and status of sol. pts
xsol = apply(zs2, zi -> (coordinates zi)#0)
<< "solutions: " << xsol << endl;
<< apply(zs2, zi-> status zi) << endl;
-- use top-level Macaulay2 homotopy continuation routines
-- (instead of built-in ones, could also use Bertini)
zs = solveSystem F2, Software=>"M2"
-- example 3: z*(z-1)*...*(z-n) = 0
n = 5
F3 = {product(n+1, i -> z-i)}
zs = solveSystem(F3)
<< apply(zs, zi -> status zi) << endl;
<< apply(zs, zi -> (coordinates zi)#0) << endl;
-- example 3 with n=10 (try also n=20 and n=30)
n = 10
F3 = {product(n, i -> z-i)}
zs = solveSystem(F3)
<< apply(zs, zi -> status zi) << endl;
xsol = apply(zs, zi -> (coordinates zi)#0)
<< "solutions " << xsol << endl;
<< "number of sols "<< (length(xsol)) << endl;
<< "degree of polynomial "<< (degree(F3#0))#0 << endl;
<< "coeffs " << apply(listForm(F3#0), yi->yi#1) << endl;
<< "exponents " << apply(listForm(F3#0), yi->(yi#0)#0) << endl;
-- but for this polynomial it seems to work ...
-- (this is the polynomial used for the startsystem)
length(solveSystem({z^(11)-1})) -- just printing nu. of solutions
-- example 4 (a nicer real polynomial with singular (double) points)
n = 21
x = sort(apply(n, j -> cos(2*pi*j/n))) -- real parts of unit roots
F4 = {product(n, i -> z-x#i)}
zs = solveSystem(F4)
xsol = sort(apply(zs, zi -> (coordinates zi)#0))
<< apply(zs, zi -> status zi) << endl;
<< "error " << norm(xsol-x)/norm(xsol) << endl;
<< "residual "<<
norm(F4#0-product(xsol,xi -> z-xi))/norm(F4#0)<yi#1) << endl;
-- check the last solution and compare with the real part of
-- the corresponding start solution -- do this for different n
<< peek zs#(n-1) << endl;
<< "start sol " << x#(n-1) << endl;
-- example 5 (random solutions)
n = 6 -- works until around n=7 ...
x = sort(apply(n, j -> random(1.0)))
F5 = {product(x, xi -> z-xi)}
zs = solveSystem(F5)
xsol = sort(apply(zs, zi->(coordinates zi)#0))
<< apply(zs, zi -> status zi) << endl;
<< "error " << norm(xsol - x)/norm(x) << endl;
<< "residual "<<
norm(F5#0-product(xsol,xi -> z-xi))/norm(F5#0)<yi#1) << endl;
-- example 6 (random coefficients)
n = 40
c = apply(n+1, j -> random(1.0))
F6 = {sum(n+1, i -> c#i * z^i)}
zs = solveSystem(F6)
xsol = sort(apply(zs, zi->(coordinates zi)#0))
-- residual is difference between original and recovered polynomial
<< "residual "<< norm(F6#0 - c#n * product(xsol, xi->z-xi))/norm(F6#0) << endl;
<< "degree of polynomial "<< (degree(F6#0))#0 << endl;
<< "number of sols "<< (length(zs)) << endl;
<< "coeffs " << apply(listForm(F6#0), yi->yi#1) << endl;