When there is no formula: trapezoids, Simpson’s parabolas, the race between rules, and throwing darts at random
Most integrals have no formula, and the ones that do are often easier to compute without it. A computer integrates by sampling the function at a handful of points and adding up a weighted sum. Which points, which weights, and how fast the answer improves as the samples multiply: that is the whole subject, and the differences between rules are dramatic.
Pick a rule and see the shape it really integrates. Raise n and compare the errors; the dashed arcs are Simpson’s parabolas, one through each three samples.
Simpson’s rule needs an even number of panels (they are used in pairs), so an odd n is rounded up. Pink marks area below the axis, which counts as negative. The polynomial preset is a cubic, which Simpson integrates exactly: its parabolas miss the curve, but the misses cancel.
With n panels of width h = (b − a)/n, the classic error estimates for a smooth f are
for some ξ in [a, b] (for the left sum, f′ is an average value). The midpoint and trapezoid errors have opposite signs and the midpoint’s is half the size, so the blend (2M + T)/3 cancels the h² terms. That blend is Simpson’s rule. It integrates every cubic exactly, one degree more than the parabolas it is built from, because the cubic’s misfit is odd about each panel pair’s centre and cancels.
The rule carries Thomas Simpson’s name from his 1743 book, but it is older. Bonaventura Cavalieri and James Gregory used it in the seventeenth century, and Johannes Kepler’s 1615 rule for gauging the volume of wine barrels is the same formula applied to a solid.
Run the race as n doubles from 2 to 1024. On these axes an error proportional to n⁻ᵖ is a straight line of slope −p. Then try the rough and the periodic functions.
The fitted order is minus the least-squares slope of log |error| against log n, using only the points shown that are above 10⁻¹³; it is blank when fewer than two points qualify. Hollow markers have hit exactly zero error in floating point and are drawn at 10⁻¹⁶. The dashed guides have slopes −1, −2 and −4.
A rule has order p if its error behaves like C·n⁻ᵖ, so that doubling n divides the error by 2ᵖ: by 2 for the left sum, 4 for the midpoint and trapezoid rules, and 16 for Simpson. On log–log axes the order is the slope. Every rule eventually meets the floor of floating-point arithmetic: numbers carry about 16 significant digits, and errors near 10⁻¹⁶ are rounding noise.
The orders are promises with conditions. Simpson’s error estimate needs a bounded fourth derivative. For √x the first derivative is already infinite at 0, and the error is dominated by the first panel, where every rule except the left sum falls to order about 1.5. A corner, as in |x − 0.3|, holds Simpson to order 2. Before the panels are fine enough to see a function’s features, as with the sharp peak of 1/(1 + 25x²), the rates do not apply at all.
Periodic functions break the rules the other way. The Euler–Maclaurin formula expands the trapezoid rule’s error in powers h², h⁴, h⁶, …, with coefficients that are differences of odd derivatives at the two ends. Over a full period those differences are all zero, so the error falls faster than any power of n, exponentially for analytic functions. Trefethen and Weideman surveyed this “exponentially convergent trapezoidal rule” in 2014. The same expansion powers Romberg integration, which repeatedly cancels the leading error term.
Throw random darts at a box and count the fraction under the curve. Watch the running estimate wander inside a funnel that narrows only like 1/√N.
Darts come from a seeded generator (Mulberry32, seed 1), so “Start again” replays exactly the same throws. The run stops at 200,000 darts. Green darts are under the curve, red ones above it. The standard error is (box area)·√(p(1 − p)/N), where p is the true fraction of the box under the curve.
With N darts and a true hit fraction p, the number of hits is binomial, so the estimate (box area)·(hits/N) has standard error (box area)·√(p(1 − p)/N). Halving the error takes four times as many darts; one more correct digit takes a hundred times as many. In one dimension that is hopeless next to Simpson’s rule, as the comparison in the demo shows.
The point is that the 1/√N rate does not depend on the dimension. A grid rule with n points per axis needs nᵈ points in d dimensions, so Simpson’s error, h⁴, becomes N^(−4/d) in terms of the total number of points N, worse than Monte Carlo once d is above 8. Integrals over hundreds of dimensions, in physics, statistics and finance, are done by random sampling or its refinements.
The method was named and developed at Los Alamos in the late 1940s by Stanisław Ulam, John von Neumann and Nicholas Metropolis; Metropolis and Ulam’s paper “The Monte Carlo Method” appeared in 1949.
Spot a mistake or have an idea?