ksb's profile - activity

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.

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:

S.<x11,x10,x9> = PolynomialRing(QQ,order='lex');

I = S.ideal(
            -x10*x9 + x11**2 + x11 - x10,
            -x10*x9 + x11**2 + x11**2 - x9**2,
            x9 + x10 + x11 - 1
            )
print(I.dimension())
B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

print(I.variety())

The solution is {1/3,1/3,1/3}

copying my comments to tmoteil here:

A less trivial problem looks like this:

R.<x0,x1,x2,x3,x4,x5,x6,x7,x8> = PolynomialRing(QQ)
F=R.fraction_field()
S.<x9,x10,x11> = PolynomialRing(F,order='lex');

I = S.ideal(
            -x7*x10*x9 + x8*x2*x11**2 + x3*x0*x11 - x4*x10,
            -x7*x10*x9 + x8*x2*x11**2 + 2*x5*x1*x11**2 - 2*x6*x9**2,
            x9 + x10 + x11 - 1
            )
print(I.dimension()) # 0
B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

print(I.variety())

When I do this, I get an empty list as the result.

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.

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:

S.<x11,x10,x9> = PolynomialRing(QQ,order='lex');

I = S.ideal(
            -x10*x9 + x11**2 + x11 - x10,
            -x10*x9 + x11**2 + x11**2 - x9**2,
            x9 + x10 + x11 - 1
            )
print(I.dimension())
B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

print(I.variety())

The solution is {1/3,1/3,1/3}

I haven't yet adapted the very large chemical networks to this framework, but a less trivial problem looks like this:

R.<x0,x1,x2,x3,x4,x5,x6,x7,x8> = PolynomialRing(QQ)
F=R.fraction_field()
S.<x9,x10,x11> = PolynomialRing(F,order='lex');

I = S.ideal(
            -x7*x10*x9 + x8*x2*x11**2 + x3*x0*x11 - x4*x10,
            -x7*x10*x9 + x8*x2*x11**2 + 2*x5*x1*x11**2 - 2*x6*x9**2,
            x9 + x10 + x11 - 1
            )
print(I.dimension()) # 0
B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

print(I.variety())

When I do this, I get an empty list as the result.

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.

From this reference, it seems that a Groebner basis with respect to lexicographic order should yield a 'triangular system': www dot scholarpedia.org/article/Groebner_basis. My system has three equations and three unknowns, which I want to transform into a series of three 1-D root finding problems.

R.<x0,x1,x2,x3,x4,x5,x6,x7,x8> = PolynomialRing(QQ)
F=R.fraction_field()
S.<x9,x10,x11> = PolynomialRing(F,order='lex');

I = S.ideal(
        -x7*x10*x9 + x8*x2*x11**2 + x3*x0*x11 - x4*x10,
        -x7*x10*x9 + x8*x2*x11**2 + 2*x5*x1*x11**2 - 2*x6*x9**2,
        2*x7*x10*x9 - 2*x8*x2*x11**2 - x3*x0*x11 + x4*x10 - 2*x5*x1*x11**2 + 2*x6*x9**2
        )

B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

This yields four equations (rather than 3) which do not seem to exhibit this triangular property (first three equations mention all three variables x9,x10,x11, and the last one involves both x9 and x10).

x9^2 + ((-x4)/(2*x6))*x10 + ((-2*x1*x5)/(2*x6))*x11^2 + x0*x3/(2*x6)*x11

x9*x10 + x4/x7*x10 + ((-x2*x8)/x7)*x11^2 + ((-x0*x3)/x7)*x11

x9*x11^2 + 2*x0*x3/(2*x2*x8)*x9*x11 + ((-x4*x7)/(2*x2*x6*x8))*x10^2 + ((-2*x1*x5*x7)/(2*x2*x6*x8))*x10*x11^2 + x0*x3*x7/(2*x2*x6*x8)*x10*x11 + 2*x4^2/(2*x2*x7*x8)*x10 + ((-2*x4)/(2*x7))*x11^2 + ((-2*x0*x3*x4)/(2*x2*x7*x8))*x11

x10^3 + 2*x1*x5/x4*x10^2*x11^2 + ((-x0*x3)/x4)*x10^2*x11 + ((-2*x4*x6)/x7^2)*x10^2 + 4*x2*x6*x8/x7^2*x10*x11^2 + 4*x0*x3*x6/x7^2*x10*x11 + ((-2*x2^2*x6*x8^2)/(x4*x7^2))*x11^4 + ((-4*x0*x2*x3*x6*x8)/(x4*x7^2))*x11^3 + ((-2*x0^2*x3^2*x6)/(x4*x7^2))*x11^2

Do I have a misunderstanding about what to expect from Sage's output? Or is there a test to see if my system of equations is pathological in some way?

If I ignore the fact that variables x1...x8 are parameters, I am surprised to get a basis of three equations which is also not triangularizable with respect to x9,x10,x11.

S.<x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11> = PolynomialRing(QQ,order='lex');

I = S.ideal(
        -x7*x10*x9 + x8*x2*x11**2 + x3*x0*x11 - x4*x10,
        -x7*x10*x9 + x8*x2*x11**2 + 2*x5*x1*x11**2 - 2*x6*x9**2,
        2*x7*x10*x9 - 2*x8*x2*x11**2 - x3*x0*x11 + x4*x10 - 2*x5*x1*x11**2 + 2*x6*x9**2
        )

B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

>>x0*x3*x6*x9^2 + 1/2*x0*x3*x7*x9*x10 - x1*x4*x5*x10*x11 - x1*x5*x7*x9*x10*x11 - 1/2*x2*x4*x8*x10*x11 + x2*x6*x8*x9^2*x11

x0*x3*x11 + x2*x8*x11^2 - x4*x10 - x7*x9*x10

x1*x5*x11^2 + 1/2*x2*x8*x11^2 - x6*x9^2 - 1/2*x7*x9*x10