[SciPy-user] Orthogonal polynomial evaluation

[Pauli Virtanen]

Mon, 13 Apr 2009 14:26:27 -0400, Anne Archibald wrote:
> 2009/4/13 Pauli Virtanen <pav at iki.fi>:
[clip]
>> We don't currently have implementations of stable evaluation
>> algorithms. So, a complete rewrite of the orthogonal polynomial
>> routines is required.
> 
> It's not just a question of polynomial evaluation schemes; the
> polynomial coefficients themselves present severe problems. I'd say that
> 1, x, ..., x**n is often a very poor basis for expressing polynomials.
> It would actually be nice to have a polynomials class that could use
> other bases - in fact families of orthogonal polynomials are natural
> bases for polynomials.

Yes, this was exactly my point: evaluating polynomials using the basic 
coefficients is not suitable for high-order polynomials.

> But let's leave aside poly1d entirely. There are algorithms, and I
> thought there was code in scipy, to evaluate orthogonal polynomials,
> their zeros, and their derivatives and integrals, in terms of recurrence
> relations. These are quite general, fairly efficient, and reasonably
> accurate - for small orders. But many orthogonal polynomials have
> special-purpose evaluation algorithms: for example the cos(n*arccos(x))
> expression for Chebyshev polynomials. Such algorithms are usually fast
> and accurate, but require separate implementation for each family that
> has one. I think making it possible to simply override the evaluation of
> a particular family of orthogonal polynomials would be a better solution
> to the OP's problem.

Perhaps I said it earlier in some strange way, but this is exactly the 
solution I proposed.

The problem is that the current code in orthogonal.py only computes the 
coefficients in the 1, x, x^2, ... basis, so we don't have specialized 
algorithms for evaluating the polynomials.

-- 
Pauli Virtanen