Where to take the samples

Adding up rectangles defined the integral as what the sum of thin rectangles approaches, and measured how slowly it approaches: on a smooth curve, the left-hand sum of $n$ rectangles is off by an amount proportional to $1/n$, so a thousand rectangles buy three decimal places and no more. Area is the undoing of slope then gave the integral a second meaning, the undoing of a derivative, and noted that it usually cannot be used: most integrands have no antiderivative that can be written down, and the sums come back.

What comes back, though, is not one sum but a family. Every method of numerical integration in use is a quadrature rule, a weighted sum

$$\int_{-1}^{1} f(x)\,dx \approx \sum_{i=1}^{n} w_i\, f(x_i),$$

and the methods differ only in where they put the sample points $x_i$ and how much weight $w_i$ they give each. This essay compares four of them by measurement, and the measurements say that where the samples go matters more than how many there are.

Four ways to choose #

The trapezoid rule takes evenly spaced samples, joins neighbouring ones by straight lines, and adds the areas of the trapezoids underneath. Its weights are all equal except at the two ends, where they are halved. It is exact for straight lines and its error shrinks like $1/n^2$. Simpson’s rule fits a parabola through each consecutive triple of samples instead; its weights alternate $\tfrac13, \tfrac43, \tfrac23, \tfrac43, \ldots$, and it is exact, by a lucky cancellation, not only for parabolas but for cubics. Its error shrinks like $1/n^4$. These are the rules taught first, and in the figure they are the two straight lines on the logarithmic scale — a fixed power of $n$.

The lucky cancellation is worth a sentence, because it is the first sign of the idea Gauss pushed to its limit. On a single panel $[-h, h]$, Simpson’s rule samples at $-h$, $0$ and $h$, and any cubic can be split into a parabola plus a multiple of $x^3$. The parabola is integrated exactly by construction, and $x^3$ is odd, so its integral over the symmetric panel is nought — and so is the rule’s value for it, since the samples at $\pm h$ cancel and the sample at $0$ is nought. A rule designed for degree two gets degree three free, because its samples are placed symmetrically. Placement, not just the number of samples, buys exactness. The midpoint rule has the same gift one level down: one sample at the centre of each panel is designed to get constants right and gets straight lines right as well.

The other two do not space their samples evenly. Clenshaw–Curtis, published by Charles Clenshaw and Alan Curtis in 1960, samples at the Chebyshev points $\cos(k\pi/m)$, which crowd towards the two ends of the interval, and chooses weights that integrate exactly the polynomial through the samples. Gauss’s rule, from 1814, chooses both the points and the weights to make the rule exact for polynomials of the highest degree possible.

The figure’s integrand, $1/(1 + 16x^2)$, is the function that made Carl Runge’s name in the points that ruin the fit: perfectly smooth on the real line, but with poles close to it at $\pm i/4$, which is what makes equally spaced interpolation of it fail. Its integral is exactly $\tfrac12 \arctan 4$, so every error can be measured. With eighty samples the trapezoid rule is still wrong in the fifth decimal place and Simpson’s in the eighth. Gauss’s rule reaches the limit of double-precision arithmetic at seventy-one samples, and Clenshaw–Curtis follows it down, a little behind.

Twice as many conditions as samples #

Gauss’s idea is a count. A rule with $n$ samples has $n$ points and $n$ weights, $2n$ numbers in all. A polynomial of degree $2n - 1$ has $2n$ coefficients. So it is at least plausible that the $2n$ numbers can be chosen to integrate every polynomial of degree up to $2n - 1$ exactly — one condition for each power $1, x, x^2, \ldots, x^{2n-1}$ — and Gauss proved that they can, in exactly one way.

