[SciPy-user] Orthogonal polynomials

[Anne Archibald]

Hi,

I have been using the orthogonal polynomials implemented in
scipy.special. They are very convenient, but I have been encountering
numerical headaches. For example:

In [127]: chebyt(100)(0.5)
Out[127]: 3468.56206844

In [128]: cos(100*arccos(0.5))
Out[128]: -0.5

In some cases, I realize that there are numerical obstacles to
computing high-degree orthogonal polynomials. Certainly their
coefficients will be a problem. But in some cases (for example for
Chebyshev polynomials of the first kind) there are more accurate,
possibly more efficient, ways to compute them. These do not
necessarily allow one to compute the weights efficiently (although
root-finding does make it possible) but they would nevertheless be
useful.

Does any code like this exist in scipy? Is there any that would be
convenient to import? How difficult is it to add special-case code to
the orthogonal polynomial objects? I am willing to supply a patch for
the Chebyshev, at least, if I can do it cleanly.

I am particularly interested in the Legendre polynomials, but I can
get away with using the Chebyshev polynomials in my application.

Thanks,
Anne M. Archibald