[SciPy-user] Legendre Polynomials
Anne Archibald
peridot.faceted@gmail....
Mon Aug 27 23:23:14 CDT 2007
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?
>
> As a test I tried to expand cos(x) over Legendre-Polynomials because I
> know the result: Coefficient 1 should be 1 and the others as small as
> possible!
>
> It works quite fine as long as l<33!
All of scipy's orthogonal polynomials are implemented using recurrence
relations. These are a fairly good way to address the problem in
general, but they can go berserk for high orders. I ran into a similar
problem when trying to evaluate high-order orthogonal polynomials of
several types (to construct my own with custom weights).
In some cases there are other approaches for evaluating them
accurately (for example the nth Chebyshev polynomial is
cos(n*arccos(X))); you might check Abramowitz and Stegun for some
helpful analytic relations. I think scipy.special might also have some
legendre functions, implemented using cephes, that might be more
reliably accurate.
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
More information about the SciPy-user
mailing list