The table makes the count visible. With six samples, the trapezoid rule is exact for degree one and gets $x^2$ wrong already; Simpson’s rule, with its parabolas, is exact to degree three; Clenshaw–Curtis, which integrates the interpolating polynomial through its samples, is exact to degree five; and Gauss’s rule is exact all the way to degree eleven. Twelve is impossible for any rule with six samples, because the polynomial that is the square of $(x - x_1)\cdots(x - x_6)$ has degree twelve, is never negative and vanishes at every sample, so every rule with those samples assigns it the value nought while its integral is positive.

The same argument says where the samples must go. Let $p(x) = (x - x_1)\cdots(x - x_n)$ be the polynomial vanishing at the nodes. For any polynomial $q$ of degree less than $n$, the product $p\,q$ has degree less than $2n$ and vanishes at every node, so an exact rule assigns it nought, and therefore

$$\int_{-1}^{1} p(x)\, q(x)\, dx = 0.$$

The node polynomial must be perpendicular to every polynomial of lower degree, in the inner product given by the integral. That condition determines it: it is the Legendre polynomial $P_n$, which is what one subtraction clears a direction produced by applying Gram–Schmidt to $1, x, x^2, x^3$ without aiming at anything in particular. The samples go at its roots.

Two properties of the result are worth reading off the picture. The nodes are not evenly spaced: they crowd towards the ends of the interval, at a density proportional to $1/\sqrt{1 - x^2}$ — the same crowding as the Chebyshev points, and for the same reason, which is that a polynomial of high degree is most able to misbehave near the ends. And every weight is positive. That is not automatic for a rule with many samples; Newton–Cotes rules with many evenly spaced samples have weights of both signs and become useless as their degree rises, because they add large numbers that cancel. Gauss’s weights are positive and add up to two, the length of the interval, so the rule is a genuine average of the samples and cannot amplify the errors in them.

The nodes are eigenvalues #

Finding the roots of $P_n$ for large $n$ is a computation in its own right, and the way it is done in practice connects the rule to a different part of linear algebra. The Legendre polynomials obey a three-term recurrence,

$$(k + 1) P_{k+1}(x) = (2k + 1)\,x\,P_k(x) - k\,P_{k-1}(x),$$

which is how the figures evaluate them, and a recurrence of that shape says that multiplication by $x$, acting on the polynomials $P_0, \ldots, P_{n-1}$, is described by a symmetric tridiagonal matrix once the polynomials are scaled to unit length. Gene Golub and John Welsch observed in 1969 that the eigenvalues of that $n \times n$ matrix are exactly the roots of $P_n$, and that the weights are twice the squares of the first components of its unit eigenvectors.

So computing a Gauss rule is an eigenvalue problem for a symmetric matrix, and the facts that make it well behaved are the facts symmetry forces a right angle established: the eigenvalues of a symmetric matrix are real, so the nodes lie on the real line, and its eigenvectors are perpendicular, which is what makes the weights come out positive. The figures here take the older route — Newton’s method on the recurrence, started from an approximation of each root by a cosine — and check each node by evaluating $P_n$ there, but the two routes are computing the same numbers.

The eigenvalue picture also explains why the nodes interlace. The roots of $P_{n-1}$ separate those of $P_n$, one between each neighbouring pair, because the smaller matrix is the larger with its last row and column removed, and the eigenvalues of a symmetric matrix and of such a submatrix always interlace. Adding one sample to a Gauss rule therefore moves every existing sample, which is the practical drawback the comparison below turns on.

Exact for polynomials is not the point #

Exactness for polynomials is how the rules are designed, but it is not what the hero figure measures. The integrand there is not a polynomial, and the reason the well-placed rules converge geometrically on it is a theorem about approximation rather than about integration. If $f$ is analytic inside an ellipse in the complex plane with foci at $\pm 1$, the best polynomial approximations of degree $m$ to $f$ on the interval have errors shrinking like $\rho^{-m}$, where $\rho$ is the sum of the ellipse’s semi-axes. That is the content of an ellipse, not a disc. A rule exact for degree $d$ makes an error no larger than a fixed multiple of the best approximation of degree $d$, so Gauss’s rule, with $d = 2n - 1$, converges like $\rho^{-2n}$ and Clenshaw–Curtis, with $d = n$, would seem to converge only like $\rho^{-n}$.

