📈Interpolation and Approximation

Chebyshev Points and Approximation

Runge's phenomenon, the node polynomial, Chebyshev points, Lebesgue constants and spectral accuracy

Hannover, 1901

Carl Runge — remembered today for the Runge–Kutta methods of Unit 9 — publishes a short paper, Über empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten: on interpolation between equally spaced ordinates. His example is innocent. The function

f(x) = 1/(1 + 25x^2)

[-1, 1]

is smooth, positive and bell-shaped, with every derivative defined everywhere. Interpolate it at 6, 11, 21, 41 equally spaced points. Near the centre the polynomial hugs the curve better and better. Near the ends it begins to swing above and below, and the swings grow: maximum errors of 0.43, 1.9, 60 and

10^5

. More data, worse answer — without any rounding error at all. This is a property of the mathematics, not the machine. And the cure turned out to be not fewer points, but different points.

The error formula of lesson 7-1 explains it. The factor

\prod(x - x_j)

— the node polynomial — depends only on where the nodes are. For equally spaced nodes it is tiny in the middle of the interval and enormous near the ends: at degree 20 its size near

\pm 1

is thousands of times its size near 0. Runge's function has derivatives that grow like

n!\,5^n

, which the

(n+1)!

cannot beat, so the large end values of the node polynomial win.

💡Chebyshev points

Choose the nodes to make

\max|\prod(x - x_j)|

as small as possible. The answer is the zeros of the Chebyshev polynomial

T_{n+1}(x) = \cos((n+1)\arccos x)

x_j = \cos\dfrac{(2j+1)\pi}{2n+2}, \qquad j = 0, \dots, n,

for which

\max_{[-1,1]}|\prod(x - x_j)| = 2^{-n}

, the smallest possible. They are equally spaced in angle on a semicircle, projected down: clustered towards the ends of the interval, exactly where equispaced points let the error escape. (The extrema

\cos(j\pi/n)

— 'Chebyshev points of the second kind' — work equally well and include the endpoints.) Any interval

[a, b]

is mapped to

[-1, 1]

first.

Derivation

Chebyshev polynomials, their zeros, and why nothing beats

2^{-n}

They are polynomials. Put

x = \cos\theta

and

T_k(x) = \cos k\theta

. The identity

\cos((k+1)\theta) + \cos((k-1)\theta) = 2\cos\theta\cos k\theta

becomes

T_{k+1}(x) = 2x\,T_k(x) - T_{k-1}(x), \qquad T_0 = 1,\ T_1 = x,

T_2 = 2x^2 - 1

T_3 = 4x^3 - 3x

, and every

T_k

is a polynomial. Each step doubles the leading coefficient:

T_{n+1}

starts with

2^n x^{n+1}

. And

|T_k(x)| = |\cos k\theta| \le 1

[-1, 1]

Its zeros.

\cos((n+1)\theta) = 0

when

(n+1)\theta = (2j+1)\pi/2

, which gives

x_j = \cos\big((2j+1)\pi/(2n+2)\big)

j = 0, \dots, n

— the Chebyshev points.

The node polynomial. With these nodes

\prod_j(x - x_j)

has the same roots as

T_{n+1}

and leading coefficient 1, so it equals

2^{-n}T_{n+1}(x)

, whose largest size on

[-1, 1]

2^{-n}

Why no other nodes do better.

2^{-n}T_{n+1}

reaches

\pm2^{-n}

alternately at

n + 2

points,

\cos(k\pi/(n+1))

. Suppose some other monic polynomial

q

of degree

n + 1

stayed strictly inside

(-2^{-n}, 2^{-n})

. Then

2^{-n}T_{n+1} - q

would be positive at the first of those points, negative at the next, and so on —

n + 1

sign changes, hence

n + 1

roots. But the leading terms cancel, so it has degree at most

n

, and must be identically zero: a contradiction.

Degree n	Equispaced nodes	Chebyshev nodes
5	0.43	0.56
10	1.9	0.11
20	60	0.015
40	1.0 × 10⁵	2.9 × 10⁻⁴

Maximum interpolation error for Runge's function 1/(1 + 25x²)

Experience never errs; it is only your judgments that err by promising themselves effects such as are not caused by your experiments.

– Leonardo da Vinci, Notebooks

