Ask Your Question

Not understandable error when solving polynomial equation

asked 2016-04-10 18:29:46 +0200

wjansen gravatar image

updated 2016-04-11 16:08:51 +0200

tmonteil gravatar image

Hi all,

when I enter command

solve(symbolic_expression(x^12 - x^11 - 12*x^10 + 11*x^9 + 54*x^8 - 43*x^7 - 113*x^6 + 71*x^5 + 110*x^4 - 46*x^3 - 40*x^2 + 8*x + 1)==0, var(x), to_poly_solve=True)

I get the expected result, but when I enter command

solve(symbolic_expression(x^10 - 10*x^8 + 35*x^6 + x^5 - 50*x^4 - 5*x^3 + 25*x^2 + 5*x - 1), var(x), to_poly_solve=True)

I get the error message

TypeError                                 Traceback (most recent call last)
<ipython-input-203-6108bea90b72> in <module>()
----> 1 solve(symbolic_expression(x**Integer(10) - Integer(10)*x**Integer(8) + Integer(35)*x**Integer(6) + x**Integer(5) - Integer(50)*x**Integer(4) - Integer(5)*x**Integer(3) + Integer(25)*x**Integer(2) + Integer(5)*x - Integer(1)),var(x),to_poly_solve=True)

/usr/local/sage-6.4.1-x86_64-Linux/local/lib/python2.7/site-packages/sage/symbolic/ in solve(f, *args, **kwds)
    732     from sage.symbolic.expression import is_Expression
    733     if is_Expression(f): # f is a single expression
--> 734         ans = f.solve(*args,**kwds)
    735         return ans

/usr/local/sage-6.4.1-x86_64-Linux/local/lib/python2.7/site-packages/sage/symbolic/ in sage.symbolic.expression.Expression.solve (build/cythonized/sage/symbolic/expression.cpp:47061)()

/usr/local/sage-6.4.1-x86_64-Linux/local/lib/python2.7/site-packages/sage/symbolic/ in sage.symbolic.expression.Expression.solve (build/cythonized/sage/symbolic/expression.cpp:46887)()

TypeError: 'sage.symbolic.expression.Expression' object does not support indexing

What happened here? The error message is totally misleading (no index in the command!) and it is not to understand why the second command fails while the first works fine.

By the way, the polynomial in question has ten simple real roots, so there should be no problem to compute the roots if symbolic evaluation is not possible.

Thanks in advance Wolfgang

edit retag flag offensive close merge delete


To display blocks of code, either indent them with 4 spaces, or select the corresponding lines and click the "code" button (the icon with '101 010'). Can you edit your question to do that?

