Abstract

The error in the trapezoidal rule quadrature formula can be attributed to discretization in the interior and non-periodicity at the boundary. Using a contour integral, we derive a unified bound for the combined error from both sources for analytic integrands. The bound gives the Euler–Maclaurin formula in one limit and the geometric convergence of the trapezoidal rule for periodic analytic functions in another. Similar results are also given for the midpoint rule.

1. Introduction

Let f be continuous on [0,1] and let n be a positive integer. The (composite) trapezoidal rule approximates the integral
1.1by the sum
1.2where the prime indicates that the terms k=0 and k=n are multiplied by 1/2. Throughout this paper, f may be real or complex, and ‘periodic’ means periodic with period 1.

The approximation of I by In has many interesting properties. One is that, if f is periodic and analytic, the convergence is geometric. This observation in some sense goes back to Poisson in the 1820s [1], though it seems to have been Davis in 1959 who first stated a theorem [2,3]. Another is that, for non-periodic f, the accuracy is O(n−2) and this can be improved to O(n−4), O(n−6) and so on by subtracting appropriate multiples of f′(1)−f′(0), f′′′(1)−f′′′(0) and so on, if f is sufficiently smooth. The latter process is described by the Euler–Maclaurin formula, published independently by Euler and Maclaurin around 1740 [4,5]. Compendia of results related to the Euler–Maclaurin formula can be found in [6,7].

A standard derivation of Davis's result involves contour integrals in the complex plane, and contour integrals can also be used to derive the Bernoulli numbers that appear in the Euler–Maclaurin formula. With these facts in mind, we have attempted to develop a unified formulation based on contour integrals that would make it possible to derive both kinds of results at once for analytic integrands. It is hard to believe that our theorem can be new, but we have been unable to find such a result in the literature. The closest we have found is [8], appendix B, which sets up the problem in the same way without deriving an explicit estimate. Another related reference is [9], §8.3, which gives more mathematical detail than [8] but considers the special case in which f is analytic in the infinite strip 0<Re z<1, , the context of the Abel–Plana formula, which can also be found discussed in earlier references such as [10], §3.14. All in all, our impression is that, whereas the techniques we apply in this paper are old ones, they may never have been combined before to derive an explicit error estimate for a function analytic in a finite strip.

The theorem is stated in §2. Various existing results are derived as corollaries in §3, and the proof of the theorem is presented in §4. Section 5 mentions the variation of the midpoint rather than trapezoidal formula, and the discussion in §6 points to connections to rational approximation and the theory of hyperfunctions.

2. Theorem

The theorem is stated in terms of the following Euler–Maclaurin correction sum. For any m≥0 and sufficiently smooth f, we define
2.1where Bk is the kth Bernoulli number (). We further define
2.2and
2.3Note that Qm,n, Δm and Δ(m+1) all depend on f, though the notation does not make this explicit.

Here is the theorem. The region of analyticity is sketched in figure 1 (see also §4).

Contour of integration for the proof of theorem 2.1. At z=0 and 1, the integral is defined in the principal value sense, as suggested by the gaps in the contour. The dots show the sample points of the trapezoidal rule, which become the poles of the characteristic function . The theorem bounds the error by the sum of three terms: Einterior corresponding to ΓT and ΓB, Eboundary corresponding to ΓL and ΓR, and Etail associated with extending the integrals along ΓL and ΓR to infinite intervals.

Theorem 2.1

Given real numbers a>0 and M≥0 and an integer m≥0, let f satisfy | f(z)|≤M and have a continuous (m+1)st derivative in the region defined by 0≤Re z≤1, −a<Im z<a, and be analytic in the interior of this region. Then
2.4with2.52.6and2.7

In words, theorem 2.1 breaks the error of the trapezoidal rule adjusted by Qm,n into three terms, the first is related to the discretization error in the interior and the others to boundary effects. Two of these are exponentially small as , and the other is algebraically small. We use the labels ‘boundary’ and ‘tail’ for reasons that will become evident in the proof.

