## 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 *I*_{n} 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 *B*_{k} is the *k*th Bernoulli number (). We further define
2.2and
2.3Note that *Q*_{m,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).

### 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.4*with*
2.5
2.6
*and*
2.7

In words, theorem 2.1 breaks the error of the trapezoidal rule adjusted by *Q*_{m,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 *Q*_{m,n}=*Q*_{m+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

*Let f be analytic and* 1-

*periodic with*|

*f*(

*z*)|≤

*M*

*in the region*0≤Re

*z*≤1, −

*a*<Im

*z*<

*a*

*for some*

*a*>0.

*Then*, 3.1

### Proof.

By (2.1)–(2.3), *Q*_{m,n}, Δ_{m} and Δ^{(m+1)} are zero when *f* is periodic. It follows from (2.6) and (2.7) that *E*_{boundary} and *E*_{tail} 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.2

*as*.

### 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 , *E*_{boundary}=*O*(*n*^{−m−2}) and both *E*_{interior} and *E*_{tail} are of asymptotically smaller order, *O*(*n*^{m} *e*^{−2πan}). ▪

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

### Corollary 3.3

*Let f be a polynomial of degree at most m*+1

*for some*

*m*≥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 *E*_{boundary}=0 in (2.4). The bounds for the other terms *E*_{interior} and *E*_{tail} 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 *E*_{tail}, Δ_{m} is fixed, so (2.7) implies that *E*_{tail} becomes arbitrarily small as . In the case of *E*_{interior}, *M* grows as , but only at a polynomial rate, so (2.5) implies that this term too becomes arbitrarily small as . Thus, *I*_{n}−*I*−*Q*_{m,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*)=*n*^{m+1}*x*^{m}; the term *n*^{m+1}/(*m*+1) is the integral *I*, and the term *n*^{m}/2 appears because *I*_{n} 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−0*i* to 1−0*i* minus half the integral of *f* from 1+0*i* to 0+0*i*. 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), *I*_{T} and *I*_{B} are each bounded by *M*/(*e*^{2πan}−1), which implies bound (2.5) with *E*_{interior} 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 *I*_{L}+*I*_{R} can be broken into the pieces
4.8with *E*_{boundary} and *E*_{tail} 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 *Q*_{m,n}, *E*_{boundary} and *E*_{tail}. The quantity *E*_{boundary} is the number inside absolute value signs on the left of (4.14), satisfying
4.15The quantity *Q*_{m,n} is the negative of the sum on the left in (4.14), but with the integrals running from 0 to instead of *a*
4.16And *E*_{tail} 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 *Q*_{m,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)=*i*^{k}( *f*^{(k)}(1)−*f*^{(k)}(0)) and *i*^{2k+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 rule*
5.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 *Q*_{m,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 *E*_{boundary} and *E*_{tail} 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 2*n*-point trapezoidal rule minus the *n*-point trapezoidal rule
5.3In view of the factor *n*^{k+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 *E*_{boundary} and *E*_{tail} 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 let* *and* *be defined by (*5.1*) and (*5.2*). Then
*5.4*with E*_{interior}*, E*_{boundary} *and E*_{tail} *satisfying the same bounds (*2.5*)–(*2.7*) as before. Corollaries 3.1–3.3 hold as before with I*_{n} *and Q*_{m,n} *replaced by* *and* .

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*)=*n*^{m+1}*x*^{m}. 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 *I*_{n}, 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 {*x*_{k}} and weights {*w*_{k}}, *I*_{n} 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 *x*_{k} and residues *w*_{k}. Writing *I* itself as a contour integral of *f* times a hyperfunction, such as *u*(*z*) in (4.3), makes it possible to estimate *I*_{n}−*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 *e*^{0.75}>2, replacing the denominator in (A 1) by e^{t} 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 *b*^{k+1}+(*b*+1)^{k} and the right-hand side is equal to *b*(*b*+1)^{k}+(*b*+1)^{k}. ▪

- Received August 27, 2013.
- Accepted October 15, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.