💡The Lebesgue constant

Interpolation is a linear map from data to polynomial. Its norm,

\Lambda_n = \max_x \sum_j |\ell_j(x)|,

the Lebesgue constant, bounds how much errors in the data — measurement noise, rounding — can be amplified, and how far the interpolant can be from the best possible polynomial approximation:

\|f - p\| \le (1 + \Lambda_n)\|f - p^*\|

. For Chebyshev points

\Lambda_n \approx \tfrac{2}{\pi}\ln n + 1

: 2.9 at

n = 20

. For equispaced points it grows exponentially: about 11,000 at

n = 20

and 6.6 million at

n = 30

. Data accurate to 0.001 at 21 equispaced points can give an interpolant wrong by 2.5.

Derivation

Where the factor

1 + \Lambda_n

comes from

Write

I

for the interpolation map: it takes the values of a function at the nodes and returns the polynomial

Ig = \sum_j g(x_j)\ell_j

. Two properties are all we need.

It amplifies by at most $\Lambda_n$ . At any point,

|(Ig)(x)| \le \max|g|\sum_j|\ell_j(x)| \le \Lambda_n\max|g|

It reproduces polynomials. If

q

has degree

\le n

, its interpolant is

q

itself (uniqueness, lesson 7-1).

Now let

p^*

be the best polynomial approximation of degree

n

f

. Since

Ip^* = p^*

f - If = (f - p^*) - I(f - p^*) \quad\Longrightarrow\quad \|f - If\| \le \|f - p^*\| + \Lambda_n\|f - p^*\| = (1 + \Lambda_n)\,\|f - p^*\|.

Interpolation is never worse than the best possible approximation by more than the factor

1 + \Lambda_n

. With Chebyshev points that factor is about 4 at

n = 20

; with equispaced points, about 11,000. And because data errors

\delta_j

produce the polynomial

\sum\delta_j\ell_j

, the same

\Lambda_n

bounds how much noise in the data is amplified.

Spectral accuracy. At Chebyshev points the interpolant is essentially the truncated Chebyshev series

\sum c_kT_k(x)

, and for smooth functions the coefficients

c_k

fall geometrically — the error drops by a constant factor with every extra degree. A smooth specific-heat curve on 300–3000 K is reproduced to

1.7 \times 10^{-2}

J/kg/K at degree 8 and

5.6 \times 10^{-5}

at degree 16, with coefficients falling from

10^3

10^{-5}

. This is how libraries evaluate special functions, how flight software stores ephemerides of the Moon and planets (JPL's DE ephemerides are Chebyshev coefficients over short intervals), and the idea behind the Chebfun software that computes with functions 'to machine precision'. A function with a kink like

|x|

converges only algebraically — polynomials should not be asked to represent corners.

💡Where you will use this

Orbital Mechanics and GNSS: planetary and lunar ephemerides (JPL DE440) are stored as Chebyshev coefficients. Propulsion: NASA's thermodynamic polynomials for species properties are fitted piecewise over temperature ranges for exactly the reasons of this lesson. Optimal Guidance: pseudospectral trajectory optimisation places its collocation points at Chebyshev or Legendre nodes. Instrumentation: sample where the function changes fastest, not uniformly by habit.

📚Further reading

Lloyd N. Trefethen, Approximation Theory and Approximation Practice (SIAM, 2013; extended ed. 2019), chapters 2–3, 13 and 15 — Chebyshev points and polynomials, Runge's phenomenon and Lebesgue constants, each a short chapter with its own experiments.
Carl Runge, 'Über empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten', Zeitschrift für Mathematik und Physik 46, 224–243 (1901) — the paper of the cold open.
Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 7 — Chebyshev points and orthogonal polynomials.

↓But there is a problem

Chebyshev points rescue polynomial interpolation when we can choose where to sample. Usually we cannot: the wind tunnel ran at the angles it ran at, the survey points are where the surveyor stood, the cam profile was measured every 30 degrees. And a single high-degree polynomial through scattered data is rarely what an engineer wants anyway. Shipwrights and draughtsmen had solved the practical problem centuries earlier. They bent a thin, flexible strip of wood — a spline — through pegs at the desired points, held it with lead weights called ducks, and traced the smooth curve it took. The mathematics of that curve went into car factories in 1959.