An ellipse, not a disc

The function $1/(1 + 25x^2)$ is as well behaved on the real line as a function can be: positive, bounded, infinitely differentiable, a smooth bump centred at nought. Its Taylor series at nought converges only for $|x| < 0.2$. The reason, found in the radius of convergence, is not on the real line at all: the function has poles at $\pm i/5$ in the complex plane, a fifth of a unit from the centre, and a power series converges on the largest disc that avoids every singularity.

On the interval $[-1, 1]$ that is a disaster. The Taylor series covers a fifth of the interval and diverges on the rest. And yet polynomials approximate this function on the whole interval without any difficulty: the best polynomial of degree sixteen is within $0.02$ of it everywhere, and the error falls steadily as the degree rises. Whatever governs how well polynomials do on an interval, it is not the disc.

The figure below is what governs it. The same two poles, the same small Taylor disc — and an ellipse, with foci at the ends of the interval, passing through the poles. The ellipse contains the whole interval. The size of that ellipse is the interval’s version of the radius of convergence, and this essay is about why.

A polynomial on an interval is a cosine series in disguise #

The tool is the substitution that made Chebyshev’s points and Chebyshev’s polynomials appear in the first place: $x = \cos\theta$. As $\theta$ runs from $0$ to $\pi$, $x$ runs from $1$ back to $-1$, so a function on the interval becomes a function of an angle. Keep going round and the function of the angle repeats: it is periodic, and even, since $\cos(-\theta) = \cos\theta$.

The Chebyshev polynomials are defined so that this substitution turns them into cosines: $T_k(\cos\theta) = \cos k\theta$. So writing a function on $[-1, 1]$ as a sum of Chebyshev polynomials,

$$f(x) = a_0 + a_1 T_1(x) + a_2 T_2(x) + \cdots,$$

is exactly the same as writing $f(\cos\theta)$ as a sum of cosines, with the same coefficients. A Chebyshev series is a Fourier series, after a change of variable, and every fact about where Fourier coefficients come from and how fast they fall transfers without alteration.

The substitution is not an arbitrary trick. Equally spaced angles on a semicircle, dropped onto the diameter, are exactly the Chebyshev points, so the natural way to sample a function on an interval — the way that makes interpolation behave — is equally spaced in $\theta$ and bunched in $x$. Once the function is written in the variable in which its good sample points are evenly spaced, the tools of evenly spaced sampling are available: the discrete cosine transform computes all the coefficients from the samples at once, just as a discrete Fourier transform does for a periodic signal.

That is worth pausing on, because it is the surprising connection this essay rests on. Polynomials on an interval and trigonometric waves on a circle look like different subjects — one is algebra, one is harmonic analysis — and the substitution $x = \cos\theta$ makes them one subject. Fourier needed his series to solve the flow of heat; the same series, bent round a semicircle, is how a computer evaluates a function.

Across the whole interval, at once #

The Chebyshev series of $1/(1 + 25x^2)$ converges everywhere on $[-1, 1]$, and the convergence is uniform: the partial sums close in on the function across the interval together, not from the centre outward.

The worst errors are $0.36$, $0.074$ and $0.0068$ at degrees four, twelve and twenty-four. Between twelve and twenty-four the error falls by a factor of eleven in twelve degrees — about $1.22$ per degree. That number is the heart of the matter. Where does $1.22$ come from, when the function’s Taylor series has radius $0.2$?

The ellipse through the poles #

Under the substitution $x = \cos\theta$, let $\theta$ become complex. The point $e^{i\theta}$ then ranges over the complex plane, and $x = \cos\theta = \tfrac12 (w + 1/w)$ with $w = e^{i\theta}$. That map sends each circle $|w| = \rho$, for $\rho > 1$, to an ellipse — with foci at $\pm 1$, semi-major axis $\tfrac12(\rho + 1/\rho)$, semi-minor axis $\tfrac12(\rho - 1/\rho)$ — and it sends the unit circle itself to the interval, traversed twice.

So the Fourier series in $\theta$, whose coefficients decay like $\rho^{-k}$ exactly when the periodic function extends analytically into the annulus $1/\rho < |w| < \rho$, becomes a Chebyshev series whose coefficients decay like $\rho^{-k}$ exactly when $f$ extends analytically into the ellipse of parameter $\rho$. That is Sergei Bernstein’s theorem, from 1912: the Chebyshev coefficients of $f$ fall geometrically at rate $\rho$, where $\rho$ is the parameter of the largest ellipse with foci $\pm 1$ inside which $f$ is analytic, and the best polynomial error at degree $n$ falls at the same rate.

For $1/(1 + 25x^2)$ the poles are at $\pm i/5$, and the ellipse through them has $\rho = \tfrac15 + \sqrt{1 + \tfrac1{25}} = 1.2198$. The first figure draws it, and the rate of $1.22$ per degree measured from the partial sums is the same number, computed a different way.

