43.09.02 · numerical-analysis / 09-numerical-quadrature

Composite rules, Euler-Maclaurin, and Romberg / adaptive quadrature

shipped3 tiersLean: none

Anchor (Master): Davis-Rabinowitz 1984 *Methods of Numerical Integration* 2e (Academic Press) Ch. 2-4, 6 (composite rules, Euler-Maclaurin, Romberg extrapolation and the periodic/trapezoidal spectral accuracy); Trefethen-Weideman 2014 *The exponentially convergent trapezoidal rule* (SIAM Review 56) (geometric convergence for analytic periodic integrands); Henrici 1977 *Applied and Computational Complex Analysis* vol. 2 (Wiley) Ch. 11 (Euler-Maclaurin, Bernoulli numbers, and the trapezoidal rule on the circle)

Intuition Beginner

A single trapezium or parabola laid across a whole interval is a crude stand-in for a curve that bends a lot. The fix is the one you would reach for instinctively: chop the interval into many narrow strips, lay a small trapezium over each strip, and add up the little areas. Over a narrow strip even a wiggly curve looks almost straight, so each small trapezium is accurate, and the total is far better than one giant trapezium. This chopping-and-summing recipe is called a composite rule, and it is the workhorse of practical integration.

Halving the width of the strips does something predictable to the error. For the composite trapezium rule, cutting the strip width in half cuts the error to about a quarter. That fixed factor of four is a gift, because it lets you guess the answer you would get with infinitely many strips before you ever compute it. If two estimates differ by a known pattern, you can combine them to cancel the leading mistake and leap to a much sharper value. Doing this combining over and over is the idea behind a method called Romberg integration.

There is a deeper reason the error shrinks by a clean factor of four, and it comes from a remarkable bookkeeping formula. The formula says the trapezium error is built entirely out of corrections measured at the two endpoints of the interval — how steep the curve is there, how its steepness is changing, and so on. Nothing from the middle of the interval contributes to the leading error. This has a startling payoff: if the curve repeats itself like a wave, so that its behaviour at the right endpoint exactly matches the left, all those endpoint corrections cancel, and the humble trapezium rule becomes spectacularly accurate.

Finally, real curves are not uniformly difficult. A function may be calm over most of its range and sharply peaked in one small region. Spacing strips evenly wastes effort on the calm parts and starves the peak. Adaptive quadrature watches the local error as it goes, and spends strips where the curve misbehaves and few where it is tame, refining only the regions that need it.

Visual Beginner

Picture the interval cut into several equal strips, with a small straight chord laid across the top of each strip. The composite trapezium estimate is the sum of all the little trapezoid areas. Doubling the number of strips halves their width and roughly quarters the total error, and the table below tracks that pattern. The second idea in the picture is the endpoint bookkeeping: the leftover error is read off only from the two ends of the whole interval, not from anything in between.

number of strips	strip width	rough error size
$N$	$h$	$E$
$2 N$	$h \div 2$	about $E \div 4$
$4 N$	$h \div 4$	about $E \div 16$

Each row halves the width and quarters the error, so the error behaves like the square of the strip width. That fixed factor of four is exactly the lever that Romberg integration pulls: knowing the error shrinks by four lets you cancel most of it by combining two estimates.

Worked example Beginner

Let us estimate the area under $f (x) = 1 \div x$ from $1$ to $2$ with the composite trapezium rule using two strips, then refine to four and watch the error pattern. The exact area is the natural logarithm of $2$ , about $0.6931$ .

Step 1. Two strips. The width is $1$ , split into two strips of width $0.5$ , with sample points $1$ , $1.5$ , $2$ . The values are $f (1) = 1$ , $f (1.5) = 0.6667$ , $f (2) = 0.5$ .

Step 2. Apply the composite trapezium rule. Count the two end values at half weight and the interior value at full weight, then multiply by the strip width: $0.5 \times (0.5 \times 1 + 0.6667 + 0.5 \times 0.5) = 0.5 \times (0.5 + 0.6667 + 0.25) = 0.5 \times 1.4167 = 0.7083$ .

Step 3. The two-strip error. The estimate $0.7083$ minus the true $0.6931$ is about $0.0152$ .

Step 4. Four strips. Width $0.25$ , sample points $1, 1.25, 1.5, 1.75, 2$ with values $1, 0.8, 0.6667, 0.5714, 0.5$ . The rule gives $0.25 \times (0.5 \times 1 + 0.8 + 0.6667 + 0.5714 + 0.5 \times 0.5) = 0.25 \times 2.7881 = 0.6970$ .

Step 5. The four-strip error. The estimate $0.6970$ minus $0.6931$ is about $0.0039$ .

What this tells us: halving the strip width took the error from $0.0152$ down to $0.0039$ , a drop by a factor close to four, exactly as the square-of-the-width rule predicts. That predictable factor is what makes the next step — combining the two estimates to cancel the error — possible.

Check your understanding Beginner

Exercise (easy, true-false).

True or false: for a smooth wave-like function that repeats so its left end matches its right end, the trapezium rule can be far more accurate than the square-of-the-width rule alone would suggest.

Hint

Recall that the leftover error is read off only from the two endpoints, and ask what happens when those endpoints match.

Answer

True. The trapezium error is built from endpoint corrections, and when the function repeats, the corrections at the two ends cancel, so the error collapses far below the usual size. Feedback-correct: this endpoint cancellation is why the trapezium rule is the method of choice for smooth periodic integrands. Feedback-wrong: if you answered false, recall that nothing from the interior contributes to the leading error, so matching endpoints removes it.

Formal definition Intermediate+

Throughout, $f : [a, b] \to R$ , the interval is partitioned into $N$ equal subintervals of width $h = (b - a) / N$ with nodes $x_{j} = a + j h$ for $j = 0, \dots, N$ , and $f_{j} = f (x_{j})$ . The space of polynomials of degree at most $k$ is $Π_{k}$ , and $f \in C^{m} [a, b]$ means $f$ has continuous derivatives through order $m$ . The single-interval trapezium and Simpson rules and their Peano-kernel error terms are those of 43.09.01.

Definition (composite trapezium and Simpson rules). The composite trapezium rule on the partition is $$ T_N(f) = h\Big(\tfrac12 f_0 + f_1 + \dots + f_{N-1} + \tfrac12 f_N\Big) = h\sum_{j=0}^{N}{}' f_j , $$ where the primed sum halves the two endpoint terms. The composite Simpson rule (for even $N$ ) groups the subintervals in pairs: $$ S_N(f) = \frac{h}{3}\Big(f_0 + 4f_1 + 2f_2 + 4f_3 + \dots + 2f_{N-2} + 4f_{N-1} + f_N\Big). $$ Both are obtained by applying the corresponding single-interval rule of 43.09.01 on each subinterval (each pair, for Simpson) and adding. All weights are positive, so neither rule amplifies data error.

