2.8.2 Trapezoid and Simpson rules

In this lesson

The trapezoid rule: a straight top on every strip, and the average of the left and right sums of lesson 2.5.1
Its error, which shrinks like $\Delta x^2$ and is about twice the midpoint error, with the opposite sign
Simpson's rule: a parabola through each pair of strips, with weights 1, 4, 2, 4, …, 4, 1
Simpson's rule as $\frac{2 \times \text{midpoint} + \text{trapezoid}}{3}$ , exact for cubics, with an error that shrinks like $\Delta x^4$
Where the rules break: a corner, a curve that blows up at an end, an odd number of strips
SciPy's trapezoid and simpson, for values you already have

You need Riemann sums from lesson 2.5.1 and the fundamental theorem from lesson 2.5.3. The error rules use second derivatives and bending (lesson 2.4.2).

The trapezoid rule

Lesson 2.5.1 gave each strip a flat top, at the height of one point. The trapezoid rule joins the heights at the two edges of the strip with a straight line instead. Each strip becomes a trapezoid, whose area is the average of its two heights times its width.

As in lesson 2.5.1, cut $[a, b]$ into $n$ strips of width $\Delta x = \frac{b - a}{n}$ , with edges $x_i = a + i\,\Delta x$ . Write $L_n$ , $R_n$ and $M_n$ for the left, right and midpoint sums, and $T_n$ for the trapezoid rule.

Rule: Trapezoid rule

\begin{gathered} T_n = \Delta x \big[\tfrac{1}{2}f(x_0) + f(x_1) + \cdots \\ + f(x_{n-1}) + \tfrac{1}{2}f(x_n)\big] \\[6pt] T_n = \tfrac{1}{2}\,(L_n + R_n) \end{gathered}

Every inner edge belongs to two strips, so its height counts in full; the two outer edges count half. The left sum uses the heights at $x_0$ to $x_{n-1}$ , and the right sum those at $x_1$ to $x_n$ . Their average weights the heights the same way.

Like any area, $T_n$ is in height units times width units. Take $x^2$ on $[0, 1]$ with 4 strips. The heights at the edges are $0, \frac{1}{16}, \frac{1}{4}, \frac{9}{16}$ and 1, so

\begin{aligned} T_4 &= \tfrac{1}{4}\left(\tfrac{0}{2} + \tfrac{1}{16} + \tfrac{4}{16} + \tfrac{9}{16} + \tfrac{1}{2}\right) \\ &= \tfrac{1}{4} \cdot \tfrac{22}{16} = 0.34375 \end{aligned}

That is the average of lesson 2.5.1's left sum 0.21875 and right sum 0.46875. The exact area is $\frac{1}{3}$ , so $T_4$ is too large by 0.0104. For $\frac{1}{x}$ on $[1, 2]$ , the left and right sums 0.759524 and 0.634524 average to 0.697024, against the exact $\ln 2 \approx 0.693147$ .

The trapezoid rule with 4 strips on $y = \sin x$ over $[0, \pi]$ . The curve bends down, so every straight top runs under it and the total, 1.8961, is short of 2.

T_4

import numpy as np def trapezoid_rule(f, a, b, n): dx = (b - a) / n y = f(a + dx * np.arange(n + 1)) return dx * (y[0] / 2 + np.sum(y[1:-1]) + y[-1] / 2) def riemann(f, a, b, n, rule): dx = (b - a) / n offset = {"left": 0, "right": 1, "mid": 0.5}[rule] return np.sum(f(a + (np.arange(n) + offset) * dx)) * dx def recip(x): return 1 / x for f, a, b in [(np.square, 0, 1), (recip, 1, 2)]: left, right = riemann(f, a, b, 4, "left"), riemann(f, a, b, 4, "right") print(round(trapezoid_rule(f, a, b, 4), 6), round(left, 6), round(right, 6), round((left + right) / 2, 6)) print(round(trapezoid_rule(np.square, 0, 1, 4) - 1 / 3, 4))

⌘/Ctrl + Enter runs the cell

The trapezoid error, and the midpoint's