Three functions, three ellipses, three rates. The function $1/(1 + x^2)$ has its poles at $\pm i$, five times further off, and its ellipse has $\rho = 1 + \sqrt 2 \approx 2.414$. The function $1/(1.2 - x)$ has its pole on the real axis just past the interval’s end, at $1.2$, and its ellipse is squeezed flat, with $\rho = 1.2 + \sqrt{0.44} \approx 1.863$. The rates fitted to sixty coefficients agree with the ellipses to four decimal places.

Notice what the ellipse does that the disc could not. A pole at distance $0.2$ from the middle of the interval gives an ellipse of parameter $1.22$, and a pole at distance $0.2$ beyond the end gives $1.86$ — the second is far easier, though it is just as close to the interval. The ellipse is thin near the ends and fat in the middle, so a singularity opposite the middle of the interval does the most damage and one beyond an end does the least. The disc treats all directions alike, because a power series is about one point; the ellipse knows the interval has ends.

The Runge function, reread #

The ellipse also finishes the story of the points that ruined the fit. That essay found that interpolating $1/(1 + 25x^2)$ at equally spaced points diverges near the ends, and that interpolating at Chebyshev points converges. The Chebyshev half now has its rate: the interpolant at $n + 1$ Chebyshev points converges at rate $\rho = 1.22$ per degree, like the series, because Chebyshev interpolation is nearly the same thing as cutting off the Chebyshev series.

The equally spaced half has a rate too, and it is the reason the failure happens. For equally spaced nodes the relevant curves are not ellipses but the level curves of a different potential, fatter near the ends of the interval, and the level curve through the interval’s ends encloses the poles at $\pm i/5$. So equally spaced interpolation “sees” the poles as lying inside the region where it needs analyticity, and diverges near the ends, where that region is widest. Carl Runge worked this out in 1901. The Chebyshev points are exactly the points whose curves are the Bernstein ellipses, and the ellipses are thinnest at the ends, where the equally spaced curves are widest.

The two failures met so far now have the same shape. A Taylor series fails outside the largest disc free of singularities. An interpolation scheme fails outside the largest level curve free of them. What differs is the family of curves, and the family is decided by where the information comes from: one point, or a spread of points on an interval.

A corner anywhere costs a power #

Geometric decay needs analyticity in some ellipse, however thin. A function with a corner, or a jump in some derivative, is not analytic in any neighbourhood of the interval at all, and its coefficients fall only like a power of $k$.

The rule is the Fourier rule, carried across by the substitution. $|x|$ has a corner; $f(\cos\theta) = |\cos\theta|$ has corners too, and the Fourier coefficients of a function with corners fall like $1/k^2$ — the same decay that makes a square wave’s ripples fall like $1/k$ one derivative lower. $|x|^3$ has two more continuous derivatives and falls like $1/k^4$.

The square root is the instructive case. $\sqrt{1 + x}$ is worse than $|x|$ in the obvious sense — its derivative is infinite at $x = -1$ — and yet its coefficients fall at the same rate, $1/k^2$. The substitution explains it: $\sqrt{1 + \cos\theta} = \sqrt2\,|\cos(\theta/2)|$, a function with an ordinary corner. The change of variable stretches the ends of the interval, since $x = \cos\theta$ moves slowly near $\theta = 0$ and $\pi$, and a singularity at an end is smoothed by the stretch into something milder. The same stretch is why Chebyshev points crowd at the ends: the variable $\theta$ is the natural one, and in $\theta$ the points are evenly spaced.

Nearly the best, for the price of a transform #

The error that keeps coming back to its worst found the best polynomial of each degree by an exchange algorithm, iterating until a lower and an upper bound met. Cutting off the Chebyshev series is far cheaper: one cosine transform of the function’s values at Chebyshev points gives every coefficient at once. The question is how much the cheapness costs.

Very little. At every degree up to sixteen, the cut-off series is within a factor of $1.67$ of the best error, and the two curves are parallel on a log scale, falling at the same rate $\rho$. In general the ratio is bounded by a constant plus a multiple of $\log n$ — the Lebesgue constant of the Chebyshev projection, which grows so slowly that at degree a thousand it is still below eight.

So in practice the best polynomial is almost never computed. A cut-off Chebyshev series, or equivalently the interpolant at Chebyshev points, gives an approximation that is within a digit of the best and costs a fast Fourier transform. The software system Chebfun, begun by Lloyd Trefethen and Zachary Battles in 2004, is built on nothing else: it represents each function as a Chebyshev series cut off where the coefficients reach rounding level, and computes with those series — adding, multiplying, integrating, finding roots — as though they were the functions themselves. The coefficient plots here are what it looks at to decide where to stop.

A point has circles, an interval has ellipses #

Why an ellipse, and why these foci? There is a way of seeing it that does not go through Fourier series at all, and it explains the disc in the same breath.

Polynomial approximation of a function near a set is controlled by how quickly a polynomial of degree $n$ can grow as it moves away from the set, given that it is bounded by one on the set itself. Near a single point, the polynomial $x^n$ is the extreme case: it is tiny near nought and grows like $r^n$ at distance $r$, so the curves of equal growth are circles about the point. Those circles are the discs of Taylor series, and moving the centre moves the circles with it.

