[SciPy-user] Orthogonal polynomial evaluation

Pauli Virtanen pav@iki...
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.

Ticket: http://projects.scipy.org/scipy/ticket/921

Pauli Virtanen

More information about the SciPy-user mailing list