[SciPy-user] Legendre Polynomials
Kurt Smith
kwmsmith@gmail....
Tue Aug 28 01:18:39 CDT 2007
On 8/27/07, Anne Archibald <peridot.faceted@gmail.com> wrote:
> On 23/08/07, rgold@lsw.uni-heidelberg.de <rgold@lsw.uni-heidelberg.de> wrote:
>
> > A first quick (but not satisfying) idea was to use the values calculated
> > for x>=0 and copy them (of course correcting for the minus sign if l is
> > odd).
> > Problems remaining:
> > -> WHY is the evaluation wrong for x close to -1?
> > -> Can I trust the routine at all?
> >
[snippage]
>
> In terms of improving scipy, I don't think there's a better wayto
> handle arbitrary families of orthogonal polynomials. But it would be
> nice to have certain classes of orthogonal polynomial - chebyshev, for
> example - implemented using special-case, more reliable, methods. I
> haven't looked to see whether the orthogonal polynomial class (which
> is in pure python) admits this kind of extensibility.
>
> Anne
I took a stab at Anne's suggestion and tried a direct implementation
of Legendre's polys in terms of a series, not a recurrence. It
certainly isn't as general as the orthogonal poly class, but it seems
to solve the OP's accuracy problems.
For implementation, see
http://mathworld.wolfram.com/LegendrePolynomial.html
equation #33. I don't claim this is the best series to implement for
the problem, just the first one I tried. Use at your own risk!!!
import numpy as np
from scipy import factorial as fac
def nCk(n,k):
return fac(n,1) / fac(n-k,1) / fac(k,1)
def P_l(n):
""" P_l(n) -> legendre poly. of order n.
returns a (vectorized) function that evaluates
the legendre poly. of order n. Able to handle
x = -1.0 faithfully for large order.
Reference:
http://mathworld.wolfram.com/LegendrePolynomial.html, eqn. 33
"""
if n < 0:
raise ValueError("n must be >= 0")
front = 1.0/2.**n
coeff_arr = np.array([nCk(n,k)**2 for k in range(n+1)])
def inner_func(x):
if not -1.0 <= x <= 1.0: #perhaps do without this check
raise ValueError("x outside bounds")
xm1,xp1 = x-1., x+1.
val_arr = np.array(
[xm1**(n-k)*xp1**k for k in range(n+1)])
return front * np.dot(coeff_arr, val_arr)
return np.vectorize(inner_func)
from scipy.special import legendre as ssl
err = [ np.abs(ssl(k)(xx)-P_l(k)(xx)).sum() for k in range(33)]
# max(err) = 0.70110163376646795 for k == 32; error small for smaller
k, diverges for larger k...
endpts = np.array([P_l(k)(-1.0) for k in range(200)])
print endpts
(endpts[::2] == 1.0).all() # True.
(endpts[1::2] == -1.0).all() # True.
More information about the SciPy-user
mailing list