ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Thu, 16 Jan 2014 14:17:21 +0100Bug in Legendre polynomial.https://ask.sagemath.org/question/10931/bug-in-legendre-polynomial/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. Thu, 16 Jan 2014 13:38:01 +0100https://ask.sagemath.org/question/10931/bug-in-legendre-polynomial/Answer by DSM for <p>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>
<pre><code>P.<t> = QQ[]
legendre_P(90,t)(0.9)
</code></pre>
<p>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. </p>
https://ask.sagemath.org/question/10931/bug-in-legendre-polynomial/?answer=15927#post-id-15927The 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')
Thu, 16 Jan 2014 14:11:24 +0100https://ask.sagemath.org/question/10931/bug-in-legendre-polynomial/?answer=15927#post-id-15927Comment by Shashank for <p>The Legendre polynomial itself looks fine. To avoid the numerical imprecisions, you could work over QQ instead:</p>
<pre><code>sage: legendre_P(90,t)(0.9)
1.56319477427853e13
sage: legendre_P(90,t)(QQ(0.9)).n()
-0.0874931768791928
</code></pre>
<p>It looks like you need 100+ bits of precision in order to evaluate the function in the original order:</p>
<pre><code>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?
</code></pre>
<p>The polynomial is very oscillatory in that there's a large amount of cancellation going on.</p>
<p>If you're really interested in the value, and only secondarily the polynomial, then you could use either <code>scipy</code> or <code>mpmath</code>, depending on how much precision you'll need:</p>
<pre><code>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')
</code></pre>
https://ask.sagemath.org/question/10931/bug-in-legendre-polynomial/?comment=16426#post-id-16426Thanks a lot. I am interested in the numerical value so scipy function works for me.Thu, 16 Jan 2014 14:17:21 +0100https://ask.sagemath.org/question/10931/bug-in-legendre-polynomial/?comment=16426#post-id-16426