Near an interval the extreme polynomial is the Chebyshev polynomial, which is bounded by one on $[-1, 1]$ and grows faster off it than any other polynomial of its degree with that bound. Its size at a complex point $z$ is about $\rho^n/2$, where $\rho$ is the parameter of the ellipse through $z$. So the ellipses are the level curves of polynomial growth away from the interval, exactly as circles are for a point. A singularity on the ellipse of parameter $\rho$ limits the approximation to an error that falls like $\rho^{-n}$, for the same reason a singularity at radius $r$ limits a power series.

In the language of potential theory, the circles and the ellipses are both level curves of the Green’s function of the set — the electrostatic potential of the set when it is a conductor charged to one volt, measured in the plane around it. For a point the potential is $\log|z|$ and its level curves are circles; for a segment it is $\log|z + \sqrt{z^2 - 1}|$ and its level curves are the ellipses with foci at the segment’s ends. Every set in the plane has such curves, and polynomial approximation on the set converges geometrically inside the largest one the function is analytic in. The disc and the ellipse are the two simplest cases of one theorem.

Entire functions, and no ellipse too large #

A function with no singularities anywhere in the plane — $e^x$, $\sin x$, any polynomial — is analytic inside every ellipse, however large. For such a function Bernstein’s theorem says the coefficients fall faster than any geometric rate, and they do. The Chebyshev coefficients of $e^x$ are $2I_k(1)$, where $I_k$ is a modified Bessel function, and they behave like $1/(2^k\,k!)$: at $k = 10$ the coefficient is about $5 \times 10^{-10}$, and at $k = 15$ it is about $5 \times 10^{-17}$, at the level of rounding. Sixteen terms give $e^x$ on $[-1, 1]$ to full double precision.

That is the same distinction one point’s worth of information drew between functions with an infinite radius of convergence and functions without, transferred from a point to an interval. And it is why the best degree-three approximation of $e^x$ in the previous essay was already accurate to half a per cent: an entire function is as easy as a function can be, and the best error falls by a larger factor at each degree than the one before.

The practical rule that falls out of all this is short. On an interval, smoothness buys speed in two currencies: a finite number of derivatives buys a power of $n$, and analyticity in an ellipse buys a geometric rate set by the ellipse. Which currency a function pays in is visible at a glance on a plot of its coefficients — a straight line on a log–log plot, or a straight line on a log plot.

Integrating what has been approximated #

A Chebyshev series is not only a way to evaluate a function. Each $T_k$ has an integral over $[-1, 1]$ that can be written down — nought for odd $k$, and $2/(1 - k^2)$ for even $k$ — so integrating the series term by term integrates the function, with an error no larger than the tail of the series. Evaluated from values at Chebyshev points, this is Clenshaw–Curtis quadrature, published in 1960, and its error falls at the rate $\rho$ of the function’s ellipse, exactly like the approximation.

The more famous rule, Gauss’s, chooses its points to integrate polynomials of twice the degree exactly and converges at rate $\rho^2$ — twice as many correct digits per point, on paper. In practice the two are close for most functions, a surprise documented by Trefethen in 2008, because both are limited by the same ellipse and the extra exactness of Gauss’s rule is spent on polynomials of high degree that a function analytic in a thin ellipse barely contains. The ellipse, not the rule, sets the pace.

What the coefficient plots cannot show #

The coefficients are computed by a cosine transform on 1,024 points (8,192 for the slowly decaying ones), which is exact for polynomials of degree below that and very accurate for the functions drawn. What the plots show is a rate, and a rate is a statement about infinitely many coefficients; sixty or two hundred are evidence for it, and the theorem is what turns the evidence into a fact.

The plots also stop at the floor of the arithmetic. Past about $10^{-16}$ the computed coefficients are rounding noise, and the rate cannot be read there — which is why the fit uses only the part of each sequence above that floor, and why the figure says where the floor is rather than letting the noise pass for mathematics.

And the ellipse picture assumes the singularities are known. For a rational function they are at the zeros of the denominator, and the figures compute $\rho$ from them directly. For a function given only by its values, the direction is reversed: the rate is measured from the coefficients and the ellipse is inferred, which is how Chebfun can report where a function it has never seen as a formula stops being analytic.

Still open: the nodes that interpolate best #

Chebyshev points are nearly optimal for interpolation: the factor by which their interpolant can magnify an error — the Lebesgue constant — grows like $\tfrac{2}{\pi}\log n$, and no set of nodes can do better than that rate. But they are not exactly optimal. At each degree there is a set of nodes with a slightly smaller constant.

Sergei Bernstein and Paul Erdős conjectured in the 1930s that the optimal nodes are the ones for which the magnification function has equal local maxima between consecutive nodes — the alternation idea once more, applied to the nodes rather than the polynomial. Thomas Kilgore, and independently Carl de Boor and Allan Pinkus, proved it in 1978. What is not known is a formula for those optimal nodes: they are characterised completely and computed numerically, degree by degree, and no closed description of them has been found, even though the Chebyshev points they so nearly equal have one of the simplest descriptions in the subject.

That is where this line of questions ends, for now. A power series lives on a disc, a polynomial approximation on an interval lives on an ellipse, and the best points to sample the interval at are an alternation problem with a known answer and no formula.