For $1/(1 + 16x^2)$ the poles at $\pm i/4$ fix the ellipse, and $\rho = \tfrac14 + \sqrt{1 + 1/16} \approx 1.28$. Gauss’s errors should then fall by a factor of about $1.28^2 \approx 1.64$ per sample, which is a decade every 4.7 samples; the hero figure’s Gauss line falls by about that much. But Clenshaw–Curtis does not lag by the factor of two that its degree suggests. Its line is nearly parallel to Gauss’s.

Gauss against Clenshaw–Curtis #

That near-equality is the subject of a 2008 paper by Nick Trefethen with the title Is Gauss quadrature better than Clenshaw–Curtis?, and the next figure repeats its comparison on four integrands.

For $e^x$, which is analytic everywhere, Gauss’s doubled degree is visible: it reaches $10^{-8}$ with five samples against eight. For Runge’s function the two lines nearly coincide, 38 samples against 39. For the two rough integrands, both rules slow to a power of $n$ and finish within a decade of each other. The factor of two in the degree of exactness buys a factor of two in the number of samples only for functions so smooth that both rules finish almost at once anyway.

Trefethen’s explanation is that the Clenshaw–Curtis error is controlled not by its degree but by how well the integrand’s Chebyshev series has converged, and for any function whose series coefficients decay slower than $\rho^{-2n}$ — any function with a singularity anywhere in the plane — the rule’s error tracks the coefficient at index $2n$ rather than $n$, the same as Gauss’s. Clenshaw–Curtis has practical advantages that this result makes decisive: its points are known in closed form, its weights come from a fast cosine transform, and the samples for $n$ points are reused when the number is doubled. Gauss’s points must be computed, and change completely when $n$ changes.

The trapezoid rule round a circle #

The trapezoid rule looks hopeless beside these, and on most integrands it is. But one class of integrand reverses the verdict entirely.

On $e^x$ over $[0, 1]$ the trapezoid rule behaves as the textbooks say: doubling the number of panels quarters the error. On $e^{\cos\theta}$ over a whole period, the same rule with the same evenly spaced samples reaches rounding error at fifteen panels. The difference is not that the second integrand is smoother — both are entire functions — but that it is periodic, and the trapezoid rule over a whole period has no ends.

The reason is clearest in Fourier’s language. Expand a periodic integrand as $\sum_k c_k e^{ik\theta}$. The exact integral over a period is $2\pi c_0$, since every other term integrates to nought. The trapezoid rule with $n$ equally spaced samples sums each term $e^{ik\theta}$ exactly to nought too — unless $k$ is a multiple of $n$, when every sample of $e^{ik\theta}$ equals one. So the rule computes $2\pi(c_0 + c_n + c_{-n} + c_{2n} + \cdots)$ and its error is the sum of the coefficients at multiples of $n$. That is aliasing, exactly as in a wave known only at its ticks: sampled $n$ times per turn, the wave of frequency $n$ is indistinguishable from a constant. For a smooth periodic function the Fourier coefficients fall off quickly, so the rule’s error is as small as the $n$-th coefficient.

The Euler–Maclaurin formula says the same thing from the other side. The trapezoid rule’s error on an interval is a series in powers of the panel width, with coefficients that involve the differences between the derivatives of $f$ at the two ends. On a whole period every one of those differences is nought, and every term of the series vanishes. What is left is smaller than any power of the panel width.

The rate is set by the nearest singularity #

How much smaller depends on how far the integrand can be continued into the complex plane before it breaks.