Definition (Bernoulli numbers and polynomials). The Bernoulli polynomials $B_{k} (x)$ are defined by the generating function $\frac{t e ^{x t}}{e ^{t} - 1} = \sum_{k = 0}^{\infty} B_{k} (x) \frac{t ^{k}}{k !}$ , and the Bernoulli numbers are $B_{k} = B_{k} (0)$ . They satisfy $B_{k}^{'} (x) = k B_{k - 1} (x)$ , $\int_{0}^{1} B_{k} (x) d x = 0$ for $k \geq 1$ , and $B_{2 k + 1} = 0$ for $k \geq 1$ ; the first even values are $B_{2} = \frac{1}{6}$ , $B_{4} = - \frac{1}{30}$ , $B_{6} = \frac{1}{42}$ . The periodic Bernoulli function $B_{k} (t) = B_{k} ({t})$ , with ${t}$ the fractional part, extends $B_{k}$ on $[0, 1)$ to a $1$ -periodic function.

Definition (Romberg T-table). Let $T_{j}^{(0)} = T_{N_{j}} (f)$ be the composite trapezium estimate on $N_{j} = 2^{j}$ subintervals, $j = 0, 1, 2, \dots$ . The Romberg table is generated by the Richardson recurrence $$ T^{(k)}_j = T^{(k-1)}_j + \frac{T^{(k-1)}j - T^{(k-1)}{j-1}}{4^{k} - 1}, \qquad k \ge 1,\ j \ge k . $$ The column $k$ has leading error $O (h^{2 k + 2})$ ; the diagonal entries $T_{k}^{(k)}$ are the Romberg estimates.

The symbols $\sum$ , $\int$ , the Bernoulli polynomials $B_{k}$ , the periodic Bernoulli function $B_{k}$ , the fractional part ${t}$ , and the primed summation $\sum^{'}$ are recorded in _meta/NOTATION.md.

Counterexamples to common slips Intermediate+

"The composite trapezium error expansion converges as a power series in $h$ ." It is an asymptotic expansion: for fixed $f$ the series in $h^{2 k}$ is generally divergent, because the Bernoulli numbers grow like $(2 k)!$ . One truncates at an optimal index and keeps the remainder term; treating it as a convergent series is wrong.
"Euler-Maclaurin makes every integrand integrable to spectral accuracy by the trapezium rule." The endpoint corrections vanish only for smooth periodic integrands (or those with matching endpoint derivatives). For a generic $f$ on $[a, b]$ the corrections $f^{(2 k - 1)} (b) - f^{(2 k - 1)} (a)$ are nonzero and the rule remains $O (h^{2})$ .
"Romberg's extrapolation factor is always $4^{k} - 1$ in the denominator." That factor is correct only because the trapezium error expansion runs in even powers $h^{2}, h^{4}, \dots$ . A rule whose error has an $h^{1}$ or odd-power term needs the factor $2^{p_{k}} - 1$ keyed to its actual error orders; applying the $4^{k}$ pattern to such a rule corrupts the result.
"Adaptive quadrature with a local tolerance $ε$ guarantees a global error below $ε$ ." The local estimates can be unreliable where $f$ is rough, and summing many local errors near tolerance can exceed the intended global bound; robust schemes apportion the global tolerance across subintervals and guard against under-resolved spikes.

Key theorem with proof Intermediate+

The signature result is the Euler-Maclaurin summation formula. It writes the difference between the composite trapezium sum and the exact integral as a finite sum of endpoint-derivative corrections weighted by Bernoulli numbers, plus an integral remainder. The expansion runs in even powers of $h$ , which both explains the clean $O (h^{2})$ leading behaviour and powers the Richardson extrapolation underlying Romberg integration. The proof is repeated integration by parts against the Bernoulli polynomials on each subinterval. This follows Süli-Mayers ^{[Süli-Mayers §7.6]}.

Theorem (Euler-Maclaurin summation formula). Let $f \in C^{2 m + 2} [a, b]$ , let $h = (b - a) / N$ , and let $T_{N} (f)$ be the composite trapezium rule. Then $$ T_N(f) - \int_a^b f(x),dx = \sum_{k=1}^{m} \frac{B_{2k}}{(2k)!},h^{2k}\Big[f^{(2k-1)}(b) - f^{(2k-1)}(a)\Big] + R_m, $$ where the remainder is $$ R_m = -\frac{h^{2m+2}}{(2m+2)!}\int_a^b \widetilde{B}{2m+2}!\Big(\tfrac{x-a}{h}\Big), f^{(2m+2)}(x),dx , $$ *and $|R_m| \le \dfrac{|B{2m+2}|}{(2m+2)!},h^{2m+2}(b-a),\max_{[a,b]}|f^{(2m+2)}|$.*

Proof. It suffices to prove the formula on a single subinterval $[0, 1]$ (after the affine rescaling $x = x_{j} + h s$ , $s \in [0, 1]$ ) and sum. On $[0, 1]$ the claim is the Euler-Maclaurin one-panel identity $$ \frac{g(0)+g(1)}{2} - \int_0^1 g(s),ds = \sum_{k=1}^{m}\frac{B_{2k}}{(2k)!}\Big[g^{(2k-1)}(1) - g^{(2k-1)}(0)\Big] - \frac{1}{(2m+2)!}\int_0^1 \widetilde{B}_{2m+2}(s),g^{(2m+2)}(s),ds, $$ for $g \in C^{2 m + 2} [0, 1]$ . Begin from $\int_{0}^{1} g (s) d s$ and integrate by parts using $B_{1} (s) = s - \frac{1}{2}$ , so that $B_{1}^{'} (s) = 1$ : $$ \int_0^1 g,ds = \big[B_1(s)g(s)\big]_0^1 - \int_0^1 B_1(s),g'(s),ds = \frac{g(1)+g(0)}{2} - \int_0^1 B_1(s),g'(s),ds , $$ using $B_{1} (1) = \frac{1}{2}$ , $B_{1} (0) = - \frac{1}{2}$ . This already gives the trapezium value minus a correction integral.

Now integrate the correction repeatedly using $B_{k}^{'} (s) = k B_{k - 1} (s)$ , i.e. $B_{k - 1} (s) = \frac{1}{k} B_{k}^{'} (s)$ , so each integration by parts raises the Bernoulli index by one and a derivative onto $g$ . At each even index $2 k$ the boundary term contributes $\frac{B _{2 k}}{( 2 k )!} [g^{(2 k - 1)} (1) - g^{(2 k - 1)} (0)]$ , using $B_{2 k} (1) = B_{2 k} (0) = B_{2 k}$ ; at each odd index the boundary term vanishes because $B_{2 k + 1} (0) = B_{2 k + 1} (1) = 0$ for $k \geq 1$ . After $2 m + 2$ integrations by parts the surviving boundary terms are exactly the displayed Bernoulli sum, and the final integral is the periodic-Bernoulli remainder. Summing the one-panel identity over the $N$ subintervals telescopes the interior endpoint terms — the derivative correction at $x_{j}$ from the right of one panel cancels the one from the left of the next — leaving only the global endpoints $a$ and $b$ , and assembling the remainder integrals into one over $[a, b]$ gives the stated formula. The bound on $R_{m}$ follows from $\int_{0}^{1} ∣ B_{2 m + 2} (s) ∣ d s \leq ∣ B_{2 m + 2} ∣$ (the periodic Bernoulli function attains its extreme value $∣ B_{2 m + 2} ∣$ at the integers) applied panel-by-panel. $□$

Corollary (composite error orders and Romberg). Taking $m = 1$ recovers the global trapezium error $T_{N} (f) - \int_{a}^{b} f = \frac{h ^{2}}{12} [f^{'} (b) - f^{'} (a)] + O (h^{4})$ , hence $O (h^{2})$ , matching the summed single-interval Peano term of 43.09.01; the composite Simpson error is $O (h^{4})$ . Because the trapezium expansion contains only even powers $h^{2}, h^{4}, \dots$ , one Richardson step on $T_{N}$ and $T_{2 N}$ eliminates the $h^{2}$ term: $\frac{4 T _{2 N} - T _{N}}{3} = \int_{a}^{b} f + O (h^{4})$ , which is exactly composite Simpson, and iterating eliminates $h^{4}, h^{6}, \dots$ in turn — the Romberg table.

Bridge. The Euler-Maclaurin formula is the foundational reason the trapezium rule's error is not a single mean-value term but a structured series in even powers of $h$ with endpoint-derivative coefficients: the interior contributions telescope away, leaving only what happens at $a$ and $b$ . This is exactly the integrated and summed form of the single-interval Peano error of 43.09.01, whose leading constant $- h^{3} f^{''} (η) /12$ becomes the $h^{2}$ coefficient $\frac{B _{2}}{2 !} = \frac{1}{12}$ here once the $N$ local errors are added. The result builds toward two distinct payoffs and appears again in 43.09.03, where the same exactness-versus-smoothness exchange reappears as the degree- $2 n + 1$ accuracy of Gauss rules. The central insight is that an error expansion in known powers of $h$ is precisely what Richardson extrapolation needs: the even-power structure is the foundational reason the Romberg denominators are $4^{k} - 1$ and not $2^{k} - 1$ , and this is exactly why one trapezium-to-Simpson step gains two orders at once. Putting these together, the endpoint-only character of the error generalises into the spectral accuracy of the periodic trapezium rule, where the corrections all cancel and the error is dual to the decay of the integrand's Fourier coefficients.

Exercises Intermediate+

Exercise 3 (medium, symbolic).

Derive the global composite trapezium error $- \frac{( b - a ) h ^{2}}{12} f^{''} (η)$ by summing the single-interval Peano errors $- \frac{h ^{3}}{12} f^{''} (η_{j})$ from 43.09.01.

Hint

Sum over $j = 0, \dots, N - 1$ , factor out $- h^{3} /12$ , write the sum of $f^{''} (η_{j})$ as $N$ times an average, and use the intermediate value theorem on the continuous $f^{''}$ to replace the average by $f^{''} (η)$ .

Answer

On subinterval $j$ the single-interval trapezium error is $E_{j} = - \frac{h ^{3}}{12} f^{''} (η_{j})$ for some $η_{j}$ in that subinterval. Summing, the total error is $\sum_{j = 0}^{N - 1} E_{j} = - \frac{h ^{3}}{12} \sum_{j} f^{''} (η_{j}) = - \frac{h ^{3}}{12} N \cdot (\frac{1}{N} \sum_{j} f^{''} (η_{j}))$ . The bracketed average lies between $min f^{''}$ and $max f^{''}$ , so by the intermediate value theorem applied to the continuous $f^{''}$ there is $η \in (a, b)$ with $\frac{1}{N} \sum_{j} f^{''} (η_{j}) = f^{''} (η)$ . With $N h = b - a$ , the error is $- \frac{h ^{3}}{12} N f^{''} (η) = - \frac{( b - a ) h ^{2}}{12} f^{''} (η)$ , the $O (h^{2})$ global error. Rubric: full credit for the sum, the average-to- $f^{''} (η)$ step, and the $N h = b - a$ substitution.

Exercise 5 (medium, symbolic).

Use the one-panel Euler-Maclaurin identity to show that the midpoint and trapezium errors on $[0, 1]$ have opposite leading signs and a $1 : 2$ magnitude ratio, recovering the relation noted in 43.09.01.

Hint

The trapezium leading error is $\frac{B _{2}}{2 !} [g^{'} (1) - g^{'} (0)] = \frac{1}{12} [g^{'} (1) - g^{'} (0)] \approx \frac{1}{12} g^{''}$ . For the midpoint rule, an analogous Euler-Maclaurin expansion gives leading coefficient $- \frac{1}{24}$ ; compare.

Answer

For the trapezium rule the leading Euler-Maclaurin term is $\frac{B _{2}}{2 !} [g^{'} (1) - g^{'} (0)] = \frac{1}{12} [g^{'} (1) - g^{'} (0)]$ , and since $g^{'} (1) - g^{'} (0) = \int_{0}^{1} g^{''} d s \approx g^{''} (η)$ , the trapezium error is $\approx + \frac{1}{12} g^{''} (η)$ relative to the integral (equivalently the rule overshoots a convex $g$ ). The midpoint rule's leading error is $- \frac{1}{24} g^{''} (η)$ , half the magnitude and opposite sign. Their ratio $\frac{1/12}{1/24} = 2$ and opposite signs are exactly the bracketing relation of 43.09.01: midpoint and trapezium straddle the true integral, and Simpson is the weighted average $\frac{1}{3}$ (trapezium) $+ \frac{2}{3}$ (midpoint) that cancels the $g^{''}$ term. Rubric: full credit for both leading coefficients, the sign, and the $1 : 2$ ratio.

Exercise 6 (medium, numeric).

For the periodic integrand $f (x) = e^{c o s (2 π x)}$ on $[0, 1]$ , explain via Euler-Maclaurin why the composite trapezium rule converges faster than $O (h^{2})$ , and state what governs the residual error.

Hint

$f$ and all its derivatives are $1$ -periodic, so $f^{(2 k - 1)} (1) = f^{(2 k - 1)} (0)$ for every $k$ . What does that do to every term in the Euler-Maclaurin sum?

Answer

Since $f$ is smooth and $1$ -periodic, every derivative satisfies $f^{(2 k - 1)} (1) = f^{(2 k - 1)} (0)$ , so every endpoint-correction term $\frac{B _{2 k}}{( 2 k )!} h^{2 k} [f^{(2 k - 1)} (1) - f^{(2 k - 1)} (0)]$ vanishes identically. The error equals only the remainder $R_{m}$ for every $m$ , which can be driven below any power of $h$ ; for an analytic periodic $f$ the convergence is geometric. The residual error is governed by the decay of $f$ 's Fourier coefficients (equivalently the width of its strip of analyticity), not by $h$ polynomially. Rubric: full credit for the vanishing of every correction term and the Fourier/analyticity statement.

Exercise 7 (hard, symbolic).

Prove the one-panel Euler-Maclaurin identity to first order: for $g \in C^{2} [0, 1]$ , $\frac{g ( 0 ) + g ( 1 )}{2} - \int_{0}^{1} g = \frac{1}{12} [g^{'} (1) - g^{'} (0)] - \frac{1}{2} \int_{0}^{1} B_{2} ({s}) g^{''} (s) d s$ , where $B_{2} (s) = s^{2} - s + \frac{1}{6}$ .

Hint

Start from $\int_{0}^{1} g = \frac{g ( 0 ) + g ( 1 )}{2} - \int_{0}^{1} B_{1} (s) g^{'} (s) d s$ with $B_{1} (s) = s - \frac{1}{2}$ . Integrate the last integral by parts using $B_{1} = \frac{1}{2} B_{2}^{'}$ and $B_{2} (0) = B_{2} (1) = \frac{1}{6}$ .

Answer

Integration by parts with $B_{1} (s) = s - \frac{1}{2}$ , $B_{1}^{'} = 1$ gives $\int_{0}^{1} g = [B_{1} g]_{0}^{1} - \int_{0}^{1} B_{1} g^{'} = \frac{g ( 1 ) + g ( 0 )}{2} - \int_{0}^{1} B_{1} g^{'}$ . Hence $\frac{g ( 0 ) + g ( 1 )}{2} - \int_{0}^{1} g = \int_{0}^{1} B_{1} (s) g^{'} (s) d s$ . Now $B_{1} = \frac{1}{2} B_{2}^{'}$ since $B_{2}^{'} = 2 B_{1}$ , so $\int_{0}^{1} B_{1} g^{'} = \frac{1}{2} \int_{0}^{1} B_{2}^{'} (s) g^{'} (s) d s = \frac{1}{2} ([B_{2} g^{'}]_{0}^{1} - \int_{0}^{1} B_{2} g^{''})$ . The boundary term is $\frac{1}{2} B_{2} (1) g^{'} (1) - \frac{1}{2} B_{2} (0) g^{'} (0) = \frac{1}{12} [g^{'} (1) - g^{'} (0)]$ using $B_{2} (0) = B_{2} (1) = \frac{1}{6}$ . Therefore $\frac{g ( 0 ) + g ( 1 )}{2} - \int_{0}^{1} g = \frac{1}{12} [g^{'} (1) - g^{'} (0)] - \frac{1}{2} \int_{0}^{1} B_{2} (s) g^{''} (s) d s$ , which is the claim with ${s} = s$ on $[0, 1)$ . Rubric: full credit for the first integration by parts, the use of $B_{2}^{'} = 2 B_{1}$ , and the boundary evaluation with $B_{2} = 1/6$ .

Exercise 8 (hard, symbolic).

Justify the Romberg denominator $4^{k} - 1$ : assuming $T_{N} = I + c_{1} h^{2} + c_{2} h^{4} + \dots$ in even powers, show that the recurrence $T_{j}^{(k)} = T_{j}^{(k - 1)} + \frac{T _{j}^{(k - 1)} - T _{j - 1}^{(k - 1)}}{4 ^{k} - 1}$ eliminates the leading surviving error term at each step.

Hint

After $k - 1$ steps the leading error is $c_{k} h^{2 k}$ . Halving $h$ multiplies this term by $(1/2)^{2 k} = 4^{- k}$ . Form a linear combination of the two estimates that cancels $h^{2 k}$ .

Answer

Suppose after $k - 1$ extrapolations the column- $(k - 1)$ entries satisfy $T_{j}^{(k - 1)} = I + c_{k} h_{j}^{2 k} + O (h_{j}^{2 k + 2})$ with $h_{j} = h_{0} / 2^{j}$ , so $h_{j} = h_{j - 1} /2$ and $h_{j}^{2 k} = 4^{- k} h_{j - 1}^{2 k}$ . The two consecutive entries give $T_{j}^{(k - 1)} - I = c_{k} h_{j}^{2 k} + \dots$ and $T_{j - 1}^{(k - 1)} - I = c_{k} h_{j - 1}^{2 k} + \dots = 4^{k} c_{k} h_{j}^{2 k} + \dots$ . Their difference is $T_{j}^{(k - 1)} - T_{j - 1}^{(k - 1)} = c_{k} h_{j}^{2 k} (1 - 4^{k}) + O (h^{2 k + 2})$ , so $\frac{T _{j}^{(k - 1)} - T _{j - 1}^{(k - 1)}}{4 ^{k} - 1} = - c_{k} h_{j}^{2 k} + O (h^{2 k + 2})$ . Adding to $T_{j}^{(k - 1)}$ cancels the $c_{k} h_{j}^{2 k}$ term: $T_{j}^{(k)} = I + O (h_{j}^{2 k + 2})$ . The denominator $4^{k} - 1$ is dictated entirely by the even-power spacing of the trapezium expansion; an odd-power term would require $2^{p} - 1$ keyed to its order. Rubric: full credit for the $4^{- k}$ scaling, the difference computation, and the cancellation.

Advanced results Master

The composite rules fix the polynomial-order picture; the master layer turns the Euler-Maclaurin expansion into the precise statement that the trapezoidal rule is the canonical quadrature for analytic periodic integrands, places Romberg integration inside the general theory of Richardson extrapolation, and bounds the asymptotic nature of the expansion.

Theorem 1 (Euler-Maclaurin, full form and the asymptotic remainder). For $f \in C^{2 m + 2} [a, b]$ the composite trapezium rule satisfies $T_{N} (f) - \int_{a}^{b} f = \sum_{k = 1}^{m} \frac{B _{2 k}}{( 2 k )!} h^{2 k} [f^{(2 k - 1)} (b) - f^{(2 k - 1)} (a)] + R_{m}$ with $R_{m} = - \frac{h ^{2 m + 2}}{( 2 m + 2 )!} \int_{a}^{b} B_{2 m + 2} (\frac{x - a}{h}) f^{(2 m + 2)} (x) d x$ ^{[Davis-Rabinowitz Ch. 4]}. The series is generally divergent because $∣ B_{2 k} ∣ \sim 2 (2 k)! / (2 π)^{2 k}$ , so the coefficients $∣ B_{2 k} ∣/ (2 k)! \sim 2/ (2 π)^{2 k}$ multiply $h^{2 k} f^{(2 k - 1)}$ whose derivative factor can grow factorially; the expansion is asymptotic, and for fixed $f$ one truncates near the index minimising $∣ R_{m} ∣$ . The remainder admits the sharp bound $∣ R_{m} ∣ \leq \frac{2}{( 2 π ) ^{2 m + 2}} h^{2 m + 2} \int_{a}^{b} ∣ f^{(2 m + 2)} ∣$ , which exhibits the $(2 π)^{- (2 m + 2)}$ gain per order.

Theorem 2 (spectral accuracy of the trapezoidal rule on periodic integrands). Let $f$ be $L$ -periodic and $C^{\infty}$ , and let $T_{N}$ be the trapezoidal rule on $N$ equal panels over one period. Then every Euler-Maclaurin endpoint correction vanishes, $f^{(2 k - 1)} (a + L) = f^{(2 k - 1)} (a)$ , so $T_{N} (f) - \int_{0}^{L} f = O (h^{p})$ for every $p$ . If $f$ extends analytically to the strip $∣ Im z ∣ < σ$ , the error decays geometrically, $∣ T_{N} (f) - \int_{0}^{L} f ∣ \leq C e^{- 2 π σ N / L}$ ^{[Trefethen-Weideman 2014]}. Equivalently, writing $f = \sum_{n} \hat{f}_{n} e^{2 π in x / L}$ , the trapezoidal rule integrates each Fourier mode exactly except for aliasing, and the error is precisely the sum of aliased modes $\sum_{n \neq = 0, N ∣ n} \hat{f}_{n}$ , whose decay is governed by the analyticity strip. This is the structural reason spectral methods, contour-integral evaluation of special functions, and inversion of Laplace transforms all rest on the trapezoidal rule.

Theorem 3 (Romberg integration and the Richardson tableau). With $T_{j}^{(0)} = T_{2^{j}} (f)$ , the Romberg recurrence $T_{j}^{(k)} = T_{j}^{(k - 1)} + (T_{j}^{(k - 1)} - T_{j - 1}^{(k - 1)}) / (4^{k} - 1)$ produces, in exact arithmetic, $T_{j}^{(k)} = \int_{a}^{b} f + O (h_{j}^{2 k + 2})$ for $f \in C^{2 k + 2}$ ^{[Davis-Rabinowitz Ch. 4]}. Each column is a polynomial extrapolation of the trapezium estimates to step $h = 0$ along the variable $h^{2}$ , exact for the truncated even-power model; the first column is composite Simpson, the second is composite Boole's rule, and the diagonal $T_{k}^{(k)}$ converges superalgebraically for analytic $f$ . The construction is the quadrature instance of the general Richardson extrapolation principle: given any quantity $A (h) = A_{0} + \sum_{i} a_{i} h^{p_{i}}$ with known exponents $p_{1} < p_{2} < \dots$ , the Neville-type tableau with denominators $2^{p_{k}} - 1$ removes the terms in order, and the even-power $p_{k} = 2 k$ of the trapezium rule specialises this to $4^{k} - 1$ .

Theorem 4 (adaptive quadrature: local estimation and refinement). Adaptive schemes estimate the local error on a subinterval $[α, β]$ by comparing a base rule $Q$ with a higher-order or subdivided rule $Q^{'}$ : the difference $Q^{'} (f) - Q (f)$ estimates the error of the cruder rule by Richardson's principle, and the interval is accepted when $∣ Q^{'} (f) - Q (f) ∣ \leq ε (β - α) / (b - a)$ (tolerance apportioned by length) and bisected otherwise, recursing on the halves ^{[Davis-Rabinowitz Ch. 6]}. For composite Simpson the standard estimate uses $S_{[α, β]}$ against $S_{[α, μ]} + S_{[μ, β]}$ at the midpoint $μ$ , giving the error estimate $\frac{1}{15} (S_{two halves} - S_{whole})$ from the $O (h^{4})$ model; the factor $1/15 = 1/ (2^{4} - 1)$ is the Richardson denominator at order four. The scheme concentrates nodes where $f^{(4)}$ is large, achieving a target accuracy with far fewer evaluations than a uniform grid when $f$ has localised features, at the cost of error estimates that can mislead on integrands rougher than the rule's model assumes.

Synthesis. The Euler-Maclaurin formula is the foundational reason quadrature on a partition is governed by a single structured object: the composite trapezium error is an asymptotic series in even powers of $h$ whose coefficients are Bernoulli numbers times endpoint-derivative jumps, all interior contributions telescoped away. The central insight is that this even-power, endpoint-only structure is the source of three phenomena at once: the leading error is $O (h^{2})$ with the same constant $1/12$ that the single-interval Peano kernel of 43.09.01 produces; the rule is spectrally accurate on periodic integrands, where the endpoint jumps vanish and the error collapses to the aliased Fourier tail; and Richardson extrapolation with denominators $4^{k} - 1$ removes the error two orders at a time, the Romberg table. This is exactly the exactness-versus-smoothness exchange of the chapter, dual to the Peano-kernel error of 43.09.01 and to the degree- $2 n + 1$ accuracy that Gauss quadrature buys in 43.09.03 by abandoning equispacing. Putting these together, the periodic trapezoidal rule's geometric convergence and the Romberg diagonal's superalgebraic convergence are one statement read at two resolutions: once smoothness removes the algebraic error terms, the residual is controlled by complex-analytic data — the strip of analyticity, the Fourier decay — not by any finite power of $h$ . The bridge from the local single-panel error to this global spectral picture is the integration-by-parts against the Bernoulli polynomials, which generalises the single mean value theorem of the Peano analysis into a full asymptotic ladder.

Full proof set Master

Proposition 1 (composite trapezium global error). For $f \in C^{2} [a, b]$ and $h = (b - a) / N$ , $T_{N} (f) - \int_{a}^{b} f = - \frac{( b - a ) h ^{2}}{12} f^{''} (η)$ for some $η \in (a, b)$ .

Proof. On subinterval $[x_{j}, x_{j + 1}]$ the single-interval trapezium error of 43.09.01 is $\int_{x_{j}}^{x_{j + 1}} f - \frac{h}{2} (f_{j} + f_{j + 1}) = - \frac{h ^{3}}{12} f^{''} (η_{j})$ for some $η_{j} \in (x_{j}, x_{j + 1})$ , since the trapezium Peano kernel $K_{2} = \frac{1}{2} (t - x_{j}) (t - x_{j + 1}) \leq 0$ is one-signed. Summing over $j = 0, \dots, N - 1$ and using that the composite rule is the sum of the panel rules, $\int_{a}^{b} f - T_{N} (f) = - \frac{h ^{3}}{12} \sum_{j = 0}^{N - 1} f^{''} (η_{j})$ . The quantity $\frac{1}{N} \sum_{j} f^{''} (η_{j})$ lies in $[min_{[a, b]} f^{''}, max_{[a, b]} f^{''}]$ , and $f^{''}$ is continuous, so by the intermediate value theorem there is $η \in (a, b)$ with $\frac{1}{N} \sum_{j} f^{''} (η_{j}) = f^{''} (η)$ . Thus $\int_{a}^{b} f - T_{N} (f) = - \frac{h ^{3}}{12} N f^{''} (η) = - \frac{( b - a ) h ^{2}}{12} f^{''} (η)$ , i.e. $T_{N} (f) - \int_{a}^{b} f = + \frac{( b - a ) h ^{2}}{12} f^{''} (η)$ after the sign convention; the displayed form reflects $T_{N}$ minus the integral. $□$

Proposition 2 (one-panel Euler-Maclaurin identity). For $g \in C^{2 m + 2} [0, 1]$ , $\frac{g ( 0 ) + g ( 1 )}{2} - \int_{0}^{1} g = \sum_{k = 1}^{m} \frac{B _{2 k}}{( 2 k )!} [g^{(2 k - 1)} (1) - g^{(2 k - 1)} (0)] - \frac{1}{( 2 m + 2 )!} \int_{0}^{1} B_{2 m + 2} (s) g^{(2 m + 2)} (s) d s$ .

Proof. By induction on $m$ . Base $m = 0$ is the single integration by parts: with $B_{1} (s) = s - \frac{1}{2}$ , $\int_{0}^{1} g = [B_{1} g]_{0}^{1} - \int_{0}^{1} B_{1} g^{'} = \frac{g ( 1 ) + g ( 0 )}{2} - \int_{0}^{1} B_{1} g^{'}$ , and $B_{1} = B_{1}$ on $[0, 1)$ , giving the claim with empty sum and remainder index $2$ after one more step. For the inductive step, take the remainder integral $- \frac{1}{( 2 k )!} \int_{0}^{1} B_{2 k} (s) g^{(2 k)} (s) d s$ and integrate by parts twice using $B_{r}^{'} = r B_{r - 1}$ . The first parts step, with $B_{2 k} = \frac{1}{2 k + 1} B_{2 k + 1}^{'}$ , produces a boundary term $\frac{1}{( 2 k + 1 )!} [B_{2 k + 1} (s) g^{(2 k)} (s)]_{0}^{1}$ which vanishes since $B_{2 k + 1} (0) = B_{2 k + 1} (1) = 0$ , and an integral $+ \frac{1}{( 2 k + 1 )!} \int_{0}^{1} B_{2 k + 1} g^{(2 k + 1)}$ . The second parts step, with $B_{2 k + 1} = \frac{1}{2 k + 2} B_{2 k + 2}^{'}$ , produces the boundary term $\frac{1}{( 2 k + 2 )!} [B_{2 k + 2} g^{(2 k + 1)}]_{0}^{1} = \frac{B _{2 k + 2}}{( 2 k + 2 )!} [g^{(2 k + 1)} (1) - g^{(2 k + 1)} (0)]$ , using $B_{2 k + 2} (0) = B_{2 k + 2} (1) = B_{2 k + 2}$ , plus the next remainder $- \frac{1}{( 2 k + 2 )!} \int_{0}^{1} B_{2 k + 2} g^{(2 k + 2)}$ . Thus each pair of integrations adds one Bernoulli term and advances the remainder index by two; after $m$ such advances the stated identity holds. $□$

Proposition 3 (Euler-Maclaurin on $[a, b]$ by panel summation). The composite-rule formula of the Key theorem holds, with the interior endpoint corrections telescoping to the global endpoints.

Proof. Rescale panel $[x_{j}, x_{j + 1}]$ by $x = x_{j} + h s$ , so $f (x) = g_{j} (s)$ with $g_{j} (s) = f (x_{j} + h s)$ and $g_{j}^{(r)} (s) = h^{r} f^{(r)} (x_{j} + h s)$ . Proposition 2 on $g_{j}$ gives $\frac{f _{j} + f _{j + 1}}{2} - \frac{1}{h} \int_{x_{j}}^{x_{j + 1}} f = \sum_{k = 1}^{m} \frac{B _{2 k}}{( 2 k )!} h^{2 k - 1} [f^{(2 k - 1)} (x_{j + 1}) - f^{(2 k - 1)} (x_{j})] - \frac{h ^{2 m + 1}}{( 2 m + 2 )!} \int_{x_{j}}^{x_{j + 1}} B_{2 m + 2} (\frac{x - x _{j}}{h}) f^{(2 m + 2)} d x$ . Multiply by $h$ and sum over $j$ : the left side becomes $T_{N} (f) - \int_{a}^{b} f$ , the Bernoulli sums telescope because $f^{(2 k - 1)} (x_{j + 1})$ from panel $j$ cancels $- f^{(2 k - 1)} (x_{j})$ from panel $j + 1$ , leaving only $f^{(2 k - 1)} (b) - f^{(2 k - 1)} (a)$ , and the remainder integrals assemble into the single integral over $[a, b]$ with the $1$ -periodic $B_{2 m + 2}$ . The bound follows from $max_{s} ∣ B_{2 m + 2} (s) ∣ = ∣ B_{2 m + 2} ∣$ for the even periodic Bernoulli function. $□$

Proposition 4 (Richardson extrapolation removes the leading order). If $A (h) = A_{0} + \sum_{i \geq 1} a_{i} h^{p_{i}}$ with $0 < p_{1} < p_{2} < \dots$ known, then $\frac{2 ^{p_{1}} A ( h /2 ) - A ( h )}{2 ^{p_{1}} - 1} = A_{0} + O (h^{p_{2}})$ ; for $p_{i} = 2 i$ the factor is $4 - 1 = 3$ , and iterating gives the Romberg recurrence.

Proof. Compute $2^{p_{1}} A (h /2) = 2^{p_{1}} A_{0} + 2^{p_{1}} \sum_{i} a_{i} (h /2)^{p_{i}} = 2^{p_{1}} A_{0} + \sum_{i} a_{i} 2^{p_{1} - p_{i}} h^{p_{i}}$ . Subtract $A (h) = A_{0} + \sum_{i} a_{i} h^{p_{i}}$ : the $A_{0}$ terms give $(2^{p_{1}} - 1) A_{0}$ , and the $i = 1$ term gives $a_{1} (2^{p_{1} - p_{1}} - 1) h^{p_{1}} = a_{1} (1 - 1) h^{p_{1}} = 0$ . Hence $2^{p_{1}} A (h /2) - A (h) = (2^{p_{1}} - 1) A_{0} + \sum_{i \geq 2} a_{i} (2^{p_{1} - p_{i}} - 1) h^{p_{i}}$ , and dividing by $2^{p_{1}} - 1$ yields $A_{0} + O (h^{p_{2}})$ . For the trapezium expansion $p_{i} = 2 i$ , $2^{p_{1}} = 4$ , the denominator is $3$ , and the same step applied to the column- $(k - 1)$ entries (whose leading exponent is $2 k$ , so $2^{2 k} = 4^{k}$ ) gives the denominator $4^{k} - 1$ , the Romberg recurrence. $□$

Proposition 5 (periodic spectral accuracy via aliasing). For $f$ smooth and $1$ -periodic with Fourier series $f (x) = \sum_{n \in Z} \hat{f}_{n} e^{2 π in x}$ , the $N$ -point trapezoidal rule $T_{N} (f) = \frac{1}{N} \sum_{j = 0}^{N - 1} f (j / N)$ satisfies $T_{N} (f) - \int_{0}^{1} f = \sum_{n \neq = 0, N ∣ n} \hat{f}_{n}$ .

Proof. By linearity it suffices to evaluate $T_{N}$ on a single mode $e_{n} (x) = e^{2 π in x}$ . Then $T_{N} (e_{n}) = \frac{1}{N} \sum_{j = 0}^{N - 1} e^{2 π inj / N}$ , a geometric sum equal to $1$ when $N ∣ n$ (every term is $1$ ) and $0$ otherwise (the $N$ -th roots of unity sum to zero). Meanwhile $\int_{0}^{1} e_{n} = 1$ if $n = 0$ and $0$ otherwise. Subtracting, $T_{N} (e_{n}) - \int_{0}^{1} e_{n}$ equals $1$ exactly when $N ∣ n$ and $n \neq = 0$ , and $0$ otherwise. Summing against the Fourier coefficients, $T_{N} (f) - \int_{0}^{1} f = \sum_{n \neq = 0, N ∣ n} \hat{f}_{n}$ . If $f$ is analytic in the strip $∣ Im z ∣ < σ$ then $∣ \hat{f}_{n} ∣ \leq C e^{- 2 π σ ∣ n ∣}$ , so the aliased tail starting at $∣ n ∣ = N$ is bounded by a geometric series with ratio $e^{- 2 π σ}$ , giving $∣ T_{N} (f) - \int_{0}^{1} f ∣ \leq C^{'} e^{- 2 π σ N}$ . This is the geometric (spectral) convergence; the Euler-Maclaurin endpoint corrections vanish precisely because periodicity forces every $f^{(2 k - 1)}$ to match at the period ends, consistent with this aliasing computation. $□$

Connections Master

The composite trapezium and Simpson rules of this unit are the summed forms of the single-interval Newton-Cotes rules of 43.09.01: the global $O (h^{2})$ and $O (h^{4})$ errors are obtained by adding the one-signed Peano-kernel errors $- h^{3} f^{''} (η_{j}) /12$ across the partition, and the Euler-Maclaurin leading coefficient $B_{2} /2! = 1/12$ is exactly the single-interval error constant $\int K_{2} = E (x^{2}) /2!$ computed there; the positivity of all composite weights is the stability that the negative-weight high-order rules of that unit lack.
Gauss quadrature in 43.09.03 is the alternative response to the same accuracy demand this unit meets by subdivision: where the composite strategy fixes the order at $1$ or $2$ and shrinks $h$ to drive the Euler-Maclaurin error down, the Gauss strategy keeps a single interval and raises the degree of exactness to $2 n + 1$ by placing nodes at orthogonal-polynomial roots, so the chapter's two repairs of the Peano-kernel ceiling appear side by side, both governed by the exactness-versus-smoothness exchange visible in the even-power error structure here.
The Romberg construction is the quadrature instance of Richardson extrapolation, the same acceleration principle that produces high-order finite-difference stencils and the deferred-correction methods used for stiff ordinary differential equations; the requirement that the underlying expansion run in known powers of $h$ , satisfied here by the even-power Euler-Maclaurin series, is exactly the hypothesis those methods need, linking this unit to the convergence-acceleration toolkit used throughout numerical analysis.
The spectral accuracy of the periodic trapezoidal rule is the foundation of Fourier spectral methods and the discretised Cauchy/contour-integral evaluation of analytic functions; the aliasing identity $T_{N} (f) - \int f = \sum_{N ∣ n, n \neq = 0} \hat{f}_{n}$ is the quadrature face of the discrete-Fourier-transform sampling theorem, tying this unit to harmonic analysis and to the analytic-function-theory machinery that controls Fourier-coefficient decay through the strip of analyticity.
The integration-by-parts ladder against the Bernoulli polynomials that proves Euler-Maclaurin is the same mechanism by which the Bernoulli numbers enter the values $ζ (2 k)$ of the Riemann zeta function and the asymptotics of factorials through Stirling's series; the Euler-Maclaurin formula was Euler's original tool for both summing slowly convergent series and evaluating $ζ (2)$ , connecting this numerical unit to analytic number theory and the theory of asymptotic expansions.

Historical & philosophical context Master

The summation formula was found independently by Leonhard Euler, who announced it in a 1738 paper in the Commentarii of the Saint Petersburg Academy as a general method for summing slowly convergent series, and by Colin Maclaurin, who gave it in his Treatise of Fluxions (1742) within the fluxional calculus ^{[Euler 1738]}. Euler used it to compute the constant now bearing the name of Lorenzo Mascheroni and to evaluate $ζ (2)$ to many digits before he found the closed form $π^{2} /6$ . The systematic appearance of the Bernoulli numbers, introduced by Jacob Bernoulli in the Ars Conjectandi (posthumous, 1713) for sums of powers, as the coefficients of the expansion, and the recognition that the resulting series is asymptotic rather than convergent, were clarified through the eighteenth and nineteenth centuries; the integral remainder in terms of the periodic Bernoulli functions is the modern rigorous form.

Romberg integration was introduced by Werner Romberg in a 1955 note, Vereinfachte numerische Integration, in the proceedings of the Norwegian academy, recasting the centuries-old Euler-Maclaurin expansion as a practical extrapolation tableau; its convergence was analysed by Friedrich Bauer, Heinz Rutishauser, and Eduard Stiefel shortly after, placing it within the general Richardson-extrapolation framework that Lewis Fry Richardson had proposed in 1911 under the name "the deferred approach to the limit." The spectral accuracy of the trapezoidal rule on periodic and analytic integrands, long known to the contour-integration community, was surveyed and sharpened by Lloyd N. Trefethen and J. A. C. Weideman in 2014, who collected the geometric-convergence estimates and their applications to special-function computation ^{[Trefethen-Weideman 2014]}. Adaptive quadrature, comparing a rule against its refinement to drive recursive subdivision, became standard with the automatic-integration routines of the 1960s and 1970s, notably those of Robert Piessens and collaborators in the QUADPACK library.

Bibliography Master

@book{sulimayers2003,
  author    = {S\"{u}li, Endre and Mayers, David F.},
  title     = {An Introduction to Numerical Analysis},
  publisher = {Cambridge University Press},
  year      = {2003}
}

@book{davisrabinowitz1984,
  author    = {Davis, Philip J. and Rabinowitz, Philip},
  title     = {Methods of Numerical Integration},
  edition   = {2},
  publisher = {Academic Press},
  year      = {1984}
}

@article{trefethenweideman2014,
  author  = {Trefethen, Lloyd N. and Weideman, J. A. C.},
  title   = {The Exponentially Convergent Trapezoidal Rule},
  journal = {SIAM Review},
  volume  = {56},
  number  = {3},
  year    = {2014},
  pages   = {385--458}
}

@article{romberg1955,
  author  = {Romberg, Werner},
  title   = {Vereinfachte numerische Integration},
  journal = {Det Kongelige Norske Videnskabers Selskabs Forhandlinger},
  volume  = {28},
  year    = {1955},
  pages   = {30--36}
}

@book{henrici1977,
  author    = {Henrici, Peter},
  title     = {Applied and Computational Complex Analysis, Volume 2},
  publisher = {Wiley},
  year      = {1977}
}

@article{euler1738,
  author  = {Euler, Leonhard},
  title   = {Methodus generalis summandi progressiones},
  journal = {Commentarii Academiae Scientiarum Petropolitanae},
  volume  = {6},
  year    = {1738},
  pages   = {68--97}
}

@book{bernoulli1713,
  author    = {Bernoulli, Jacob},
  title     = {Ars Conjectandi},
  publisher = {Thurnisius, Basel},
  year      = {1713}
}

Prerequisites

43.09.01

Tier anchors

beginner: Süli-Mayers 2003 *An Introduction to Numerical Analysis* (Cambridge) §7.5-7.7 (composite trapezium and Simpson rules and their global error, at the level of a first numerical-methods course); Burden-Faires 2011 *Numerical Analysis* 9e (Brooks/Cole) §4.4-4.6 for the composite-rule bookkeeping and adaptive Simpson on worked integrands
intermediate: Süli-Mayers 2003 *An Introduction to Numerical Analysis* (Cambridge) §7.5-7.7 (composite rules, the Euler-Maclaurin expansion of the trapezium error in even powers of h with Bernoulli-number coefficients, Romberg integration by Richardson extrapolation); Stoer-Bulirsch 2002 *Introduction to Numerical Analysis* 3e (Springer) §3.2-3.4 (Euler-Maclaurin summation, the Romberg T-table, and adaptive subdivision)
master: Davis-Rabinowitz 1984 *Methods of Numerical Integration* 2e (Academic Press) Ch. 2-4, 6 (composite rules, Euler-Maclaurin, Romberg extrapolation and the periodic/trapezoidal spectral accuracy); Trefethen-Weideman 2014 *The exponentially convergent trapezoidal rule* (SIAM Review 56) (geometric convergence for analytic periodic integrands); Henrici 1977 *Applied and Computational Complex Analysis* vol. 2 (Wiley) Ch. 11 (Euler-Maclaurin, Bernoulli numbers, and the trapezoidal rule on the circle)

References

Süli, E. & Mayers, D. F. — An Introduction to Numerical Analysis · Cambridge University Press (2003). Chapter 7 (§7.5-7.7) builds the composite rules by partitioning [a,b] into N equal subintervals of width h=(b-a)/N and applying the basic rule on each: the composite trapezium rule T_N(f)=h[f_0/2 + f_1 + ... + f_{N-1} + f_N/2] with global error -(b-a)h^2 f''(eta)/12 = O(h^2), the composite Simpson rule with global error O(h^4). §7.6 derives the Euler-Maclaurin summation formula, expressing the difference between the composite trapezium sum and the exact integral as an asymptotic series in even powers of h, int_a^b f - T_N(f) = -sum_{k=1}^m B_{2k}/(2k)! h^{2k}[f^{(2k-1)}(b)-f^{(2k-1)}(a)] - remainder, with B_{2k} the Bernoulli numbers; the endpoint-derivative corrections are the only obstruction to high accuracy, so for a smooth periodic integrand (all periodic derivatives matching at the endpoints) the trapezium rule converges faster than any power of h. §7.7 uses the even-power structure to justify Romberg integration: repeated Richardson extrapolation of T_N on a sequence of halved step sizes eliminates the h^2, h^4, ... error terms in turn, the T-table, and discusses adaptive subdivision driven by a local error estimate.
Davis, P. J. & Rabinowitz, P. — Methods of Numerical Integration · Academic Press, 2nd edition (1984). Chapters 2-4 and 6 give the comprehensive composite-rule and extrapolation account: the Euler-Maclaurin formula with the Bernoulli-number coefficients and the Bernoulli-polynomial remainder term, the demonstration that the trapezoidal rule is the optimal rule for smooth periodic integrands (the periodic Euler-Maclaurin endpoint corrections all cancel, leaving an error governed only by the smoothness/decay of the Fourier coefficients), the Romberg T-table construction R_{j,k}=R_{j,k-1}+(R_{j,k-1}-R_{j-1,k-1})/(4^k-1) as repeated Richardson extrapolation against the known h^{2k} expansion, the convergence of the Romberg columns and the diagonal, and adaptive quadrature schemes that compare a rule against its subdivided refinement to estimate the local error and recursively bisect intervals where a per-interval tolerance is exceeded, with global tolerance apportioned across subintervals.
Trefethen, L. N. & Weideman, J. A. C. — The exponentially convergent trapezoidal rule · SIAM Review 56(3), 385-458 (2014). A survey establishing that for a function analytic and periodic on the real line (equivalently analytic in a strip around the real axis), the N-point trapezoidal rule over one period converges geometrically, the error decaying like exp(-2 pi a N / L) where a is the half-width of the strip of analyticity and L the period; the paper traces this to the Euler-Maclaurin endpoint corrections all vanishing under periodicity, so that the error equals the sum of the omitted (aliased) Fourier modes, and connects the result to contour integration, the discretised Cauchy integral, inverse Laplace transforms, and the computation of special functions, with the trapezoidal rule emerging as the natural and near-optimal quadrature for analytic periodic integrands.
Euler, L. — Methodus generalis summandi progressiones (and Maclaurin, C. — A Treatise of Fluxions) · Euler, Commentarii Academiae Scientiarum Petropolitanae 6 (1738), 68-97; Maclaurin, A Treatise of Fluxions (Edinburgh, 1742), Book II. The summation formula relating a finite sum sum_{k} f(k) to the integral int f plus a series of endpoint-derivative corrections weighted by the Bernoulli numbers B_{2k}/(2k)!, discovered independently by Euler (for accelerating the summation of slowly convergent series and evaluating zeta(2)) and Maclaurin (in the fluxional calculus); the asymptotic, generally divergent, character of the series and the integral remainder term involving the periodic Bernoulli functions were made precise in the nineteenth and twentieth centuries.

Estimated time

beginner: 20m
intermediate: 50m
master: 90m