Numerical Integration and the Redemption of the Trapezoidal Rule

Numerical integration and the redemption of the trapezoidal rule 1 Numerical integration (“quadrature”) Freshman calculus revolves around differentiation and integration. Unfortunately, while you can almost always differentiate functions by hand (if the derivative exists at all), most functions cannot be integrated by hand in closed form [e.g. try integrating sin(x + cos(x))]. Instead, they must be integrated approximately on a computer, a process known as numerical integration or quadrature. (Historically, “quadrature” was a synonym for integration in general—literally, converting areas into equivalent squares— but in modern usage “quadrature” almost exclusively refers to computational algorithms.) In particular, suppose we are computing the following definite integral, over an interval [0, π] (chosen for convenience below):1 Z π I= f (x)dx 0 for some given function f (x) that is a “black box” (we may only know how to evaluate it given any x, not how to manipulate it symbolically). In numerical integration, we want to approximate this integral I by a sum IN : N I≈ ∑ wn f (xn ) = IN n=0 for N + 1 quadrature points xn ∈ [0, π] and corresponding quadrature weights wn . The central problem in numerical integration is this: • How can we choose the points xn and the weights wn so that the error |IN − I| goes to zero as rapidly as possible as we increase N? This problem has a long history, dating back to ancient approximations for π by approximating the area or circumference of a circle with polygons, and involves beautiful and 1 Note that there is no loss of generality in the choice of the [0, π] interval: if we are computing some R R arbitrary integral ab g(y)dy, we can convert it back to the [0, π] form by a change of variables: ab g(y)dy = R R |b−a| π |b−a| π x x 0 g( π [b − a] + a)dx = 0 f (x)dx for f (x) = π g( π [b − a] + a). π 1 (a) (b) f(x) f(x) errors error f(0) error f( ) trapezoid area = [f(0)+f( )]/2 ~ integral total trapezoid areas ~ integral x 0 0 x N=5 Figure 1: Illustration of (a) the trapezoidal rule and (b) the composite trapezoidal rule for integrating f (x) on [0, π]. In each case, we approximate the area under f (x) by the area of (a) one or (b) N trapezoids. That is, we evaluate f (x) at N + 1 points xn = nπ/N for n = 0, 1, . . . , N, connect the points by straight lines, and approximate the integral by the integral of this interpolated piecewise-linear function. sometimes surprising mathematics—the exponential blowup of polynomial (“Newton– Cotes”) approximations (Runge phenomena), orthogonal bases of polynomials (Gaussian quadrature, Chebyshev approximation, etc.), and deep forays into Fourier analysis (e.g., Clenshaw–Curtis quadrature). Quadrature is closely related to other important numerical algorithms such as numerical linear algebra (e.g. Lanczos iterations), approximation theory, and fast Fourier transform algorithms (FFTs, which themselves encompose a host of group theory, number theory, polynomial algebras, and other fascinating topics). For higher-dimensional numerical integration (cubature), the story becomes even more intricate, ranging from statistics (Monte–Carlo integration) and number theory (low-discrepancy sequences and quasi-Monte–Carlo methods) to fractals (sparse grids). Here, we will just dip our toes into the problem, beginning by analyzing one of the simplest—deceptively simple!—quadrature methods, the trapezoidal rule. Not only does a full analysis of the accuracy of this method lead us directly into the farreaching topic of Fourier series, but we also find that a simple transformation turns the lowly trapezoidal rule from one of the crudest quadrature schemes into one of the best, Clenshaw–Curtis quadrature. 2 The trapezoidal rule The trapezoidal rule, in its most basic form, connects the endpoints (0, f (0)) and (π, f (π)) by a straight line and approximates the area by the area of a trapezoid: f (0) + f (π) , 2 as shown in figure 1(a). Of course this approximation is rather crude, so we refine it by increasing the number of trapezoids: by “trapezoidal rule” one usually means a I≈π 2 x2 f(x) = (linear) + f'' f'' (quadratic) error area ~ x3 f'' x= /N Figure 2: Schematic estimate of the local error of the trapezoidal rule. In a small interval ∆x = π/N, the function f (x) can be approximated by a Taylor expansion, and the lowest-order correction to the trapezoidal rule’s linear approximation is represented by the quadratic term. Correspondingly, the maximum deviation of f (x) from the trapezoid is ∼ ∆x2 within the interval (multiplied by f 00 at some point). The resulting error area (red shaded region) is therefore proportional to ∆x2 × ∆x = ∆x3 . composite trapezoidal rule: divide [0, π] into N intervals and apply the trapezoidal rule to each one, as shown in figure 1(b). In the common case of equal intervals of width ∆x = π/N, summing these trapezoid areas yields the following approximate integral, also called the Euler–Maclaurin formula: " # π f (0) + f (π) N−1 IN = + ∑ f (nπ/N) . N 2 n=1 Note that the 1/2 factors cancelled except for the first and last points. Clearly, as N → ∞ we must have IN → I (at least, for any Riemann-integrable function). The question now is, how fast does the error |I − IN | decrease with N? 2.1 A simple, pessimistic error estimate A crude, but perhaps too pessimistic, upper bound on the error is as follows. The trapezoidal rule corresponds to approximating f (x) by a straight line on each interval ∆x = π/N. If we look at the Taylor expansion of f (x), the lowest-order deviation from a straight line is the quadratic ( f 00 ) term, and this term means that f deviates from a straight line by at most ∼ ∆x2 within the interval (multiplied by some coefficient proportional to f 00 ), as depicted in figure 2. The corresponding error area is then proportional to ∆x2 ×∆x = ∆x3 ∼ 1/N 3 . This is the local error from a single interval. As there are N such intervals, the total error should be bounded above by N × 1/N 3 = 1/N 2 : the error of the trapezoidal rule decreases at worst proportional to 1/N 2 (for continuous integrands). This is an upper bound on the error, because we have neglected the possibility that the errors from different intervals will be of opposite signs and mostly cancel if we are very lucky. At first glance, however, it seems unlikely that such cancellations will occur to such an extent that the error will decrease faster than 1/N 2 , and indeed this turns out to be the case—for most f (x), the trapezoidal-rule error is exactly proportional to 1/N 2 3 2 2 esin[(x+1) + 2 cos(4x+1)] esin[sin(x+1) + 2 cos(4x+1)] 3 1 0.5 0 -4 0 10 0.2 0.4 0.6 0.8 1 x/ -6 10 1/N 2 -8 10 ref e ren ce -10 fractional error in trapezoidal quadrature 1.5 -2 10 2.5 10 2 f(x) fractional error in trapezoidal quadrature 0 2.5 10 10 2 f(x) 0 3 -2 10 1.5 1 0.5 -4 10 0 0 0.2 0.4 0.6 0.8 1 x/ -6 10 -8 10 -10 10 -12 10 1 10 2 10 3 10 4 0 10 number N of trapezoids 10 20 30 40 50 number N of trapezoids Figure 3: Left: Fractional error |IN − I|/|I| verus N for the trapezoidal rule IN when 2 integrating the function esin[(x+1) +2 cos(4x+1)] (inset). For reference, the dashed line shows 1/N 2 dependence, demonstrating that the IN errors indeed decrease asymptotically at this rate. Right: the same fractional error, but this time for the integrand 2 esin[sin (x+1)+2 cos(4x+1)] (inset). Note that this is a semi-log scale: the IN errors appear to be decreasing at least exponentially fast with N—a miracle has occurred in the trapezoidal rule? for large N. However, it turns out there are some functions that do much, much better; if we can understand these special cases, and in general understand the error more rigorously, we can try to rearrange the computation so that this improvement occurs almost always. 2.2 A numerical experiment and a Miracle In figure 3(left) is plotted the fractional error |IN − I|/|I| in the trapezoidal-rule approximation versus N on a log–log scale for an arbitrarily chosen nasty-looking function 2 esin[(x+1) +2 cos(4x+1)] shown in the inset. As expected from the crude analysis above, the errors decrease asymptotically at a rate that is almost exactly proportional to 1/N 2 (an exact 1/N 2 dependence is shown as a dashed line for reference). If we want eight decimal places of accuracy we need around 104 function evaluations, but this is not too bad on a computer (unless we need to evaluate millions of such integrals, or unless our integrand is a much nastier function like the output of a planetary climate simulation!). 2 For “fun,” we try a slightly different integrand in the right panel of figure 3: esin[sin (x+1)+2 cos(4x+1)] (inset). Again plotting error versus N, it seems again to be roughly a straight line (a little more wiggly than before)—but wait, this is no longer a log–log scale, this is a log–linear scale. A straight line on such a scale indicates an exponential dependence with N. How has the ordinary trapezoidal rule, piecewise-linear approximations, man4 aged to achieve apparently exponential accuracy (or better2 ) with N? Some miracle of cancellation seems to have occurred, but why? If we could somehow manage to obtain this miracle for more integrands, what a wonderful quadrature method we would have! A clue to the miracle is that the integrand in figure 3(right) does have one special property that is obvious with a little thought: it is periodic [ f (x + π) = f (x)]. Why should periodicity help so much, though? That is probably not clear yet, but it points us towards an analytical technique where periodicity is central: Fourier analysis. 2.3 Error analysis by the Fourier cosine series A beautiful, powerful, and far-reaching way to rigorously analyze the approximation errors in the trapezoidal rule is to use Fourier analysis, which relies on an amazing fact: any “reasonable” function can be expressed as an infinite series of sines and/or cosines, and sines/cosines are much easier to analyze than arbitrary functions. In particular, for this problem it turns out to be convenient to use a Fourier cosine series: write f (x) as f (x) = ∞ a0 + ∑ ak cos(kx) 2 k=1 for coefficients ak given by: ak = 2 π Z π f (x) cos(kx) dx. 0 For any f (x) that is continuous and non-singular on [0, π], this series converges to f (x) at every x, and ak → 0 for k → ∞. It turns out that how fast the Fourier series converges (how fast ak → 0) is determined entirely by the smoothness of f (x) and by its endpoint derivatives at x = 0, π. Moreover, we will be able to describe the convergence rate of the trapezoidal rule entirely in terms of the convergence rate of the Fourier cosine series! 2.3.1 The exact integral I Plugging the cosine series into the exact integral I, we obtain a simple result: Z π I = 0 = ∞ π  a0 π  f (x)dx = + ∑ akcos(kx)dx 2 0 k=1 Z a0 π , 2 because all of the k > 0 terms integrate to zero. 2 In fact, the rate of decrease here may in fact be faster than exponential, a consequence of the niceness of this integrand everywhere in the complex x plane, but it is hard to tell the difference without plotting over a much larger range of errors (impossible here because of the finite precision of computer arithmetic). 5 2.3.2 The trapezoidal rule IN Plugging the cosine series into the trapezoidal rule IN , the constant a0 term just gets multiplied by N (it appears in N + 1 terms, two of which have coefficient 1/2), while the other terms give a nasty-looking sum: " !# ∞ π a0 cos(0kπ/N) + cos(Nkπ/N) N−1 N + ∑ ak + ∑ cos(nkπ/N) IN = N 2 2 n=1 k=1 =I+ π ∞ ∑ ak Sk , N k=1 where the a0 term is just the exact integral I and Sk is the sum:

