Ideal.variety() when working with symbolic rings
I have a system of multivariate polynomials I would like to solve (the equations come from chemical reaction networks). I have many symbolic constants, but I've managed to generate a Groebner basis with the following code (using a toy example with solutions {{x: x0, y:0},{x: 0,y: x0}):
x0 = var('x0') # In practice, there are many constants
S.<x,y> = PolynomialRing(SR,order='lex'); # in practice, there are 3-10 variables
I = S.ideal(
x**2 + y**2 - x0**2,
x + y - x0
)
print(I.dimension()) # dimension should always be zero
# >> 0
print(I.groebner_basis())
# >> [x + y - x0,
# >> y^2 + (-x0)*y]
In practice, I get many resulting equations which are extremely long, and it is numerically challenging to sequentially find the roots of these polynomials (in downstream code, I will get values for the constant parameters and then need to return solutions to these equations on the fly). Using the I.variety() method would be great to symbolically isolate all possible solutions, but I get confusing internal errors when I try to evaluate it.
print(I.variety()) # should be equivalent to: I.variety(ring=SR)
#/Users/ksb/SageMath/local/lib/python2.7/site-packages/sage/rings/polynomial/multi_polynomial_ideal.pyc in _variety(T, V, v)
# 2340
# 2341 variable = f.variable(0)
#-> 2342 roots = f.univariate_polynomial().roots(ring=ring, multiplicities=False)
# 2343
# 2344 for root in roots:
#
#/Users/ksb/SageMath/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_element.pyx in sage.rings.polynomial.polynomial_element.Polynomial.roots (build/cythonized/sage/rings/polynomial/polynomial_element.c:68811)()
# 7714 return l
# 7715 else:
#-> 7716 return [val for val,m in l]
# 7717 vname = 'do_not_use_this_name_in_a_polynomial_coefficient'
# 7718 var = SR(vname)
# TypeError: 'sage.symbolic.expression.Expression' object is not iterable
Is this expected behavior? When I try f.univariate_polynomial().roots(ring=ring, multiplicities=False) (setting multiplicities=True is what I.variety() does and leads to the error) I get [1/2*x0,2] for roots, which makes little sense (it seems like the roots should be [0,-x0]).
I would be happy to convert the MPolynomial_polydict (or the sage.rings.polynomial.polynomial_ring.PolynomialRing_field_with_category.element_class generated by univariate_polynomial()) into a symbolic expression that I could just set it equal to zero and use solve on it, but I don't know how to do this type conversion.
copying my comments to tmoteil here:
A less trivial problem looks like this:
When I do this, I get an empty list as the result.
I suppose the problem may be that there are no solutions that are valid over all possible values of the constants? Here's an example of there existing a solution if we substitute 1.0 (or 0.5) for all the constants:
The solution is {1/3,1/3,1/3}
Secretly I know that the value for every constant will be a positive number. And, while there may be lots of possible (x9,x10,x11) triples that satisfy the equations, the only ones that are "physically meaningful" are those with all three values between 0 and 1 (inclusive), possibly irrational. The fact that I cannot find any physically meaningful solutions after taking the Groebner basis and doing the rootfinding myself is what alerted me to a problem.
I'm not sure if these physically-motivated assumptions of mine are at the root of my confusion.