Skip to main content

Section 22.2 Chebyshev polynomials and interpolation

In Chapter 15 we studied Legendre polynomials \(P_n\) which are orthogonal on the interval \([-1, 1]\) in the sense that \(\int_{-1}^1 P_n(x)P_k(x)\,dx = 0\) when \(n\ne k\text{.}\) We also use Laguerre polynomials \(L_n\) which are orthogonal on \([0, \infty)\) with the weight function \(e^{-x}\text{,}\) meaning that \(\int_{-1}^1 L_n(x)L_k(x)\,e^{-x}\,dx = 0\) when \(n\ne k\text{.}\) The Chebyshev polynomials \(T_n\) are orthogonal on the interval \([-1, 1]\) with the weight function \(1/\sqrt{1-x^2}\text{,}\) meaning that

\begin{equation*} \int_{-1}^1 T_n(x)T_k(x)\,\frac{dx}{\sqrt{1-x^2}} = 0 \quad \text{when }n\ne k \end{equation*}

As for other families of orthogonal polynomials, we have a recursive formula for Chebyshev polynomials: starting with \(T_0(x)=1\) and \(T_1(x) =x\) we can compute the rest as

\begin{equation} T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x)\label{eq-chebyshev-recursion}\tag{22.2.1} \end{equation}

The recursion is simpler than for \(P_n\) or \(L_n\text{:}\) there is no division, so all polynomials \(T_n\) have integer coefficients.

Recursively find the coefficients of Chebyshev polynomials \(T_n\) for \(n\le 6\text{.}\) Plot all of them on the interval \([-1, 1]\text{.}\)

Answer

This is a straightforward modification of Example 15.4.1.

p = [1];
q = [1 0];
x = linspace(-1, 1, 1000);
hold on
plot(x, polyval(p, x), x, polyval(q, x));

for n = 1:5
    r = 2*[q 0] - [0 0 p];
    p = q;
    q = r;
    plot(x, polyval(r, x))
end
hold off

Example 22.2.1 indicates that Chebyshev polynomials oscillate between \(-1\) and \(1\) on the interval \([-1, 1]\text{.}\) The following remarkable identity confirms this observation:

\begin{equation} T_n(\cos\theta) = \cos (n\theta) \label{eq-chebyshev-cosine}\tag{22.2.2} \end{equation}

for all \(\theta\in\mathbb R\text{.}\) The proof is by induction. Base of induction is \(n=0, 1\text{:}\) for \(T_0(x)=1\) and \(T_1(x)=x\) the validity of (22.2.2) is obvious. For the inductive step, assume (22.2.2) holds for \(T_0, \dots, T_n\) and use (22.2.1):

\begin{equation*} T_{n+1}(\cos \theta) = 2\cos \theta \cos (n\theta) - \cos ((n-1)\theta) = \cos ((n+1)\theta) \end{equation*}

where the second step involves the identity

\begin{equation*} 2\cos \alpha \cos\beta = \cos (\alpha-\beta) + \cos(\alpha+\beta) \end{equation*}

The right hand side of (22.2.2) is zero when \(n\theta\) is an odd multiple of \(\pi/2\text{,}\) hence \(\theta = (2k-1)\frac{\pi}{2n}\text{.}\) Plugging \(k=1, \dots, n\) we find that

\begin{equation} T_n(x) = 0 \quad \text{when } x = \cos \left((2k-1)\frac{\pi}{2n}\right), \ k=1, \dots, n\label{eq-chebyshev-roots}\tag{22.2.3} \end{equation}

and since the polynomial \(T_n\) can have at most \(n\) real roots, these are all of its roots. They are easy to visualize: draw a semi-circle with interval \([a, b]\) as diameter, divide it into \(n\) equal arcs, and project the midpoint of each arc back to the interval.

Figure 22.2.2. Chebyshev points of the first kind

The recursive formula shows that the leading coefficient of \(T_n\) is \(2^{n-1}\text{.}\) Therefore,

\begin{equation*} T_n(x) = 2^{n-1}\prod_{k=1}^n (x-x_k) \end{equation*}

where \(x_1, \dots, x_n\) are the roots (22.2.3). It follows that the absolute value of the product \(\prod_{k=1}^n (x-x_k)\) is bounded by \(2^{1-n}\text{.}\) It can be proved that \(2^{1-n}\) is as small as one can get for any choice of \(n\) points on the interval \([-1, 1]\text{.}\) This suggests that Chebyshev points should work well for polynomial interpolation, and they do. Let us re-do Example 21.2.2 with random data to illustrate this.

Let \(x_1, \dots, x_{10}\) be the roots of \(T_{10}\text{.}\) Choose \(y_k\) randomly from the standard normal distribution, and draw the interpolating polynomial through the points \((x_k, y_k)\text{.}\)

Answer

The only change to Example 21.2.2 is replacing the x-coordinates.

n = 10;
x = cos((2*(1:n)-1)*pi/(2*n)); 
y = randn(1, n);
w = ones(1, n);
for k = 1:n
    for j = 1:n
        if j ~= k
            w(k) = w(k)/(x(k)-x(j));
        end
    end
end

p = @(t) (y.*w*(t - x').^(-1)) ./ (w*(t - x').^(-1));
t = linspace(min(x)-0.01, max(x)+0.01, 1000);
plot(t, p(t), x, y, 'r*')

Running the code Example 22.2.3 several times, we see that the interpolating polynomial behaves reasonably, without unnatural oscillations that come from interpolating at equally spaced points.