N−1   kπ cos 0 kπ kπ N + cos N N Sk = + ∑ cos n . 2 N n=1 Amazingly (you will prove this for homework, below), this sum Sk is almost always zero. The only time it is nonzero is when k = 2mN where m is any integer (i.e., when k is an even multiple of N), in which case:   2mNπ cos (0) + cos (2mNπ) N−1 + ∑ cos n S2mN = 2 N n=1 = 1 + 1 N−1 + ∑ 1 = N, 2 n=1 using the fact that cos(2π`) = 1 for any integer `. This N factor cancels with π/N in IN Therefore, the nasty IN sum reduces to the remarkably simple expression: ∞ IN = I + π ∑ a2mN . m=1 2.3.3 Convergence rate The error in the trapezoidal rule is therefore exactly: ∞ IN − I = π ∑ a2mN . m=1 As N increases, this depends on ak for larger and larger k, but ak → 0 as k → ∞, and so the error goes to zero. The convergence rate of the Fourier cosine series determines the convergence rate of the trapezoidal rule! The faster ak → 0 with increasing k, the faster we will have IN − I → 0 with N. Fortunately, there is a long history of analysis of the convergence rate of the Fourier series which we can now bring to bear on the trapezoidal rule, and very general theorems are known. In the attached handout on cosine series, I go through one of the simplest analysis of convergence rates, which only requires repeated integration by parts on the a` formula. 6 In particular, for an arbitrary continuous f (x) that has no particularly special behavior at the endpoints x = 0, π, the coefficients ak approach an ak ≈ #/k2 dependence for large k, with some coefficient #. This means that, for large N, the trapezoidal-rule error approaches: π# ∞ 1 1 # = ∼ 2, ∑ 2 2 2 4N m=1 m N m=1 (2mN) ∞ IN − I ≈ π ∑ precisely the 1/N 2 dependence we predicted from a crude analysis above. (In fact, the sum of 1/m2 , the famous “Basel problem” solved by Euler, is exactly π 2 /6, so the error approaches π 3 #/24N 2 .) However, if f 0 vanishes at the endpoints, then ak ∼ 1/k4 and hence IN − I ∼ 1/N 4 . And if f 000 vanishes at the endpoints, then ak ∼ 1/k6 and IN − I ∼ 1/N 6 . In fact, we don’t actually care what ak does for odd k, only for even k, and for even k it turns out to be enough for f 0 etcetera to be periodic, i.e. f 0 (0) = f 0 (π) and so on, to obtain this faster convergence. And if all of the odd derivatives of f are periodic at the endpoints, e.g. if f (x) itself is periodic (and smooth, i.e. infinitely differentiable) or if f (x) is even around both boundaries (which makes all the odd derivatives vanish), then IN − I vanishes faster than any power of 1/N (typically exponentially fast).3 3 Clenshaw–Curtis quadrature Even if f (x) is a nice, smooth function inside the integration interval, it is often not even or periodic at its endpoints, a fact which is responsible for reducing the error of the trapezoidal rule to the ∼ 1/N 2 worst case. However, there is a simple change of variables which allows us to obtain the miraculous exponential convergence all the time (at least for smooth functions), and this technique is called Clenshaw–Curtis quadrature. For this section, it is more convenient to start with an integral over [−1, 1] rather than [0, π], and then make a change of variables x = cos θ to transform the former into the latter: Z Z 1 I= π f (x)dx = −1 f (cos θ ) sin θ dθ . 0 By construction, f (cos θ ) = f (cos[−θ ]) and f (cos[π −θ ]) = f (cos[π +θ ]) so f (cos θ ) is an even function around both x = 0 and x = π, and hence all of its odd-order derivatives with θ vanish at the endpoints. This is precisely one of the conditions under which the Fourier cosine series converges super-fast, and hence the trapezoidal rule converges super-fast for f (cos θ ). Unfortunately, this is spoiled in I thanks to our multiplication with sin θ , which is odd. We fix this problem by an additional trick: we first expand f (cos θ ) in a cosine series, and then do the I integrals analytically term-by-term: f (cos θ ) = ∞ a0 + ∑ ak cos(kθ ), 2 k=1 3 A more detailed analysis of exactly how fast the convergence is for periodic/even functions requires a fair amount of complex analysis (contour integration and so on). For example, if f (x) is an “analytic” function in a finite-width strip around the real-x axis, then one can show by contour integration that the error converges at least exponentially fast. See the references at the end. 7 ak = Z π I= 0 " 2 π Z π f (cos θ ) cos(kθ ) dθ , 0 # ∞ ∞ a0 2a2k . + ∑ ak cos(kθ ) sin θ dθ = a0 + ∑ 2 2 k=1 k=1 1 − (2k) Now, we have performed the I integral analytically, but how do we perform the integrals to obtain the Fourier coeffients ak ? We do these numerically by the trapezoidal rule! That is: "  # 2 π f (cos 0) + f (cos π) N−1  knπ nπ  ak ≈ × cos cos . + ∑ f cos π N 2 N N n=1 But how is this an improvement? We have just swapped one trapezoidal rule summation for a zillion of them! Three points save us: • The a` integrand, f (cos θ ) cos(`θ ), is even around its endpoints and so the trapezoidal rule for ak converges super-fast (usually exponentially fast). So, we only need a small N. • f (cos θ ) is even around its endpoints, so the coefficents ak go to zero super-fast (usually exponentially fast). So, we only need a small number of ak coeffients. In fact, the convergence rate for ak is about the same as the convergence rate of the trapezoidal rule, so we usually only compute N coefficients ak . • When computing ak , we always evaluate f (x) at the same points f (cos nπ N ) independent of k. So, we only need to evaluate f (x) at N + 1 points once. There are actually additional nice properties. The trapezoidal rule for ak turns out to be a very special summation called a discrete cosine transform, allowing us to factorize the sum in a way that shares computations for different k via a “fast Fourier transform” (FFT). All in all, it turns out that one can accurately compute the integral in ∼ N log N operations. The key idea was that, by changing variables x = cos θ and doing a cosine-series expansion of f (cos θ ), we can transform any integral into integrals purely of even functions, allowing us to apply the miraculous convergence of the trapezoidal rule for such functions to any integral. (The only thing that will spoil this is if f has other discontinuities or singularities within the integration interval.) 3.1 Chebyshev polynomials I would be remiss not to mention a connection between the above analysis and another beautiful branch of mathematics: Chebyshev approximation. Look back at the cosine series for f (cos θ ) but write it in terms of the original variable x = cos θ via θ = cos−1 x: ∞ ∞ a0 a0 f (x) = + ∑ ak cos(k cos−1 x) = + ∑ ak Tk (x), 2 k=1 2 k=1 8 which is expanding f (x) in the nasty-looking functions Tk (x) = cos(k cos−1 x). What a horrible way to disguise a cosine series! But then a miracle occurs: the functions Tk (x) are not nasty at all, they are actually just polynomials called Chebyshev polynomials: T0 (x) = cos(0 cos−1 x) = 1, T1 = cos(cos−1 x) = x, T2 = cos(2 cos−1 x) = 2 cos2 (cos−1 x) − 1 = 2x2 − 1, Tk+1 (x) = 2xTk (x) − Tk−1 (x), where the T2 and Tk+1 expressions use angle-addition identities. So, expanding f (cos θ ) in a cosine series is exactly equivalent to expanding f (x) in a series of polynomials. But they are special, wonderful, beautiful polynomials that have all the amazing properties of a cosine series. The coefficients ak of these polynomials converge exponentially fast for smooth f (x). These polynomials Tk (x) stay in between −1 and 1 for x in [−1, 1], and all of their k roots are in the interval [−1, 1]. We can compute the coefficients ak quickly by FFT algorithms. And, of course, once we have expanded f (x) in polynomials, almost anything that we might want to do (integration, differentiation, finding roots, etcetera) becomes easy. 4 Further reading Another review of the convergence analysis for the trapezoidal rule via the Fourier series using complex exponentials rather than cosine series, also handling the case of discontinuities within the integration interval, can be found online at S. G. Johnson, “Notes on the convergence of trapezoidal-rule quadrature,” http://math.mit. edu/~stevenj/trapezoidal.pdf. Much more information about the convergence rates of Fourier series and related numerical methods can be found e.g. in the book Chebyshev and Fourier Spectral Methods by John P. Boyd (available online at http: //tinyurl.com/24mepgc). A wonderful practical realization of Clenshaw–Curtis quadrature and a host of related numerical techniques can be found in the chebfun package by Trefethen et al., which available as free software for Matlab along with a series of very readable review-style publications at http://www2.maths.ox.ac.uk/ chebfun/. Many numerical-analysis textbooks instead analyze the error of the trapezoidal rule via something called Bernoulli polynomials; this approach can be found (with references) in the Wikipedia article on the Euler–Maclaurin formula or in (e.g.) An Introduction to Numerical Analysis by Süli and Mayers or An Introduction to Numerical Analysis by Stoer and Bulirsch. However, I find it much nicer to relate the problem to Fourier series, which are not only much more widely known and more widely applicable, but also point the way towards more accurate methods such as Clenshaw–Curtis. I’ve also put some notes and links on Clenshaw–Curtis techniques up on Wikipedia (“Clenshaw–Curtis quadrature,” Wikipedia, accurate as of 5 January 2011). 9