finding real zeros of polynomials [was "numerical value"]

This post is a wiki. Anyone with karma >750 is welcome to improve it.

Code:

for n in range (1,11):
print "n =", n
y = solve(gen_legendre_P(n, 0, x) == 0, x)
print "Zeros: {",
for i, eq in enumerate(y):
print eq.rhs().n(digits=2),
print ",",
print "}"


for n = 9 we get: Zeros: { -0.96 - 7.5e-6 * I , 0.96 + 7.5e-6 * I , -0.83 - 0.000039 * I , 0.83 + 0.000039 * I , -0.62 + 0.000052 * I , 0.62 - 0.000052 * I , -0.33 + 0.000024 * I , 0.33 - 0.000024 * I , 0.00 , }

and for n = 10 we get:

Zeros: {
Traceback (click to the left of this block for traceback)
...
TypeError: cannot evaluate symbolic expression numerically


1) How to get rid of for example. (-0.96 - 7.5e-6 * I) and get real numbers like (-0.96816) ?

2) Why n = 10 don't work?

edit retag close merge delete

Sort by » oldest newest most voted

This post is a wiki. Anyone with karma >750 is welcome to improve it.

"2) Why n = 10 don't work?"

Look at the equation it's trying to solve:

sage: solve(gen_legendre_P(10, 0, x) == 0, x)
[0 == 46189*x^10 - 109395*x^8 + 90090*x^6 - 30030*x^4 + 3465*x^2 - 63]


after substituting x^2->y, this is a 5th-order polynomial in y, and in general you can't solve an arbitrary polynomial higher than fourth order by radicals. This result is known as the Abel-Ruffini theorem.

The good news is that if you only care about numerical values, you can work in a polynomial ring over a real or complex field and get numerical results:

R.<x> = RealField(100)[]
for n in [1..10]:
p = gen_legendre_P(n, 0, x)
print n, p.real_roots()

1 [0.00000000000000000000000000000]
2 [-0.57735026918962576450914878050, 0.57735026918962576450914878050]
3 [-0.77459666924148337703585307996, 0.00000000000000000000000000000, 0.77459666924148337703585307996]
4 [-0.86113631159405257522394648889, -0.33998104358485626480266575910, 0.33998104358485626480266575910, 0.86113631159405257522394648889]
5 [-0.90617984593866399279762687830, -0.53846931010568309103631442070, 0.00000000000000000000000000000, 0.53846931010568309103631442070, 0.90617984593866399279762687830]
6 [-0.93246951420315202781230155449, -0.66120938646626451366139959502, -0.23861918608319690863050172168, 0.23861918608319690863050172168, 0.66120938646626451366139959502, 0.93246951420315202781230155449]
7 [-0.94910791234275852452618968405, -0.74153118559939443986386477328, -0.40584515137739716690660641208, 0.00000000000000000000000000000, 0.40584515137739716690660641208, 0.74153118559939443986386477328, 0.94910791234275852452618968405]
8 [-0.96028985649753623168356086857, -0.79666647741362673959155393647, -0.52553240991632898581773904919, -0.18343464249564980493947614236, 0.18343464249564980493947614236, 0.52553240991632898581773904919, 0.79666647741362673959155393647, 0.96028985649753623168356086857]
9 [-0.96816023950762608983557620291, -0.83603110732663579429942978806, -0.61337143270059039730870203934, -0.32425342340380892903853801464, 0.00000000000000000000000000000, 0.32425342340380892903853801464, 0.61337143270059039730870203934, 0.83603110732663579429942978806, 0.96816023950762608983557620291]
10 [-0.97390652851717172007796401209, -0.86506336668898451073209668842, -0.67940956829902440623432736512, -0.43339539412924719079926594316, -0.14887433898163121088482600113, 0.14887433898163121088482600113, 0.43339539412924719079926594316, 0.67940956829902440623432736512, 0.86506336668898451073209668842, 0.97390652851717172007796401209]


The line R.<x> = RealField(100)[] generates a one-variable polynomial ring called R over (an approximation to) the real field, in this case with 100 bits of precision, and defines a variable "x" which lives there:

sage: R.<x> = RealField(100)[]
sage: R
Univariate Polynomial Ring in x over Real Field with 100 bits of precision
sage: x
x
sage: parent(x)
Univariate Polynomial Ring in x over Real Field with 100 bits of precision

more

This post is a wiki. Anyone with karma >750 is welcome to improve it.

the first one is easy replace

print eq.rhs().n(digits=2)


by

print eq.rhs().n(digits=2).real()


or

print eq.rhs().n(digits=2).abs()


The second one I don't know for some reason "solve" is not able to solve the equation.

more

I wonder why

( 2011-04-03 10:08:11 -0500 )edit