Revision history [back]
 | 1 | initial version |
Here is the complete working code!
# Import necessary modules
from sage.all import *
# Define the polynomial ring
R.<x, y, t> = PolynomialRing(QQ, order='lex')
# Define the target system and start system as polynomials
f = [x^2 + 4*y^2 - 4, 2*y^2 - x] # target system
g = [x^2 - 1, y^2 - 1] # start system
# Define the homotopy
h = [t*fi + (1-t)*gi for fi, gi in zip(f, g)]
# Expand the homotopy
eh = [expand(hi) for hi in h]
# Compute the Jacobian matrix
jh = matrix([[diff(ehi, xi) for xi in [x, y]] for ehi in eh])
# Discriminant system: homotopy equations + determinant of the Jacobian matrix
sys = eh + [jh.det()]
# Convert the system to polynomials in the polynomial ring
sys_polynomials = [R(equation) for equation in sys]
# Compute the Groebner basis with lexicographic ordering
I = Ideal(sys_polynomials)
gb = I.groebner_basis()
# Extract the discriminant polynomial
discriminant_poly = gb[-1]
# Convert the discriminant polynomial to a univariate polynomial by eliminating other variables
univariate_poly = discriminant_poly.univariate_polynomial()
print(univariate_poly)
# Solve the univariate polynomial
univariate_poly.roots()