The trapezoid error shrinks like $\Delta x^2$ , as the midpoint error does (lesson 2.5.1). With 10 strips, against the exact integrals of lesson 2.5.1:

integral	midpoint error	trapezoid error
$x^2$ on $[0, 1]$	$-8.33 \times 10^{-4}$	$1.67 \times 10^{-3}$
$\frac{1}{x}$ on $[1, 2]$	$-3.12 \times 10^{-4}$	$6.24 \times 10^{-4}$
$\sin x$ on $[0, \pi]$	$8.25 \times 10^{-3}$	$-1.65 \times 10^{-2}$
$e^x$ on $[0, 1]$	$-7.16 \times 10^{-4}$	$1.43 \times 10^{-3}$

In every row the midpoint error is about half the trapezoid error, with the opposite sign. With 100 strips each error is a hundred times smaller.

The sign comes from the bend (lesson 2.4.2). Where the curve bends up, each straight top lies above it, so the trapezoid rule is too large.

The midpoint rectangle has the same area as the strip under the tangent at the midpoint. That tangent lies below the curve, so the midpoint sum is too small.

Rule: Trapezoid and midpoint errors

\begin{gathered} T_n - I \approx \tfrac{1}{12}\,\Delta x^2\,B \\[4pt] M_n - I \approx -\tfrac{1}{24}\,\Delta x^2\,B \\[4pt] B = f'(b) - f'(a) \end{gathered}

$I$ is the exact integral. $B$ is the integral of $f''$ (lesson 2.5.3): the total bend over $[a, b]$ . Bending up makes it positive, so the trapezoid rule is too large and the midpoint sum too small, by half as much.

For $x^2$ on $[0, 1]$ , $B = f'(1) - f'(0) = 2$ , so the rule gives $\frac{2}{12}\Delta x^2 = \frac{1}{96} \approx 0.0104$ for 4 strips: exactly the error found above. For $\sin x$ on $[0, \pi]$ the total bend is $\cos\pi - \cos 0 = -2$ . So the trapezoid rule is too small there, as the figure above shows, and the midpoint sum too large.

Whether $f$ rises or falls does not matter, unlike for the left and right sums: only the bend does. A straight line has no bend, so both rules are exact for it.

B = f'(b) - f'(a)

import math import numpy as np def trapezoid_rule(f, a, b, n): dx = (b - a) / n y = f(a + dx * np.arange(n + 1)) return dx * (y[0] / 2 + np.sum(y[1:-1]) + y[-1] / 2) def midpoint(f, a, b, n): dx = (b - a) / n return dx * np.sum(f(a + dx * (np.arange(n) + 0.5))) examples = [ ("x^2", np.square, 0, 1, 1 / 3, 2), ("1/x", lambda x: 1 / x, 1, 2, math.log(2), -1 / 4 + 1), ("sin", np.sin, 0, math.pi, 2, -2), ("e^x", np.exp, 0, 1, math.e - 1, math.e - 1), ] for n in [10, 100]: for name, f, a, b, exact, bend in examples: dx = (b - a) / n m, t = midpoint(f, a, b, n) - exact, trapezoid_rule(f, a, b, n) - exact print(n, name, f"{m:.2e}", f"{t:.2e}", f"{dx**2 / 12 * bend:.2e}")

⌘/Ctrl + Enter runs the cell

Simpson's rule

Straight tops cannot follow a bend. Simpson's rule replaces them by curves: over each pair of strips, the parabola through the three heights.

Over one pair, from $x_0$ to $x_2$ , that parabola has area $\frac{\Delta x}{3}\big(f(x_0) + 4f(x_1) + f(x_2)\big)$ , where $\Delta x$ is the width of one strip. Adding the pairs gives the weights. The middle of each pair counts 4, an edge shared by two pairs counts $1 + 1 = 2$ , and the two ends count 1.

Simpson's rule with one pair of strips on $y = \sin x$ : $\frac{\pi/2}{3}(0 + 4 + 0) = 2.0944$ . Three points already give the area to within 5%.

Rule: Simpson's rule

