Ask Your Question
0

Bug in Legendre polynomial.

asked 2014-01-16 06:38:01 -0500

Shashank gravatar image

updated 2014-01-16 06:38:33 -0500

I was trying to evaluate Legendre polynomial $P_{n}(x)$ for a very large n, and it seems for n greater than 50 or so sage gives wrong answer. For eg.

P.<t> = QQ[]
legendre_P(90,t)(0.9)

gives 1.56319477427853e13 whereas wolfram alpha gives -0.08749. Is there an alternate way to calculate Legendre polynomials in sage, which would give the right answer.

edit retag flag offensive close merge delete

1 answer

Sort by ยป oldest newest most voted
0

answered 2014-01-16 07:11:24 -0500

DSM gravatar image

updated 2014-01-16 07:12:48 -0500

The Legendre polynomial itself looks fine. To avoid the numerical imprecisions, you could work over QQ instead:

sage: legendre_P(90,t)(0.9)
1.56319477427853e13
sage: legendre_P(90,t)(QQ(0.9)).n()
-0.0874931768791928

It looks like you need 100+ bits of precision in order to evaluate the function in the original order:

sage: legendre_P(90, t)(RealIntervalField(100)("0.9"))
0.?e1
sage: legendre_P(90, t)(RealIntervalField(120)("0.9"))
-0.08749?
sage: legendre_P(90, t)(RealIntervalField(150)("0.9"))
-0.08749317687919?

The polynomial is very oscillatory in that there's a large amount of cancellation going on.

If you're really interested in the value, and only secondarily the polynomial, then you could use either scipy or mpmath, depending on how much precision you'll need:

sage: import scipy.special
sage: scipy.special.eval_legendre(90, 0.9)
-0.08749317687919235
sage: import mpmath
sage: mpmath.legendre(90, 0.9)
mpf('-0.08749317687919235')
edit flag offensive delete link more

Comments

Thanks a lot. I am interested in the numerical value so scipy function works for me.

Shashank gravatar imageShashank ( 2014-01-16 07:17:21 -0500 )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

Stats

Asked: 2014-01-16 06:38:01 -0500

Seen: 93 times

Last updated: Jan 16 '14