You need \int_0^1 e^{-x^2}\,dx. Reach for the fundamental theorem of calculus and you hit a wall: e^{-x^2} has no antiderivative expressible in elementary functions. Neither does \frac{\sin x}{x}, nor \sqrt{1 + \cos^2 x} (the arc length of a sine wave), nor most integrands a scientist actually meets. The clean symbolic answer simply does not exist. And yet the area is a perfectly definite number — the bell curve encloses a real amount of ground. We just have to compute it.

This is quadrature: approximating a definite integral \int_a^b f(x)\,dx by sampling the integrand at a handful of points and combining those samples with clever weights. The name is ancient — "quadrature" meant constructing a square of equal area to a curved region, the very problem that stumped Greek geometers ("squaring the circle"). Today it is the engine inside every physics simulation, every financial model, every renderer that computes how much light reaches a pixel. This page teaches the one central idea — fit a shape you can integrate exactly, then add up the pieces — and shows why one particular choice of shape, the humble parabola, is astonishingly good.

From rectangles to better shapes

You already know the crudest quadrature rule. In Riemann sums we chopped [a, b] into n thin strips of width h = \tfrac{b-a}{n} and stood a rectangle on each. A rectangle approximates f on a strip by a constant — the flattest possible fit. It works, but it is wasteful: on a strip where the curve is visibly sloping, the flat top leaves an obvious triangular sliver of error.

So raise the ambition. Instead of fitting a constant on each panel, fit a straight line through the two endpoints — that is the trapezoid rule. Better still, fit a parabola through three points — that is Simpson's rule. Each step up matches more of the curve's shape, so each panel's error shrinks dramatically. This is the whole ladder of Newton–Cotes rules: sample f at equally spaced points, pass a polynomial of degree k through them (exactly the interpolating polynomial), and integrate that polynomial — which we can do exactly, since polynomials always have clean antiderivatives.

The trapezoid rule: fit a line on each panel

Take a single panel from x_0 to x_1 = x_0 + h. Join the two endpoints (x_0, f_0) and (x_1, f_1) with a straight line. The region under that line is a trapezoid — a rectangle of height f_0 with a triangle perched on top — and school geometry gives its area as width times the average of the two heights:

Now lay n panels side by side, with nodes x_0, x_1, \dots, x_n and f_i = f(x_i). Add the panels. Every interior node is shared by two neighbouring trapezoids, so it gets counted twice (weight 1), while the two endpoints appear only once (weight \tfrac12). Collecting terms gives the composite trapezoid rule:

It is the midpoint's cousin from the Riemann page, but with the ends counted at half weight — and that small correction turns a first-order rule into a second-order one. The error of the composite trapezoid rule behaves like

for some point \xi in [a,b]. Two consequences fall straight out of that h^2. First, if f'' = 0 — that is, f is a straight line — the error is exactly zero: the trapezoid rule integrates linear functions perfectly, which makes sense because a line fits a line with nothing left over. Second, halving the step h multiplies h^2 by \tfrac14, so each doubling of n cuts the error by about a factor of four. Let us watch that happen.

Read the error column top to bottom: each time n doubles the error falls to roughly a quarter of the line above (2^{-2} per doubling). That factor of four is the O(h^2) promise made visible.

Simpson's rule: fit a parabola over pairs of panels

A line matches a curve's height and slope but ignores its bend. Fit a parabola instead and you capture the curvature too. A parabola needs three points, so Simpson's rule works over a pair of panels at a time: pass the unique parabola through (x_0, f_0), (x_1, f_1), (x_2, f_2) and integrate it. A short calculation (integrate the interpolating parabola over [x_0, x_2]) collapses to a beautifully simple weighting — 1, 4, 1:

Chain these two-panel blocks along the whole interval (so n must be even). Interior nodes at odd index are always the middle of a block and pick up weight 4; interior nodes at even index sit where two blocks meet and collect 2 = 1 + 1; the two ends keep weight 1. That is the composite Simpson's rule:

Its error is a whole two orders better than the trapezoid rule:

The h^4 is the headline. Halving h now multiplies the error by \left(\tfrac12\right)^4 = \tfrac{1}{16}: each doubling of n cuts the error by about a factor of sixteen, not four. And here is the pleasant surprise buried in that fourth derivative.

Worked example: Simpson on one block, by hand

Estimate \int_0^1 e^{x}\,dx (true value e - 1 \approx 1.718282) with a single Simpson block, n = 2. The nodes are x_0 = 0, x_1 = \tfrac12, x_2 = 1, so h = \tfrac12 and the heights are

Apply the 1, 4, 1 weighting:

With just three function evaluations the error is only about 6\times10^{-4} — the trapezoid rule needs dozens of panels to match that. Now the same integral in code, both rules side by side so you can watch Simpson pull ahead.

Each doubling of n slices the error to about a sixteenth of the previous row — the O(h^4) rate, plain as day.

Order of accuracy: the race, side by side

"Order" is the exponent p in E \approx C h^p. It decides how fast you buy accuracy by spending function evaluations. Trapezoid has p = 2; Simpson has p = 4. That gap sounds modest until you compound it: to gain one more decimal digit (\times\tfrac{1}{10} error), trapezoid must shrink h by \sqrt{10}\approx 3.2, while Simpson needs only 10^{1/4}\approx 1.8. Over many refinements the two curves diverge wildly. Put them in one table for \int_0^\pi \sin x\,dx = 2 and count the digits.

The ratio columns are the punchline: trapezoid's error shrinks by about 4 per row, Simpson's by about 16. By n = 64 Simpson is already near the floor of double-precision arithmetic while trapezoid still limps along with a handful of correct digits. The chart below plots exactly this race on a log scale.

Beyond Newton–Cotes: smarter points, smarter panels

Newton–Cotes rules are handed their sample points — equally spaced, take it or leave it. But where you sample is a free choice, and spending it wisely buys a huge return. That is the idea of Gaussian quadrature: choose both the nodes and the weights to make the rule exact for the highest-degree polynomial possible. With m cleverly placed points — bunched toward the ends of the interval, at the roots of certain special (Legendre) polynomials — a Gauss rule integrates every polynomial up to degree 2m - 1 exactly, roughly double the degree a Newton–Cotes rule reaches with the same number of samples. For smooth integrands it is spectacularly efficient: a handful of well-chosen points can beat thousands of trapezoids.

The other big idea is adaptive quadrature. Not all of an integrand is equally hard: a function might be nearly flat across most of [a,b] and then spike sharply in one narrow region. Spreading panels uniformly wastes effort on the flat part and starves the spike. Adaptive schemes instead estimate the local error on each subinterval (typically by comparing a coarse and a fine Simpson estimate) and recursively subdivide only where the error is too large, packing panels densely where the action is and leaving them sparse where the curve is tame. This is what a library routine like scipy.integrate.quad actually does under the hood — a Gauss rule wrapped in adaptive error control.