\begin{gathered} S_n = \frac{\Delta x}{3}\big[f(x_0) + 4f(x_1) \\ + 2f(x_2) + 4f(x_3) + \cdots \\ + 4f(x_{n-1}) + f(x_n)\big] \end{gathered}

Weights 1, 4, 2, 4, …, 2, 4, 1, times $\frac{\Delta x}{3}$ . The number of strips $n$ must be even, so that the strips pair up.

For $\frac{1}{x}$ on $[1, 2]$ with 4 strips, $\Delta x = 0.25$ and the heights at 1, 1.25, 1.5, 1.75 and 2 are 1, 0.8, 0.666667, 0.571429 and 0.5. So

\begin{aligned} S_4 &= \tfrac{0.25}{3}\big(1 + 3.2 + 1.333333 \\ &\quad + 2.285714 + 0.5\big) \\ &= 0.693254 \end{aligned}

It is off by $1.07 \times 10^{-4}$ . Lesson 2.5.1's midpoint sum with the same 4 strips, 0.69122, is off by $-1.93 \times 10^{-3}$ : eighteen times as much, the other way.

Try the rules side by side below, starting from these four strips of $\frac{1}{x}$ . Switch between midpoints, trapezoids and Simpson's parabolas, then raise $n$ to 10 and 100 and compare how fast each error shrinks.

Explore

Add strips and change how each strip's top is drawn: flat, a straight line (trapezoids) or a parabola through each pair (Simpson). Compare how fast each error shrinks.

Strips n = 4

Width of each: 1 ÷ 4 = 0.25
Total of 4 strips: 0.69325397
Exact area: ln 2 = 0.693147
Error: 0.000107 (an overestimate)

f_0

import numpy as np import sympy as sp x, h, f0, f1, f2 = sp.symbols("x h f0 f1 f2") parabola = sp.interpolate([(0, f0), (h, f1), (2 * h, f2)], x) print(sp.factor(sp.integrate(parabola, (x, 0, 2 * h)))) def simpson_rule(f, a, b, n): if n % 2: raise ValueError(f"Simpson's rule needs an even number of strips, not {n}") dx = (b - a) / n w = np.ones(n + 1) w[1:-1:2], w[2:-1:2] = 4, 2 return dx / 3 * np.sum(w * f(a + dx * np.arange(n + 1))) s4, m4 = simpson_rule(lambda t: 1 / t, 1, 2, 4), 0.69122 print(round(s4, 6), f"{s4 - np.log(2):.2e}", f"{m4 - np.log(2):.2e}", round((m4 - np.log(2)) / (s4 - np.log(2)), 1))

⌘/Ctrl + Enter runs the cell

Simpson's rule from the other two

The midpoint error is about $-\frac{1}{2}$ times the trapezoid error. So in $2M_m + T_m$ the two errors cancel, to leading order.

Rule: Simpson from midpoint and trapezoid

S_{2m} = \frac{2M_m + T_m}{3}

With $m$ strips, the trapezoid rule uses the edges and the midpoint sum the middles. Together they are the $2m + 1$ points of Simpson's rule with $2m$ strips, and the weights work out to 1, 4, 2, …, 4, 1.

For $\frac{1}{x}$ on $[1, 2]$ with $m = 4$ : $\frac{2 \times 0.69122 + 0.697024}{3} = 0.693155$ , which is $S_8$ , off by only $7.4 \times 10^{-6}$ .

Combining two estimates so that their leading errors cancel is called Richardson extrapolation. Halving the strips cuts the trapezoid error about 4 times, so $\frac{4T_8 - T_4}{3}$ gives $S_8$ as well.

Its error shrinks like $\Delta x^4$

Simpson's rule beats the midpoint sum on lesson 2.5.1's examples, and by more as the strips get thinner:

