Numerical evaluation of polynomials in Chebyshev basis

Question

I have high order (15 and higher) polynomials defined in Chebyshev basis and need to evaluate them (for plotting) on some intervals inside the canonical interval $[1,\,-1]$. A good accuracy near 1 and -1, where Chebyshev polynomials change rapidly, is also required.

I know that there exists Clenshaw algorithm for evaluation of Chebyshev polynomials, which is somewhat similar to the Horner scheme.

And this is all that I know...

I also saw the question about similar problematics. The proposed solution is to use higher precision.

I wonder if there are some other methods for accurate evaluation of polynomials (particularly in Chebyshev basis) on the given intervals that don't heavily rely on extended precision calculations? Is it possible at least to improve the problem as much as possible, so that extended precision is needed only from really high polynomial orders?

I also thought about interpolation of a high order polynomial by lower order ones on the intervals of interest. However, I don't know if there's any systematic procedure for this approach.

Answer 1

Evaluating polynomials of arbitrarily large degree in a Chebyshev basis is practical, and provably numerically stable, using a barycentric interpolation formula. In this case, extended precision isn't needed, even for order 1,000,000 polynomials. See the first section of this paper and the references, or here (Myth #2) for more details. I'll summarize briefly. Let's say you have a polynomial $f$ in a Chebyshev basis, and you know its values at the Chebyshev nodes $$ f_j = f(x_j) $$ $$ x_j = \cos\left(\frac{\pi j}{N}\right), \;\; 0\leq j\leq N $$

Then for any $x$ which isn't one of the Chebyshev nodes $x_j$, we have $$ f(x) = \frac{\displaystyle \sum_{j=0}^N \frac{w_j}{x-x_j}f_j}{\displaystyle \sum_{j=0}^N \frac{w_j}{x-x_j}}, $$ where $$ w_j = \left\{ \begin{array}{cc} (-1)^j/2, & j=0\text{ or }j=N,\\ (-1)^j, & \text{otherwise} \end{array} \right. $$

I believe a naive implementation of the above for various values of $x$ provides a stable algorithm that does not suffer the numerical difficulties encountered when trying to sum up the Chebyshev polynomials directly. The key is to work with the representation of the function by its values, not by its coefficients in a Chebyshev basis.

Answer 2

You can use Java code performing arbitrary calculations in Chebyshev (also in Legendre, Laguerre, Hermite) basis, see Appendix A of a paper http://arxiv.org/abs/1510.05510 the source code is available from Ref [31] link in the paper. There is little documentation available, but class Chebyshev.java has a unit test, from which usage should be more or less clear.

Chebyshev B=new Chebyshev();

final double [] c=new double []{coef_T0,coef_T1,...}; // chebyshev basis coefs

for(double x=-1;x<1.0001;x+=0.1) { System.out.println("x="+x+" F="+B.calculateSumD0(c.length-1,c,x)); }