Although theorem 2.1 is valid for any m≥0, one would not normally apply it for odd values of m, as Qm,n=Qm+1,n and Δm=Δm+1 when m is odd. This means that if m is odd, then increasing it to the next even number yields the same bound except with O(n−m−2Δ(m+1)) improved to O(n−m−3Δ(m+2)), assuming f(m+2) exists and is continuous.

3. Corollaries

By considering special cases of theorem 2.1, we obtain various familiar results. The first is Davis's theorem for periodic integrands; see [[2];[3], §4].

Corollary 3.1

Proof.

By (2.1)–(2.3), Qm,n, Δm and Δ(m+1) are zero when f is periodic. It follows from (2.6) and (2.7) that Eboundary and Etail are zero in this case too. The bound (3.1) now follows from (2.4) and (2.5). ▪

The second corollary is one version of the Euler–Maclaurin formula.

Corollary 3.2

Let f be analytic on [0,1]. Then for any m≥0,
3.2as.

Proof.

If f is analytic on [0,1], it is analytic and bounded in the strip around [0,1] of half-width a for some a>0. The result now follows from (2.4) to (2.7) because, as , Eboundary=O(n−m−2) and both Einterior and Etail are of asymptotically smaller order, O(nme−2πan). ▪

If f is a polynomial of degree at most m+1, the mth Euler–Maclaurin approximation is exact.

Corollary 3.3

Let f be a polynomial of degree at most m+1 for somem≥0. Then for any n,
3.3

Proof.

If f is a polynomial of degree at most m+1, then f(m+1) is a constant, implying Δ(m+1)=0, so (2.6) implies that Eboundary=0 in (2.4). The bounds for the other terms Einterior and Etail contain the factor e−2πan, where n is fixed but a can be taken as large as we want as a polynomial is an entire function. In the case of Etail, Δm is fixed, so (2.7) implies that Etail becomes arbitrarily small as . In the case of Einterior, M grows as , but only at a polynomial rate, so (2.5) implies that this term too becomes arbitrarily small as . Thus, In−I−Qm,n must be equal to zero. ▪

Corollary 3.3 leads to the identity known as Faulhaber's formula.

Corollary 3.4

For any n≥1 and m≥1,
3.4

Proof.

This is equation (3.3) in the special case f(x)=nm+1xm; the term nm+1/(m+1) is the integral I, and the term nm/2 appears because In is defined in (1.2) with a factor 1/2 multiplying the term k=n. The reason for the assumption m≥1 is that, in the (trivial) case m=0, the sum on the left is missing a non-zero contribution 1/2 corresponding to the k=0 term in the trapezoidal sum. ▪

4. Proof

We now prove theorem 2.1.

As sketched in figure 1, let Γ be the boundary of the rectangle of analyticity of f, oriented in the positive sense, and let ΓL, ΓR, ΓT and ΓB be the left, right, top and bottom boundary segments, respectively. In the following argument, we suppose for simplicity that f extends continuously to ΓT and ΓB. If it does not, then the required result can be obtained by replacing ΓT by ΓT−εi and ΓB by ΓB+εi and considering ε→0.

The ‘characteristic function’ has simple poles at z=k/n for each integer k with residue (2πin)−1. It follows from residue calculus that trapezoidal approximation (1.2) can be represented by the contour integral
4.1where the integrals over ΓL and ΓR are taken in the principal value sense so as to introduce the necessary factors of 1/2.

The true integral I can also be represented by a contour integral over Γ4.2where u is defined by
4.3We can derive this formula by noting that I can be regarded as half the integral of f from 0−0i to 1−0i minus half the integral of f from 1+0i to 0+0i. As f is analytic in the rectangular region and u is analytic in the upper and lower half-planes, these two contours of integration can be deformed to the upper and lower halves of Γ without changing the value of the integral.

