# Chebyshev approximation by projection vs interpolation

Suppose we want to approximate a function $f: [a, b] \rightarrow \Re$ with a Chebyshev series:

$$f(x) \approx \sum_{k=0}^n c_k \, T_k\left( \frac{2x-b-a}{b-a} \right)$$

where $T_k(x) = \cos(k\, \cos^{-1}x)$ are Chebyshev polynomials of the first kind. There are two methods of finding the coefficients $c_k$ that make this approximation good.

### Projection / Truncation method

This uses the fact that Chebyshev polynomials are orthogonal. Let [A] = 1 if statement A is true and [A] = 0 otherwise:

\begin{align} \frac{\pi}{1 + [k \neq 0]} \, c_k &= \int_0^\pi f\left( \frac{b-a}{2} \cos\theta + \frac{b+a}{2} \right) \cos(k\theta) \, d\theta \\[1ex] &= \int_{-1}^1 f\left( \frac{b-a}{2}x + \frac{b+a}{2} \right) T_k(x) \, \frac{1}{\sqrt{1-x^2}} \, dx \\[1ex] &= \int_a^b f(x) \; T_k\left( \frac{2x-b-a}{b-a} \right) \frac{2}{\sqrt{(b-a)^2 - (2x-b-a)^2}} \, dx \end{align}

### Interpolation / Collocation method

The Chebyshev nodes in [–1, 1] are $x_k = \cos \displaystyle{\frac{k\pi}{n}}$ where k = 0, 1, ..., n. We can force the Chebyshev series to cut $f(x)$ at those nodes by solving the linear system

$$f\left( \frac{b-a}{2} \cos \left( \frac{i\pi}{n} \right) + \frac{b+a}{2} \right) = \sum_{j=0}^n c_j \, \cos \left( \frac{ij\pi}{n} \right) \qquad \forall \, i = 0, \dotsm, n$$

### Questions

1. Which method is better? In terms of accuracy, speed or any relevant criteria.

2. Which method is implemented in NumPy and Chebfun?

• Similar earlier question: scicomp.stackexchange.com/questions/23389 Also, your interpolation method doesn't require solving a linear system in the usual sense: it is, in fact, a discrete cosine transform. Jan 8 '17 at 18:25
• Which methods is more accurate depends on how you measure accuracy. Clearly, the $L_2$ project has the least $L_2$ error, but it may or may not have the smaller maximum norm, for example. Jan 9 '17 at 0:34
• @WolfgangBangerth Since the series is truncated, the error can only be roughly as small as the truncation error, which is independent of the method. Jan 9 '17 at 0:53
• Well, maybe indeed the right thing is to increase $n$ -- but that wasn't part of the question :-) Jan 11 '17 at 2:46
• @visitor: Why not just look at the source code for Numpy and Chebfun to find out?
– Paul
Jan 13 '17 at 16:53

Please don't downvote this answer just because it's incomplete. My intention is to let whoever answering my question build on it, rather than write from scratch. If your answer is more comprehensive than mine, then I'll mark yours as the answer.

Before answering your question, be aware of two kinds of Chebyshev nodes:

Roots of $T_{n+1}(x)$ are $\displaystyle{x_k = \cos \left( \frac{k+1/2}{n+1}\pi \right)}$ where k = 0, 1, ..., n

Extrema of $T_n(x)$ in [–1, 1] are $\displaystyle{x_k = \cos \left( \frac{k\pi}{n} \right)}$ where k = 0, 1, ..., n

The former may be necessary when the endpoints are problematic, such as when integrating a function that goes to infinity.

### 1. Which method is better?

The projection method can be made faster by using the fact that roots of $T_{n+1}(x)$ satisfy discrete orthogonality relations:

$$\sum_{k=0}^n T_r(x_k) T_s(x_k) = \left\{ \begin{array}{ll} 0 & \text{if } \, r \neq s \\ (n+1)/2 & \text{if } \, r = s \neq 0 \\ n+1 & \text{if } \, r = s = 0 \end{array} \right.$$

which can be proven by considering

\begin{align} S_n(\theta) &= \sum_{k=0}^n \cos (k+1/2) \theta = \frac{1}{2} \csc \frac{\theta}{2} \sin (n+1) \theta \\ S_n(0) &= n+1 \\ S_n \left( \frac{k\pi}{n+1} \right) &= 0 \qquad \qquad \forall \, k = \pm 1, \dotsm, \pm (2n+1) \\ S_n(\pm 2\pi) &= -(n+1) \end{align}

Therefore

\begin{align} \sum_{k=0}^n T_r(x_k) T_s(x_k) &= \sum_{k=0}^n \cos \left( \frac{k+1/2}{n+1} r\pi \right) \cos \left( \frac{k+1/2}{n+1} s\pi \right) \\ &= \frac{1}{2} \sum_{k=0}^n \left[ \cos \left( \frac{k+1/2}{n+1} (r+s) \pi \right) + \cos \left( \frac{k+1/2}{n+1} (r-s) \pi \right) \right] \\ &= \frac{1}{2} \left[ S_n \left( \frac{r+s}{n+1} \pi \right) + S_n \left( \frac{r-s}{n+1} \pi \right) \right] \end{align}

So to answer your question directly, computing the coefficients $c_r$ as follows will be faster than solving a linear system as in the interpolation method:

$$\frac{n+1}{1 + [r\neq 0]} \, c_r = \sum_{k=0}^n f\left( \frac{b-a}{2}x_k + \frac{b+a}{2} \right) T_r(x_k)$$

### 2. Which method is implemented in NumPy and Chebfun?

• The bit about "computing $c_r$ will be faster" I think is wrong: both the formula for $c_r$ and the interpolation formula in the question are examples of discrete cosine transform (DCT), and are evaluated using FFT in practice (e.g., in ApproxFun.jl). Evaluating the sums for $c_r$ directly would be unnecessarily slow, $O(n^2)$ instead of $O(n\log n)$, as would solving the linear system with LU/GE in time $O(n^3)$. This is also why, if you compute $c_r$ the way you've described, the methods produce identical results. Jan 14 '17 at 2:47
• Also, I think the $c_r$ in this answer are not quite the same as the integrals in the question, you seem to be ignoring higher-order terms in the Chebyshev series for $f$, which makes the different $c$'s not exactly equal. I think $f=T_{n+1}$ is a counterexample. Jan 14 '17 at 3:05