This is the same principle that sets a power series’ radius of convergence in one point’s worth of information: how fast an expansion converges is decided by the nearest place the function breaks, even when that place is off the line where the function is being used. The integrand $1/(a - \cos\theta)$ has poles where $\cos\theta = a$, at imaginary distance $\operatorname{arccosh} a$ from the real axis, and its Fourier coefficients decay like $(a + \sqrt{a^2 - 1})^{-|k|}$. The aliasing argument then predicts that each extra panel divides the error by $a + \sqrt{a^2 - 1}$, and the measured slopes match it: factors of 1.558, 2.618 and 5.828 per panel for the three values of $a$. Even the integrand with a pole close to the line, at $a = 1.1$, reaches full precision with eighty panels — a geometric rate, only a slow one.

This is the periodic counterpart of the ellipse for Gauss’s rule, and the two are connected by a substitution. Putting $x = \cos\theta$ turns an integral over $[-1, 1]$ into an integral of a periodic function over a half-turn, and turns the Chebyshev points into equally spaced angles. Clenshaw–Curtis is, essentially, the trapezoid rule applied after that substitution — which is why it converges geometrically, and why its points crowd towards the ends of the interval: they are equally spaced round the circle and crowd only when projected onto the diameter.

A corner that no choice of samples gets round #

Every rule above was rescued by smoothness. The last figure removes it.

The integrand $|x - 0.3|$ has a corner, and none of the four rules can integrate it faster than about $1/n^2$. Gauss’s rule, which was exact to degree $2n - 1$, is no better than the trapezoid rule, which was exact to degree one. The polynomial approximation that drove the geometric rates is now algebraically slow, because no polynomial can follow a corner closely, and the rules inherit that slowness whatever their degree. The errors also jump about as $n$ changes, because they depend on where the corner falls between the nearest samples, which is why the figure plots the worst error over each stretch rather than every value.

The remedy is not a cleverer placement but knowledge of the integrand: split the interval at $0.3$, and each half is a straight line, integrated exactly by any of the rules. Practical integration software does this adaptively, estimating the error on each piece and subdividing where it is large, which finds corners without being told where they are.

What these measurements do not settle #

Every integrand here has an exact integral, so every error is measured rather than estimated. That is not how quadrature is used. A real integrand comes without its answer, and a rule’s error must be estimated from the samples themselves — usually by comparing two rules, such as Gauss with $n$ points and a related rule with more, and treating their difference as the error of the less accurate one. Such estimates can be fooled. An integrand with a narrow spike between samples looks flat to every rule here, and both rules in the comparison agree on the wrong answer.

The comparisons are also one-dimensional. In two or three dimensions the product of $n$-point rules needs $n^2$ or $n^3$ samples, and in high dimensions no grid of any kind is affordable, which is where the error that does not care how many dimensions takes over: random samples, with an error of $1/\sqrt n$ whatever the dimension. The geometric rates in these figures are a one-dimensional luxury.

Still open: the best points for a class of functions #

Gauss’s rule is optimal for polynomials, and the comparison with Clenshaw–Curtis shows that being optimal for polynomials is not the same as being best for the functions people integrate. For a given class of functions — analytic in a given ellipse, or in a strip, or with a given number of derivatives — the rule with $n$ samples that minimises the worst-case error over the class is a well-defined object, and for some classes it is known: for functions analytic in a strip, the trapezoid rule after a suitable change of variables is close to optimal, which is the basis of the double-exponential rules of Hidetosi Takahasi and Masatake Mori. For most classes the optimal points have no formula, the same predicament as the best nodes have no formula for interpolation, and they are found, when they are found, by numerical optimisation.

A sharper question concerns the comparison itself. Trefethen’s analysis shows that Clenshaw–Curtis loses at most a factor of two in the number of samples, and usually nothing; whether there is a natural family of rules, cheap to construct like Clenshaw–Curtis, that matches Gauss’s rule on entire functions while keeping Clenshaw–Curtis’s advantages elsewhere has been studied in several forms, and no single answer has displaced the two old ones.