integral	$M_4$ error	$S_4$ error	$M_{100}$ error	$S_{100}$ error
$x^2$ on $[0, 1]$	$-5.2 \times 10^{-3}$	0	$-8.3 \times 10^{-6}$	0
$\frac{1}{x}$ on $[1, 2]$	$-1.9 \times 10^{-3}$	$1.1 \times 10^{-4}$	$-3.1 \times 10^{-6}$	$3.1 \times 10^{-10}$
$\sin x$ on $[0, \pi]$	$5.2 \times 10^{-2}$	$4.6 \times 10^{-3}$	$8.2 \times 10^{-5}$	$1.1 \times 10^{-8}$
$e^x$ on $[0, 1]$	$-4.5 \times 10^{-3}$	$3.7 \times 10^{-5}$	$-7.2 \times 10^{-6}$	$9.5 \times 10^{-11}$

From 10 to 100 strips, Simpson's error for $\frac{1}{x}$ falls from $3.05 \times 10^{-6}$ to $3.12 \times 10^{-10}$ : ten times as many strips, ten thousand times smaller. On a log-log chart that is a slope of $-4$ , against $-2$ for the other two rules.

Rule: Simpson's error

\begin{gathered} S_n - I \approx \tfrac{1}{180}\,\Delta x^4 \\ \times \big(f'''(b) - f'''(a)\big) \end{gathered}