To display inline code, surround it within "backticks" or "backquotes" `.

slelievre gravatar imageslelievre ( 2016-04-11 13:14:14 +0200 )edit

2 Answers

Sort by » oldest newest most voted

answered 2016-04-11 16:05:21 +0200

B r u n o gravatar image

updated 2016-04-13 15:58:32 +0200

My answer has two parts. First, I show you what happens (there is bug that should be corrected! → this is now #20436). I am not able to fully explain why this happens (maybe a quite recent change in Maxima, but maybe not...). Then I give you a workaround (that is better to use I think even when the bug will be removed).

Let me denote by $p$ your degree-$12$ polynomial, and by $q$ the degree-$10$ polynomial. The command you typed makes a call to the software Maxima. It happens that Maxima is able to find the approximate roots of $p$, and returns something like [[x == -1.99194847020934], [x == -1.927925665494726], ..., [x == 1.967859308671922]] (though in its own language), and Sage correctly transforms it to its representation. Everything's fine! For $q$, Maxima is unable to solve it (I do not know why) and returns something like [0 == 2*x^5 - 10*x^3 + 10*x + sqrt(5) + 1, 0 == 2*x^5 - 10*x^3 + 10*x - sqrt(5) + 1]. But the problem is as follows: Sage expects that Maxima raises an error when it is unable to find explicit solutions. This is apparently not the case! (Here, this may come from a change in Maxima, I have to investigate further.) So since there is no error from Maxima, Sage thinks that Maxima returned explicit solutions and tries to convert. This is where the IndexError appears since the answer from Maxima has not the required form.

For your problem, there are alternative methods to use, that are better to my mind – that is faster, safer, with better error guarantees. You first declare a polynomial ring in x over the integers and the polynomial you want to work with:

sage: R.<x> = ZZ[]
sage: q = x^10 - 10*x^8 + 35*x^6 + x^5 - 50*x^4 - 5*x^3 + 25*x^2 + 5*x - 1

Then you can ask for the roots of your polynomials in RR which is the floating-point representation of the real numbers, or in AA which is the field of (exactly-represented) algebraic real numbers:

sage: q.roots(RR)
[(-1.98422940262896, 1),
 (-1.85955297177650, 1),
 (-1.27484797949738, 1),
 (-0.851558583130145, 1),
 (-0.374762629171449, 1),
 (0.125581039058627, 1),
 (1.07165358995799, 1),
 (1.45793725484282, 1),
 (1.75261336008773, 1),
 (1.93716632225726, 1)]
sage: q.roots(AA, multiplicities=False)

Note that though you only see a finite number of decimals in the case of AA, roots are actually exactly represented as algebraic numbers, meaning that you can ask for as many decimals as you want and they will always be correct:

sage: r = q.roots(AA, multiplicities=False)[0]
sage: r.n() # 53 bits of precision
sage: r.n(500) # 500 bits of precision
edit flag offensive delete link more


Nice sleuthing - if you can open a Trac ticket with at least the information at the top that would be awesome.

kcrisman gravatar imagekcrisman ( 2016-04-12 16:26:53 +0200 )edit

I'll do that.

B r u n o gravatar imageB r u n o ( 2016-04-12 17:19:14 +0200 )edit

Thanks for the answer. Now, I understand the error's background. I was not aware of calls like "some_polynomial.roots(Some_structure)". This way, things become really simple, in particular, since the polynomials had originally been generated in QQ['x'].

wjansen gravatar imagewjansen ( 2016-04-13 09:59:35 +0200 )edit

Note that you can explore such things. After you have defined a polynomial q, type q.roots? and you will get the documentation for this method, which will tell you about the optional arguments.

slelievre gravatar imageslelievre ( 2016-04-14 13:46:31 +0200 )edit

answered 2016-04-11 16:05:01 +0200

tmonteil gravatar image

updated 2016-04-11 16:08:19 +0200

I agree that there is some issue here, and it should somehow be translated into a bug report, though i am not sure of where is the error coming from. Anyway, thanks for reporting !

That said, if you want to work with polynomials, you whould work in a well defined polynomial ring, instead of a symbolic blob.

So first define the polynomial ring you want to work in:

sage: R.<x> = PolynomialRing(QQ)

This defines both the ring R and the undetermined x. Then you can define the polynomial P, and look for its roots in the rational field QQ (and check that no root is rational):

sage: P = x^10 - 10*x^8 + 35*x^6 + x^5 - 50*x^4 - 5*x^3 + 25*x^2 + 5*x - 1
sage: P.roots()

Or in the field of algebraic numbers (to get all of them as algebraic numbers):

sage: P.roots(QQbar)
[(-1.984229402628956?, 1),
 (-1.859552971776503?, 1),
 (-1.274847979497380?, 1),
 (-0.8515585831301453?, 1),
 (-0.3747626291714493?, 1),
 (0.1255810390586268?, 1),
 (1.071653589957994?, 1),
 (1.457937254842823?, 1),
 (1.752613360087727?, 1),
 (1.937166322257263?, 1)]

You can see that all roots are simple (the second number 1 stands for the multiplicity). If you don't care about multiplicities, you can ge rid of them as follows:

sage: P.roots(QQbar, multiplicities=False)

You can also ask for real roots as numerical floating-point approximations:

sage: P.roots(RDF)
[(-1.984229402628939, 1),
 (-1.85955297177652, 1),
 (-1.274847979497381, 1),
 (-0.8515585831301442, 1),
 (-0.3747626291714496, 1),
 (0.1255810390586268, 1),
 (1.0716535899579935, 1),
 (1.4579372548428182, 1),
 (1.7526133600877463, 1),
 (1.9371663222572437, 1)]

And you can check that in this case all roots are actually real (w.r.t. complex):

sage: len(P.roots(QQbar)) == len(P.roots(RDF))
edit flag offensive delete link more


Thanks for the answer. I was not aware of calls like "some_polynomial.roots(Some_structure)".

wjansen gravatar imagewjansen ( 2016-04-13 10:14:24 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2016-04-10 18:29:46 +0200

Seen: 735 times

Last updated: Apr 13 '16