Combining (4.1)–(4.3), we find
4.4with , which simplifies to
4.5We now establish (2.4) by breaking (4.4) into four terms,
4.6corresponding to the integrals of f(z)S(z) over the four segments of the boundary (figure 1). We note immediately that, by (4.4) and (4.5), IT and IB are each bounded by M/(e2πan−1), which implies bound (2.5) with Einterior defined by
4.7To complete the derivation of (2.4)–(2.7), by (2.4) and (4.6) and (4.7), we must show that IL+IR can be broken into the pieces
4.8with Eboundary and Etail satisfying (2.6) and (2.7).

For y∈(−a,a), define
4.9This definition simplifies (2.2) and (2.3) to
4.10and
4.11As S(1+iy)=S(iy), our task is to estimate
4.12where the integral is taken in the principal value sense. As g has a continuous (m+1)st derivative on (−a,a), one of the standard forms of Taylor's theorem with remainder gives
4.13for y∈(−a,a) [11], theorem 1.36. In (4.12), we note that S(iy) is an odd function of y. Therefore, when the sum on the left-hand side of (4.13) is inserted in (4.12) as an approximation to g(y), the contributions from even values of k vanish (in the case k=0, we use the fact that it is a principal value integral). The result is
4.14because for y>0.

We can now identify the pieces Qm,n, Eboundary and Etail. The quantity Eboundary is the number inside absolute value signs on the left of (4.14), satisfying
4.15The quantity Qm,n is the negative of the sum on the left in (4.14), but with the integrals running from 0 to instead of a4.16And Etail is the error just introduced in extending those integrals to 4.17These definitions ensure that (4.8) holds as required; what remains is to derive (2.6) from (4.15) and (2.7) from (4.17) and to confirm that (4.16) matches the definition of Qm,n given in (2.1).

That (4.16) matches (2.1) follows from an identity that goes back to the late nineteenth century [12,13],
4.18as g(k)(0)=ik( f(k)(1)−f(k)(0)) and i2k+2=1 for k odd.

To derive (2.6) from (4.15), we note that the change of variables t=2πny and extension of the upper limit of integration to in (4.15) gives
4.19Bound (2.6) follows from this together with an identity closely related to (4.18) [7], (25.5.1)
4.20where ζ is the Riemann zeta function. The numbers ζ(m+2) decrease monotonically from their maximum ζ(2)=π2/6 in the case m=0.

To derive (2.7) from (4.17), we note that, by (4.10) and (4.17),
4.21or, after the change of variables t=2πny,
4.22Applying lemma A.1 of appendix A with b=2πna and using the inequalities (2πn)−k−1≤(2π)−2 and (b+1)k≤(b+1)m, we obtain (2.7).

5. Midpoint rule variant

A close cousin of trapezoidal rule (1.2) is the midpoint rule5.1where now no prime is needed in the sum as all the terms have the same weight. All the arguments of this paper can be carried through for this case with minor changes. The term in (4.1) becomes , and the correction sum Qm,n of (2.1) is adjusted slightly to become
5.2The constant −1 becomes +1 in the denominators of equations (4.14)–(4.17) and (4.19)–(4.22), altering Eboundary and Etail in a manner that leaves them still bounded as before, though slightly tighter bounds could now be derived. In the case of (4.18), 1 becomes −1 in the denominator on the left, and the right-hand side is consequently multiplied by the factor 2−k−1, explaining the appearance of this factor in the definition of above.

However, it is not necessary to carry out the estimates just described, because the n-point midpoint rule can be analysed in an elementary way as a linear combination of twice the 2n-point trapezoidal rule minus the n-point trapezoidal rule
5.3In view of the factor nk+1 in the denominator of (2.1), this explains instantly the appearance of the factor 2−k−1 in (5.2).

Here is the analogue of theorem 2.1, and corollaries 3.1–3.3, for the midpoint rule. As we have mentioned, in this case the bounds on Eboundary and Etail could be further sharpened somewhat. For a statement of the midpoint rule variant of the Euler–Maclaurin formula in the more standard form of an asymptotic series, see [6], eqn (2.9.25), and for a more detailed treatment, see [14].

