[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?
> >


> 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
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.

    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