[SciPy-user] Orthogonal polynomial evaluation
Mon Apr 13 12:35:19 CDT 2009
Mon, 13 Apr 2009 06:52:16 -0700, Lou Pecora wrote:
> I recall that stable routines to calculate Bessel functions use
> descending recurrence relations (descending in indices), rather than
> ascending recurrence which is unstable. Only trick is that the values
> have to be renormalized at the end (or as you go to keep the numbers in
> floating point range, I think) since descending recurrence does not set
> the scale initially. I wonder if this is the case for other recurrence
> relations. That is, if one recurrence is unstable (e.g. ascending),
> then the opposite (descending) will be stable. Perhaps the scipy
> routines can be easily reworked, if so. Just a thought.
Yes, the errors in polynomial evaluation come from a source somewhat
similar to the reason why Bessel recurrence relations have to be
evaluated in a certaion direction: loss of precision.
The problem here is that how the routines in Scipy currently work is that
they compute the polynomial coefficients and produce numpy.poly1d
objects. All poly1d objects evaluate the polynomials using the *same*
Horner recurrence scheme, but it is not numerically stable for high-order
polynomials. The problem is that each of the orthogonal polynomial would
need to be evaluated using a specialized recurrence relation.
We don't currently have implementations of stable evaluation algorithms.
So, a complete rewrite of the orthogonal polynomial routines is required.
More information about the SciPy-user