Fri Feb 15 17:16:12 CST 2008
Fri, 15 Feb 2008 15:09:50 -0600, Ryan May wrote:
> Ryan May wrote:
> In talking with __pv on #scipy, it seems that the problem is likely just
> due to accumulated round off error in using the recurrence relations to
> implement the polynomial. Since the N-point window function needs an
> (N-1)th degree polynomial, you can see how this begins to really show
> the error. Is there any reason not to just implement chebyt as above
> within scipy.special? It would seem to me that a formula would always
> give superior results, but I could (admittedly) be naive about this.
I think the special.chebyt function in scipy does not use recurrence
relations for evaluating the polynomials at some point (which probably
would be numerically stable), but instead creates a numpy.poly1d object
containing the coefficients of the polynomial and uses them to calculate
results. As the magnitude of the coefficients increases as polynomial
degree increases, error from cancellations in floating-point arithmetic
becomes of the order of 1 for n > 60.
This limitation would be useful to mention in the docstrings of functions
A possible fix could be to implement a way to evaluate polynomials stably
in numpy.poly1d. (Horner scheme? Using the roots which are anyway known?
I guess there are several options.)
More information about the SciPy-user