Theorem 5.1

Let f, a, M and m be as in theorem 2.1, and letandbe defined by (5.1) and (5.2). Then
5.4with Einterior, Eboundaryand Etailsatisfying the same bounds (2.5)–(2.7) as before. Corollaries 3.1–3.3 hold as before with Inand Qm,nreplaced byand.

In analogy to corollary 3.4, theorem 5.1 implies the following Faulhaber-like formula for sums of odd powers of integers. This can be derived by applying the analogue of corollary 3.3, , to the function f(x)=nm+1xm. Alternatively and more simply, it follows by combining corollary 3.4 and (5.3), that is, regarding a sum of powers of odd integers as a sum of powers of all integers minus a sum of powers of even integers. Again the key difference is the appearance of the factor 2−k−1.

Corollary 5.2

For any n≥1 and m≥1,
5.5

6. Discussion

In the theory of hyperfunctions, delta functions and other distributions are realized not by the test functions and linear functionals of real analysis, but by methods of complex analysis. Specifically, a hyperfunction on a real interval is defined as a difference of analytic functions in the upper and lower half-planes, or, more precisely, an equivalence class of such differences [15–17]. Our arguments have exactly this flavour, and, in particular, the function S(z) is expressed in (4.5) in hyperfunction form. The reason hyperfunction theory is relevant is that it provides a convenient framework in which to compare the integral I with the trapezoidal or midpoint approximations In, which are regarded essentially as integrals whose integrands are strings of delta functions. This paper is a contribution towards a longer term goal of strengthening the links between hyperfunction theory and numerical analysis.

Going beyond the trapezoidal and midpoint rules, it may be noted that, whenever an integral I of an analytic function f is approximated by a quadrature formula defined by nodes {xk} and weights {wk}, In can be written as a contour integral involving the product r(z) f(z), where r is the type (n,n+1) rational function with poles xk and residues wk. Writing I itself as a contour integral of f times a hyperfunction, such as u(z) in (4.3), makes it possible to estimate In−I by contour integrals. This technique was pioneered by Takahasi & Mori [18], who had the vision of connecting numerical analysis and hyperfunctions long before we did, and it was applied to the comparison of Gauss and Clenshaw–Curtis quadrature formulae in [19]. Such analyses highlight the fact that every quadrature formula implicitly makes use of a rational approximation, and the properties of these rational approximations are investigated in [[3], §14; [20],[21]].

Funding statement

This work was supported by the European Research Council under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement no. 291068. The views expressed in this article are not those of the ERC or the European Commission, and the European Union is not liable for any use that may be made of the information contained here.

Acknowledgements

We have benefited by the good advice from Endre Süli and André Weideman.

Appendix A

The following inequality was used in the proof of §4.

Lemma .1

For any real number b≥0 and integer k≥1,
1

Proof.

As observed below (4.20), the left-hand side of (.1) is ≤π2/6≈1.64, whereas it can be verified that the right-hand side is greater than this value for b≤0.75. To complete the proof, we may accordingly assume b>0.75. As e0.75>2, replacing the denominator in (A 1) by et decreases the integral by less than a factor of 2, so it is enough to show
2for b>0.75. We can do this by induction in k. For k=1, the inequality holds as an equality, as can be verified by integration by parts. Assume then that it holds for some k≥1 and consider the case k+1. Integration by parts gives
by the inductive hypothesis. Cancelling the common factor of e−b, this leaves us with the problem of establishing
which follows because the left-hand side is less than bk+1+(b+1)k and the right-hand side is equal to b(b+1)k+(b+1)k. ▪

References

. On the numerical integration of periodic analytic functions. In On numerical integration: Proceedings of a Symposium, Madison, 21–23 April 1958 (ed. LangerRE), pp. 45–59. Mathematics Research Center.Wisconsin, WI: University of Wisconsin.