The error shrinks like $\Delta x^4$ . For a cubic, $f'''$ is a constant, so $f'''(b) - f'''(a) = 0$ and Simpson's rule is exact.

A parabola cannot follow a cubic, yet the errors on the two halves of each pair cancel. For $\int_0^2 x^3\,dx = 4$ , two strips give $\frac{1}{3}(0 + 4 \times 1 + 8) = 4$ exactly. For $x^4$ they give 6.666667 against the exact 6.4, off by $\frac{1}{180} \times 48 = 0.2667$ , just as the rule says.

The size of each rule's error for $\int_1^2 \frac{1}{x}\,dx$ . The trapezoid and midpoint lines are parallel, a factor of 2 apart; Simpson's line is twice as steep, until rounding stops it near $10^{-16}$ .

\frac{2M + T}{3}

import math import numpy as np def trapezoid_rule(f, a, b, n): dx = (b - a) / n y = f(a + dx * np.arange(n + 1)) return dx * (y[0] / 2 + np.sum(y[1:-1]) + y[-1] / 2) def midpoint(f, a, b, n): dx = (b - a) / n return dx * np.sum(f(a + dx * (np.arange(n) + 0.5))) def simpson_rule(f, a, b, n): dx = (b - a) / n w = np.ones(n + 1) w[1:-1:2], w[2:-1:2] = 4, 2 return dx / 3 * np.sum(w * f(a + dx * np.arange(n + 1))) recip = lambda x: 1 / x m4, t4, t8 = midpoint(recip, 1, 2, 4), trapezoid_rule(recip, 1, 2, 4), trapezoid_rule(recip, 1, 2, 8) s8 = simpson_rule(recip, 1, 2, 8) print(round((2 * m4 + t4) / 3, 6), round((4 * t8 - t4) / 3, 6), round(s8, 6), f"{s8 - math.log(2):.1e}") examples = [(np.square, 0, 1, 1 / 3), (recip, 1, 2, math.log(2)), (np.sin, 0, math.pi, 2), (np.exp, 0, 1, math.e - 1)] for f, a, b, exact in examples: print([f"{rule(f, a, b, n) - exact:.1e}" for n in [4, 100] for rule in [midpoint, simpson_rule]]) print(f"{simpson_rule(recip, 1, 2, 10) - math.log(2):.2e}", f"{simpson_rule(recip, 1, 2, 100) - math.log(2):.2e}") print(simpson_rule(lambda x: x**3, 0, 2, 2), round(simpson_rule(lambda x: x**4, 0, 2, 2), 6))

⌘/Ctrl + Enter runs the cell

Where the rules break

The error rules assume that $f$ is smooth on the whole of $[a, b]$ , and Simpson's rule needs pairs of strips. Three things break them.

A corner

$e^{|x - 0.37|}$ has a corner at 0.37, and its exact integral over $[0, 1]$ is $(e^{0.37} - 1) + (e^{0.63} - 1) \approx 1.325345$ . A parabola cannot bend sharply, so the pair of strips holding the corner gets an error of order $\Delta x^2$ :

strips	midpoint error	trapezoid error	Simpson error
10	$-1.5 \times 10^{-3}$	$3.2 \times 10^{-3}$	$1.1 \times 10^{-3}$
20	$-5.4 \times 10^{-4}$	$8.8 \times 10^{-4}$	$1.0 \times 10^{-4}$
100	$-5.5 \times 10^{-6}$	$1.1 \times 10^{-5}$	$-3.3 \times 10^{-5}$

With 100 strips Simpson's rule is worse than both of the others. Its error jumps about with the position of the corner inside its pair.

The fix is to split the integral at the corner and apply the rule to each smooth piece. With 100 strips on each piece, the error is $8.1 \times 10^{-12}$ .

Case study 2.6.5 split an option-pricing integral at the corner of its payoff for this reason.

A curve that blows up at an end

$\int_0^1 \frac{1}{\sqrt{x}}\,dx = 2$ is an improper integral (lesson 2.6.3). The trapezoid and Simpson rules both need the height at $x = 0$ , and $\frac{1}{\sqrt{0}}$ does not exist: Python raises ZeroDivisionError.

The midpoint sum never touches the ends, so it runs, but slowly. Its error is $-6.05 \times 10^{-2}$ with 100 strips and $-6.05 \times 10^{-3}$ with 10,000: a hundred times the strips for ten times the accuracy.

Remove the blow-up first, with a substitution such as $x = u^2$ (lesson 2.6.1), which turns this integral into $\int_0^1 2\,du$ . Or use quad, which handles such ends.

An odd number of strips

Simpson's weights need pairs. Applied blindly to 3 strips of $x^2$ on $[0, 1]$ , the pattern gives weights 1, 4, 2, 1 and the answer $\frac{7}{27} \approx 0.259259$ against $\frac{1}{3}$ . No error message warns you.

With data at an even number of points, which makes an odd number of strips, use one of the fixes in the next section.

x^2

import math import numpy as np def trapezoid_rule(f, a, b, n): dx = (b - a) / n y = f(a + dx * np.arange(n + 1)) return dx * (y[0] / 2 + np.sum(y[1:-1]) + y[-1] / 2) def midpoint(f, a, b, n): dx = (b - a) / n return dx * np.sum(f(a + dx * (np.arange(n) + 0.5))) def simpson_rule(f, a, b, n): dx = (b - a) / n w = np.ones(n + 1) w[1:-1:2], w[2:-1:2] = 4, 2 return dx / 3 * np.sum(w * f(a + dx * np.arange(n + 1))) c = 0.37 corner = lambda x: np.exp(np.abs(x - c)) exact = (math.exp(c) - 1) + (math.exp(1 - c) - 1) for n in [10, 20, 100]: print(n, *[f"{rule(corner, 0, 1, n) - exact:.1e}" for rule in [midpoint, trapezoid_rule, simpson_rule]]) print(f"{simpson_rule(corner, 0, c, 100) + simpson_rule(corner, c, 1, 100) - exact:.1e}") try: print(1 / math.sqrt(0)) except ZeroDivisionError as err: print("ZeroDivisionError:", err) for n in [100, 10_000]: print(n, f"{midpoint(lambda x: 1 / np.sqrt(x), 0, 1, n) - 2:.2e}") dx = 1 / 3 print(round(dx / 3 * (0**2 + 4 * (1 / 3)**2 + 2 * (2 / 3)**2 + 1**2), 6))

⌘/Ctrl + Enter runs the cell

SciPy's `trapezoid` and `simpson`

scipy.integrate has both rules. They take the heights, not the function: trapezoid(y, x=x) and simpson(y, x=x), where x holds the points and y the values there. That suits data you already have, such as a table measured at fixed times.

The points need not be evenly spaced. simpson also accepts an odd number of strips: it treats the last strip with a special formula, so its answer stays accurate.

For a function you can evaluate anywhere, quad (lesson 2.5.3) is usually the better tool: it picks its own points and reports an error estimate.

\frac{1}{x}

import math import numpy as np from scipy.integrate import quad, simpson, trapezoid x = np.linspace(1, 2, 11) y = 1 / x print(round(trapezoid(y, x=x), 6), round(simpson(y, x=x), 10)) x = np.linspace(1, 2, 10) print(f"{simpson(1 / x, x=x) - math.log(2):.1e}") print(quad(lambda t: 1 / t, 1, 2)[0], math.log(2))

⌘/Ctrl + Enter runs the cell

Common mistakes

The trapezoid rule halves the two end heights. Counting them in full adds $\frac{1}{2}\Delta x\,\big(f(a) + f(b)\big)$ too much.
Simpson's weights are 1, 4, 2, …, 4, 1, times $\frac{\Delta x}{3}$ , where $\Delta x$ is the width of one strip, not of a pair.
Simpson's rule needs an even number of strips. The pattern applied to an odd number gives a wrong answer with no warning.
$\frac{M + T}{2}$ is not Simpson's rule. The midpoint error is half the trapezoid's, so the weights are 2 and 1.
A corner or a blow-up spoils the error rules. Split at a corner; remove a blow-up at an end, or use quad.
trapezoid and simpson take values, not a function. Evaluate $f$ at the points first.

Trapezoid and Simpson rules

0 of 11 correct

1.
Use the trapezoid rule with 2 strips for $x^2$ on $[0, 2]$
2.
For $\frac{1}{x}$ on $[1, 3]$ with 2 strips, the left sum is $\frac{3}{2}$ and the right sum $\frac{5}{6}$ . What is the trapezoid rule?
3.
Use the rule $\frac{1}{12}\Delta x^2\,\big(f'(b) - f'(a)\big)$ to find the trapezoid error for $x^2$ on $[0, 3]$ with $\Delta x = 0.5$
4.
With some strips, a trapezoid sum is too large by 0.006. About how far off is the midpoint sum with the same strips? Include the sign
5.
A trapezoid error is 0.004 with 10 strips. About what is it with 20 strips?
6.
Use Simpson's rule with 2 strips for $\frac{1}{x}$ on $[1, 3]$
7.
For $x^2$ on $[0, 1]$ with 1 strip, $M_1 = \frac{1}{4}$ and $T_1 = \frac{1}{2}$ . Find $S_2 = \frac{2M_1 + T_1}{3}$
8.
By what factor does Simpson's error shrink when the number of strips is multiplied by 10?
9.
Use Simpson's rule with 2 strips for $x^3 - x$ on $[0, 2]$ . Simpson is exact for cubics
10.
The Simpson pattern applied blindly to 3 strips of $x^2$ on $[0, 1]$ uses weights 1, 4, 2, 1 times $\frac{1}{9}$ . What does it give?
11.
Find the midpoint sum for $\frac{1}{\sqrt{x}}$ on $[0, 1]$ with 1 strip (type sqrt for a square root)

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

T_n = \Delta x\big[\frac{1}{2}f(x_0) + f(x_1) + \cdots + f(x_{n-1}) + \frac{1}{2}f(x_n)\big] = \frac{1}{2}(L_n + R_n)

Trapezoid and Simpson rules

The trapezoid rule

Rule: Trapezoid rule

The trapezoid error, and the midpoint's

Rule: Trapezoid and midpoint errors

Simpson's rule

Rule: Simpson's rule

Simpson's rule from the other two

Rule: Simpson from midpoint and trapezoid

Its error shrinks like $\Delta x^4$

Rule: Simpson's error

Where the rules break

A corner

A curve that blows up at an end

An odd number of strips

SciPy's `trapezoid` and `simpson`

Common mistakes

Trapezoid and Simpson rules

Summary card

Free further reading

The trapezoid rule

Rule: Trapezoid rule

The trapezoid error, and the midpoint's

Rule: Trapezoid and midpoint errors

Simpson's rule

Rule: Simpson's rule

Simpson's rule from the other two

Rule: Simpson from midpoint and trapezoid

Its error shrinks like Δx4\Delta x^4Δx4

Rule: Simpson's error

Where the rules break

A corner

A curve that blows up at an end

An odd number of strips

SciPy's trapezoid and simpson

Common mistakes

Summary card

Free further reading

Its error shrinks like $\Delta x^4$

SciPy's `trapezoid` and `simpson`