First time here? Check out the FAQ!

Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

answered 11 years ago

DSM gravatar image

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')
click to hide/show revision 2
No.2 Revision

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')