Numerical Quadrature

A definite integral is a linear functional: hand it a function and it returns a number. Numerical quadrature approximates that functional by something a computer can actually evaluate — a finite weighted sum of the integrand sampled at a handful of points, $\int_a^b f(x)\,\mathrm{d}x \approx \sum_i w_i f(x_i)$ . Everything interesting lives in the two choices this formula exposes: where to place the nodes $x_i$ , and how heavily to weight each sample $w_i$ . Fix the nodes on an even grid and integrate the interpolating polynomial and you get the Newton–Cotes family — the trapezoidal rule, Simpson's rule, and their cousins. Let the nodes move too and demand the highest-degree polynomial be integrated exactly, and the optimal positions turn out to be the roots of an orthogonal polynomial: Gaussian quadrature, which doubles the accuracy for the same number of samples. Two threads run through the story. The first is degree of exactness — the largest polynomial degree a rule integrates without error — because Taylor's theorem guarantees a rule exact for high-degree polynomials is accurate on any smooth integrand. The second is error structure: Euler–Maclaurin shows the trapezoidal error is a clean series in even powers of the step size, and that one fact powers both Romberg extrapolation and the spectral accuracy of the trapezoidal rule on periodic integrands. For a physicist, quadrature is the benchmark that validates every asymptotic method (Laplace, stationary phase, steepest descent) and the engine that evaluates WKB phase integrals and Bohr–Sommerfeld actions when no closed form exists.

Key takeaways

Numerical quadrature replaces a definite integral $\int_a^b f\,\mathrm{d}x$ by a weighted sum $\sum_i w_i f(x_i)$ of integrand values at chosen nodes; the whole design problem is choosing the $x_i$ and $w_i$ .
Newton–Cotes rules fix equally spaced nodes and integrate the Lagrange interpolant; the trapezoidal rule is second order, $E=-(b-a)h^2 f''/12$ , and Simpson's rule is fourth order, $E=-(b-a)h^4 f^{(4)}/180$ .
Gaussian quadrature frees the nodes too, reaching degree of exactness $2n-1$ with only $n$ points by placing them at the roots of orthogonal polynomials — for $n=2$ the nodes are $\pm 1/\sqrt3$ with weights $1$ .
On the benchmark $\int_0^1 e^x\,\mathrm{d}x=1.71828183$ , four Gauss nodes give an error of $9\times10^{-10}$ , while the trapezoidal rule needs thousands of panels to match.
Composite rules, Romberg extrapolation of the trapezoidal rule, and adaptive subdivision are the practical tools; for smooth periodic integrands the trapezoidal rule alone is spectrally accurate.
Key idea: replace the integral — a linear functional — by a weighted sum of samples whose nodes and weights are tuned to integrate polynomials exactly to the highest possible degree; $n$ nodes reach degree $n-1$ for Newton–Cotes but $2n-1$ for Gauss.

The core insight

Replace the integral, a linear functional, by a weighted sum of integrand values whose nodes and weights are chosen to integrate polynomials exactly up to the highest possible degree; for $n$ nodes Newton–Cotes reaches degree $n-1$ (or $n$ by symmetry) while Gauss reaches $2n-1$ .

Where this lives today

Quadrature is quietly everywhere in computational physics. The `quad` routine in every scientific library — SciPy, GSL, MATLAB — is built on QUADPACK, whose adaptive engine pairs a Gauss rule with a higher-order Gauss–Kronrod rule reusing the same nodes, so their difference gives a nearly free local error estimate that drives where the interval is subdivided. For extreme precision the double-exponential (tanh–sinh) transformation maps endpoint singularities to a smooth, exponentially decaying integrand on which the trapezoidal rule is spectrally accurate — the route to evaluating constants and lattice sums to thousands of digits. In semiclassical work, adaptive Gauss–Legendre quadrature evaluates the WKB phase integral and the Bohr–Sommerfeld action $\oint p\,\mathrm{d}q$ between turning points, a Gauss–Chebyshev rule taming the square-root singularity of $\sqrt{E-V(q)}$ . High-dimensional integrals — partition functions, Feynman loop integrals, lattice-QCD path integrals — are the open frontier: tensor-product Gauss rules and Smolyak sparse grids reach moderate dimension before the curse of dimensionality forces the switch to Monte Carlo. And because a good quadrature answer is essentially exact, it remains the gold standard against which Laplace's method, Watson's lemma, and stationary phase are validated in their large-parameter limits.

The physical picture

A quadrature rule is a recipe for weighing an integrand. You may look at $f$ only at a few sample points; each sample gets a weight saying how much "area" that reading is worth, and the weighted total estimates the integral. The trapezoidal rule fits a straight chord across each panel and adds the trapezoid areas — cheap, robust, second order. Simpson fits a little parabola through each triple and does dramatically better for the same spacing. The more faithfully the interpolant tracks the curve, the thinner the slivers you miss — and those slivers are the error.

Figure

The composite trapezoidal rule fits a straight chord across each panel of width h and adds the trapezoid areas; interior nodes carry weight h while the two endpoints carry h/2, and the slivers between chord and curve give the error −(b−a)h²f″(ξ)/12

Gaussian quadrature is the realization that equally spaced nodes throw away half your freedom. Fixing $n$ nodes on a grid leaves only the $n$ weights to tune — enough for degree $n-1$ . Let the nodes move too and you hold $2n$ free numbers, enough for degree $2n-1$ . The optimal positions are not arbitrary: they are exactly the roots of the degree- $n$ orthogonal polynomial for the interval's weight function — a deep bridge between integration and orthogonality. The payoff is enormous, as the second figure shows on one smooth benchmark: spend more evaluations and the trapezoidal rule crawls down as a power law, Simpson falls faster, and Gauss collapses fastest of all.

Figure

Error versus function evaluations for ∫₀¹ eˣ dx = 1.71828183: the trapezoidal rule falls as O(h²), Simpson as O(h⁴), and Gauss–Legendre collapses fastest — four Gauss nodes already reach an error of 9×10⁻¹⁰, which the trapezoidal rule needs thousands of panels to match

Why does smoothness matter so much? A rule of polynomial degree $p$ has an error proportional to a high derivative of $f$ — trapezoidal on $f''$ , Simpson on $f^{(4)}$ , $n$ -point Gauss on $f^{(2n)}$ . If those derivatives are bounded the error decays like a power of the step size; but a kink, a singularity, or rapid oscillation poisons the high derivatives and stalls convergence. That is the physicist's cue to switch tools: Gauss–Laguerre for a decaying tail, Gauss–Chebyshev for a square-root endpoint singularity, Filon's rule for a fast oscillation — the rule must match the integrand's analytic character, not just its shape.

The core equations

Every quadrature rule is one master formula with different nodes and weights. It approximates the integral, a linear functional, by a finite linear combination of samples, and its quality is the largest polynomial degree it reproduces exactly.

\int_a^b f(x)\,\mathrm{d}x \;\approx\; \sum_{i=1}^{n} w_i\, f(x_i) \qquad (1)

where: $x_i$ — quadrature nodes in $[a,b]$ ; $w_i$ — quadrature weights; $n$ — number of nodes; $f$ — integrand, sampled only at the nodes. Reading it: designing a rule means choosing the $2n$ numbers $\{x_i, w_i\}$ , and its degree of exactness is the largest polynomial degree for which $\approx$ becomes $=$ . If the nodes are fixed, the natural weights come from replacing $f$ by the unique degree- $(n-1)$ polynomial through the samples and integrating it — each weight is the integral of a Lagrange basis function, computed once and reused.

w_i = \int_a^b \prod_{\substack{j=1\\ j\neq i}}^{n} \frac{x-x_j}{x_i-x_j}\,\mathrm{d}x \qquad (2)

where: the product is the $i$ -th Lagrange basis polynomial $\ell_i(x)$ , with $\ell_i(x_j)=\delta_{ij}$ . Reading it: the weights depend only on the node positions, never on $f$ — which is what makes Newton–Cotes tables universal. The workhorse composite rules apply a low-order formula on every panel of a subdivided interval. Trapezoidal weights each interior node equally and halves the endpoints; Simpson alternates $4$ and $2$ from the local parabolic fits and needs an even number of panels.

\int_a^b f\,\mathrm{d}x \approx \frac{h}{2}\Big[f_0 + 2\!\sum_{i=1}^{N-1}\! f_i + f_N\Big],\qquad \int_a^b f\,\mathrm{d}x \approx \frac{h}{3}\Big[f_0 + 4\!\!\sum_{i\,\mathrm{odd}}\!\! f_i + 2\!\!\sum_{i\,\mathrm{even}}\!\! f_i + f_N\Big] \qquad (3)

where: $h=(b-a)/N$ — panel width; $N$ — number of equal panels ( $N$ even for Simpson); $f_i=f(a+ih)$ . Reading it: both cost one function value per node, yet Simpson buys two extra orders from the same samples — the parabola is accidentally exact for cubics. Gauss–Legendre abandons the grid: on $[-1,1]$ it places its $n$ nodes at the roots of the Legendre polynomial $P_n$ and assigns the weights below, for exactness up to degree $2n-1$ .

\int_{-1}^{1} f(x)\,\mathrm{d}x \approx \sum_{i=1}^{n} w_i\, f(x_i),\qquad P_n(x_i)=0,\qquad w_i=\frac{2}{(1-x_i^2)\,[P_n'(x_i)]^2} \qquad (4)

where: $P_n$ — Legendre polynomial of degree $n$ ; $x_i$ — its $n$ roots, all inside $(-1,1)$ ; $w_i$ — the Gauss weights, all strictly positive. Reading it: weight positivity is a theorem, not luck, and it is what keeps the rule stable at high order where Newton–Cotes fails. The trapezoidal error is not a shapeless remainder — the Euler–Maclaurin formula exhibits it as a series in even powers of $h$ with coefficients built from boundary derivatives and Bernoulli numbers. This clean structure is the hinge of the whole subject.

T(h) - \int_a^b f\,\mathrm{d}x = \sum_{k=1}^{\infty}\frac{B_{2k}}{(2k)!}\,h^{2k}\big[f^{(2k-1)}(b)-f^{(2k-1)}(a)\big] \qquad (5)

where: $T(h)$ — trapezoidal estimate at step $h$ ; $B_{2k}$ — Bernoulli numbers; the bracket is a boundary-derivative difference. Reading it: only even powers appear, so extrapolation cancels them two orders at a time — and for a smooth periodic integrand every boundary difference vanishes, giving spectral accuracy. Romberg integration turns that even-power structure into an algorithm: a triangular table whose first column is trapezoidal estimates on halved grids, then Richardson-extrapolated across columns, each step killing the next even power of $h$ .

R_{j,k} = \frac{4^{k}R_{j,k-1}-R_{j-1,k-1}}{4^{k}-1} \qquad (6)

where: $R_{j,0}$ — trapezoidal estimate with $2^{j}$ panels; $R_{j,k}$ — the $O(h^{2k+2})$ -accurate extrapolant; $4^{k}-1$ — the Richardson denominator cancelling the $h^{2k}$ term. Reading it: the first extrapolation ( $k=1$ ) is algebraically identical to composite Simpson, and each further column converges faster for smooth $f$ . (The Gauss error bound quantifying the degree- $2n-1$ advantage — its error rides on $f^{(2n)}$ with a factorially shrinking prefactor — is the last row of the limiting-cases table below.)

Derivation AThe interpolatory weights from a single integration

The goal is to show that once the nodes are fixed there is a canonical, integrand-independent choice of weights that makes the rule exact for every polynomial of degree $n-1$ , and that this choice is nothing more than the integral of each Lagrange basis function — so Newton–Cotes weights are computed once and stored, never recomputed per problem. We replace $f$ by its interpolant, integrate, and read off the weights.

Step 1. Interpolate: the unique polynomial of degree $\le n-1$ agreeing with $f$ at the nodes is $p(x)=\sum_{i=1}^n f(x_i)\,\ell_i(x)$ , with $\ell_i(x)=\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}$ and $\ell_i(x_j)=\delta_{ij}$ . Why: replacing $f$ by $p$ is the only step that introduces any error — and none at all when $f$ is a polynomial of degree $\le n-1$ .

Step 2. Integrate term by term: by linearity $\int_a^b f\,\mathrm{d}x \approx \sum_i f(x_i)\int_a^b \ell_i(x)\,\mathrm{d}x$ . Why: the constants $f(x_i)$ slide outside, leaving only the fixed basis functions to integrate.

Step 3. Read off the weights $w_i=\int_a^b \ell_i(x)\,\mathrm{d}x$ , giving equation (1). Why: they depend only on the node placement, so for a fixed grid they are universal — the trapezoidal $(h/2,h,\dots,h/2)$ and Simpson $(h/3)(1,4,2,\dots,1)$ patterns are just these integrals.

Step 4. Check exactness: if $f$ has degree $\le n-1$ then $p=f$ and the rule is exact. Why: this fixes the degree of exactness at (at least) $n-1$ , which node symmetry can lift by one, as Simpson shows. *Verified: for $n=2$ this gives $w_1=w_2=(b-a)/2$ (trapezoidal); for three equally spaced nodes, $(h/3)(1,4,1)$ (basic Simpson).*

Go deeper · Derivation B · 1 minSimpson's rule and its fourth-order error

Derivation BSimpson's rule and its fourth-order error

The aim is to derive the three-point Simpson formula by integrating a quadratic interpolant across two panels, and then to explain the pleasant surprise that its error is one order higher than the quadratic fit would suggest — the rule is exact for cubics too. We centre the algebra at the midpoint, integrate the local Taylor polynomial, watch the odd-power terms die by symmetry, and identify the surviving quartic term as the leading error.

Step 1. Place symmetric nodes $x_0 = x_1 - h$ , $x_2 = x_1 + h$ about $x_1$ and Taylor-expand about $x_1$ : $f=f_1+f_1'(x-x_1)+\tfrac12 f_1''(x-x_1)^2+\tfrac16 f_1'''(x-x_1)^3+\tfrac1{24}f_1^{(4)}(x-x_1)^4+\cdots$ . Why: centring makes $(x-x_1)$ odd over the symmetric interval, so half the terms vanish on integration.

Step 2. Integrate across $[x_0, x_2]$ : every odd power integrates to zero, leaving $\int_{x_0}^{x_2} f\,\mathrm{d}x = 2h f_1 + \tfrac{h^3}{3}f_1'' + \tfrac{h^5}{60}f_1^{(4)}+\cdots$ . Why: the vanishing of the odd terms, cubic included, is exactly why Simpson turns out exact for cubics.

Step 3. Eliminate $f_1''$ with the central second difference $f_1''=(f_0-2f_1+f_2)/h^2 - \tfrac{h^2}{12}f_1^{(4)}+\cdots$ , collapsing $2hf_1+\tfrac{h^3}{3}f_1''$ to $\tfrac{h}{3}(f_0+4f_1+f_2)$ . Why: this trades the unknown derivative for sampled values, giving the $1,4,1$ pattern plus its own $O(h^4)$ correction.

Step 4. Collect the leftover: the surviving quartic terms combine to $E = \tfrac{h^5}{60}f^{(4)} - \tfrac{h^5}{36}f^{(4)} = -\tfrac{h^5}{90}f^{(4)}(\xi)$ for the basic rule. Why: the local error is $O(h^5)$ , a full order beyond the naive quadratic expectation.

Step 5. Composite it: summing $N/2$ blocks multiplies the per-block error by $\sim(b-a)/(2h)$ , giving $E=-\dfrac{(b-a)h^4}{180}f^{(4)}(\xi)$ . Why: fourth-order accuracy for zero extra samples — why Simpson is the default low-order rule. *Verified: $-h^5/90$ basic and $-(b-a)h^4/180$ composite match Arfken–Weber Ch. 10 and Numerical Recipes §4.1; on $\int_0^1 e^x\,\mathrm{d}x$ Simpson with $N=4$ gives $1.71831884$ (error $3.7\times10^{-5}$ ).*

Go deeper · Derivation C · 1 minWhy the optimal Gauss nodes are roots of an orthogonal polynomial

Derivation CWhy the optimal Gauss nodes are roots of an orthogonal polynomial

The goal is to prove the central theorem of Gaussian quadrature: an $n$ -point rule is exact for all polynomials up to degree $2n-1$ — twice what a fixed grid allows — precisely when the nodes are the zeros of the degree- $n$ orthogonal polynomial for the interval. We divide an arbitrary high-degree polynomial against $P_n$ , use orthogonality to annihilate one piece under the integral, and choose the nodes so the same piece vanishes in the sum, leaving a low-degree remainder the interpolatory weights integrate exactly.

Step 1. Divide: any $p(x)$ of degree $\le 2n-1$ writes as $p=q\,P_n+r$ with $q,r$ of degree $\le n-1$ , where $P_n$ is the Legendre polynomial. Why: division by a degree- $n$ divisor leaves a remainder below degree $n$ .

Step 2. Kill the quotient under the integral: since $P_n$ is orthogonal to every lower-degree polynomial, $\int_{-1}^1 q\,P_n\,\mathrm{d}x = 0$ , so $\int_{-1}^1 p\,\mathrm{d}x = \int_{-1}^1 r\,\mathrm{d}x$ . Why: orthogonality says the quotient piece contributes nothing to the true integral.

Step 3. Kill it in the sum too: take the nodes $x_i$ as the roots of $P_n$ , so $P_n(x_i)=0$ and $p(x_i)=r(x_i)$ . Why: at the zeros of $P_n$ the quotient term drops out pointwise.

Step 4. Close the loop: the $n$ interpolatory weights integrate the degree- $(n-1)$ remainder exactly, so $\sum_i w_i p(x_i)=\int_{-1}^1 r\,\mathrm{d}x=\int_{-1}^1 p\,\mathrm{d}x$ — exact to degree $2n-1$ . Why: this doubles exactness for free and forces the weights positive, the source of Gauss stability. *Verified: for $n=2$ , $P_2(x)=\tfrac12(3x^2-1)=0$ gives $x=\pm 1/\sqrt3=\pm 0.5773502692$ and $w=1$ , matching Abramowitz & Stegun Table 25.4.*

Limiting cases

How the error of the common rules scales, and which tool each regime calls for:

Case	Result	Meaning
Composite trapezoidal	$E=-\dfrac{(b-a)h^{2}}{12}\,f''(\xi)$	Second order $O(h^2)$ ; cheapest and most robust, and spectrally accurate for smooth periodic integrands.
Composite Simpson	$E=-\dfrac{(b-a)h^{4}}{180}\,f^{(4)}(\xi)$	Fourth order $O(h^4)$ at no extra samples; the default general-purpose low-order rule.
Gauss–Legendre ( $n$ nodes)	$E=\dfrac{(b-a)^{2n+1}(n!)^{4}}{(2n+1)\,[(2n)!]^{3}}\,f^{(2n)}(\xi)$	Exact to degree $2n-1$ ; superb for smooth non-periodic integrands, weights all positive.
Smooth periodic, full period	$E$ decays faster than any $h^{p}$	Euler–Maclaurin boundary terms all cancel — the trapezoidal rule becomes spectrally accurate (the FFT's foundation).
Singular / semi-infinite	$\displaystyle\int_0^\infty e^{-x}g(x)\,\mathrm{d}x\approx\sum_i w_i\,g(x_i)$	Use Gauss–Laguerre (weight $e^{-x}$ ), Gauss–Hermite ( $e^{-x^2}$ ), or Gauss–Chebyshev ( $1/\sqrt{1-x^2}$ ) to absorb the awkward factor analytically.

Common trapHigher-order Newton–Cotes rules are always more accurate

They are not. A single high-order Newton–Cotes formula over a wide interval is a trap: above degree roughly $8$ the weights grow large and alternate in sign, so rounding errors are amplified and the underlying interpolating polynomial oscillates wildly near the endpoints — Runge's phenomenon. Accuracy and stability are both lost. The fix is never to raise the order of a single-interval rule; instead subdivide and apply a low-order composite rule on each panel, or switch to Gaussian quadrature, whose weights stay positive at every order. Two further foot-guns: Simpson's rule needs an even number of panels, and subtracting two nearly equal trapezoidal estimates (as in Romberg) loses precision unless you carry enough digits.

Worked examples

1 · One integral, three rules (fully worked)

Evaluate the benchmark $\int_0^1 e^x\,\mathrm{d}x$ , whose exact value is $e-1 = 1.71828183$ , by composite trapezoidal ( $N=4$ ), composite Simpson ( $N=4$ ), and Gauss–Legendre, and compare the error per function evaluation.

Step 1. Trapezoidal, $N=4$ , $h=0.25$ . Sample $f(x)=e^x$ at $x=0,0.25,0.5,0.75,1$ : $1,\,1.28402542,\,1.64872127,\,2.11700002,\,2.71828183$ . Then $T=\tfrac{h}{2}[f_0+2(f_1+f_2+f_3)+f_4]=0.125\,[1+2(5.04974671)+2.71828183]=\mathbf{1.72722190}$ , an error of $\mathbf{8.9\times10^{-3}}$ from five evaluations.

Step 2. Simpson, $N=4$ , same five samples. Apply the $1,4,2,4,1$ pattern: $S=\tfrac{0.25}{3}[1+4(1.28402542)+2(1.64872127)+4(2.11700002)+2.71828183]=\mathbf{1.71831884}$ , an error of $\mathbf{3.7\times10^{-5}}$ from the same five evaluations — two extra orders for nothing.

Step 3. Two-point Gauss. Map $[0,1]$ to $[-1,1]$ by $x=\tfrac12(1+t)$ , so nodes $t=\pm 1/\sqrt3$ land at $x=0.21132487,\,0.78867513$ . With weights $1$ and Jacobian $\tfrac12$ : $\tfrac12[e^{0.21132487}+e^{0.78867513}]=\mathbf{1.71789638}$ , an error of $\mathbf{3.9\times10^{-4}}$ from just two evaluations — already beating trapezoidal's five.

Step 4. Three-point Gauss. Nodes $t=0,\pm\sqrt{3/5}$ map to $x=0.11270167,\,0.5,\,0.88729833$ with weights $\tfrac59,\tfrac89,\tfrac59$ : $\tfrac12[\tfrac59(e^{0.11270167}+e^{0.88729833})+\tfrac89 e^{0.5}]=\mathbf{1.71828100}$ , error $\mathbf{8.2\times10^{-7}}$ from three evaluations; a fourth node reaches $9\times10^{-10}$ . Verdict: on a smooth integrand Gauss extracts far more accuracy per sample than any fixed-grid rule.

2 · The Romberg table (fully worked)

Build the Romberg triangle for the same integral from trapezoidal estimates on $1,2,4,8$ panels, and watch it converge to $1.71828183$ .

Step 1. First column — trapezoidal estimates $R_{j,0}=T(2^{j}\text{ panels})$ : $R_{0,0}=1.85914091$ , $R_{1,0}=1.75393109$ , $R_{2,0}=1.72722190$ , $R_{3,0}=1.72051859$ . Only second order; the coarsest is off in the second decimal.

Step 2. First extrapolation $R_{j,1}=(4R_{j,0}-R_{j-1,0})/3$ — algebraically Simpson: $R_{1,1}=1.71886115$ , $R_{2,1}=1.71831884$ , $R_{3,1}=1.71828415$ . The $O(h^2)$ term is gone; these are $O(h^4)$ and the error drops a hundredfold.

Step 3. Second extrapolation, $R_{j,2}=(16R_{j,1}-R_{j-1,1})/15$ , killing the $h^4$ term: $R_{2,2}=1.71828269$ , $R_{3,2}=1.71828184$ . Third extrapolation, $R_{3,3}=(64R_{3,2}-R_{2,2})/63=\mathbf{1.71828183}$ .

Step 4. Read the corner — the full triangle is

\begin{array}{llll} 1.85914091 & & & \\ 1.75393109 & 1.71886115 & & \\ 1.72722190 & 1.71831884 & 1.71828269 & \\ 1.72051859 & 1.71828415 & 1.71828184 & 1.71828183 \end{array}

The bottom-right entry $R_{3,3}=1.71828183$ matches $e-1$ to eight decimals, from only nine function evaluations — the trapezoidal column alone would need on the order of $10^4$ panels for the same accuracy. Verified: each column's error falls by the predicted factor ( $\approx 16$ for column one, faster thereafter).

3 · Gauss exactness and its failure mode

Confirm that the two-point Gauss rule integrates a cubic exactly on $[-1,1]$ , then quantify exactly how it fails on the first polynomial beyond its reach, $x^4$ .

Step 1. Cubic test. For $f(x)=x^3+x^2+x+1$ the exact integral is $\int_{-1}^1 f\,\mathrm{d}x=0+\tfrac23+0+2=\tfrac83$ . The rule gives $f(-1/\sqrt3)+f(1/\sqrt3)$ ; the odd powers cancel between the symmetric nodes, leaving $2(\tfrac13)+2=\tfrac83$ . Exact, as the degree- $2n-1=3$ theorem promises.

Step 2. Quartic test. $\int_{-1}^1 x^4\,\mathrm{d}x=\tfrac25=0.4$ , but the rule gives $(1/\sqrt3)^4+(-1/\sqrt3)^4=2\cdot\tfrac19=\tfrac29=0.22222$ , an error of $\mathbf{0.17778}$ .

Step 3. Match the error bound. The Gauss error formula (limiting-cases table) with $n=2$ , $b-a=2$ , and $f^{(4)}=24$ (constant, so the bound is attained) gives $\dfrac{2^{5}(2!)^4}{5\,(4!)^3}\cdot 24=\dfrac{512}{69120}\cdot 24=\mathbf{0.17778}$ . The predicted and actual errors agree to every digit, confirming the error formula and the sharp cutoff at degree $2n-1$ .

Physics of the portrait

Numerical integration is a thread running from Kepler's wine casks to modern lattice field theory, but its decisive idea — that moving the nodes doubles the accuracy — belongs to one person.

Origins · 1 minCarl Friedrich Gauss and how this was discovered

Carl Friedrich Gauss

1777–1855 · optimal quadrature and orthogonal polynomials

Carl Friedrich Gauss, already famous for number theory, least squares, and the orbit of Ceres, turned to numerical integration in an 1814 memoir with a disarmingly simple question: given a fixed budget of $n$ function evaluations, where should you place them? The Newton–Cotes tradition, inherited from Newton and Roger Cotes a century earlier, had always spaced the nodes evenly and solved for the weights, reaching degree of exactness $n-1$ . Gauss saw that the node positions were themselves free parameters, and that spending all $2n$ degrees of freedom optimally lets an $n$ -point rule integrate polynomials up to degree $2n-1$ — a factor-of-two leap that no cleverness with fixed nodes can match. His own derivation used continued fractions and hypergeometric series, but the modern reading is cleaner: the optimal nodes are the roots of the Legendre polynomials, and the rule's exactness follows directly from their orthogonality. Carl Gustav Jacobi made that connection explicit in 1826, and the framework generalizes seamlessly — Laguerre, Hermite, and Chebyshev weight functions each yield their own Gauss rule for the integrals of quantum mechanics, statistical physics, and scattering theory. It is a characteristic Gauss result: a practical recipe resting on a deep structural fact, so complete that two centuries later it remains the default engine inside every serious numerical library.

History

1722

— Cotes's tables of the equal-node interpolatory weights are published (Newton–Cotes)

1743

— Simpson popularizes the three-point parabolic rule

1814

— Gauss introduces optimal-node quadrature of degree 2n−1

1826

— Jacobi recasts Gauss's nodes as the roots of orthogonal polynomials

1955

— Romberg systematizes Richardson extrapolation of the trapezoidal rule

Going further

Texts. Stone & Goldbart, Mathematics for Physics (CUP 2009), Ch. 8 · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 10 · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed., Ch. 27 · Press, Teukolsky, Vetterling & Flannery, Numerical Recipes, 3rd ed., §4 · Davis & Rabinowitz, Methods of Numerical Integration, 2nd ed.

Provenance. The trapezoidal error $-(b-a)h^2 f''/12$ , the basic and composite Simpson errors $-h^5 f^{(4)}/90$ and $-(b-a)h^4 f^{(4)}/180$ , and the two-point Gauss nodes and weights ( $\pm 1/\sqrt3$ , $w=1$ ) are standard and cross-checked against Arfken–Weber Ch. 10, Riley–Hobson–Bence Ch. 27, and Numerical Recipes §4. A verifier should recompute: (1) the $n=2$ Gauss moments to recover $x^2=\tfrac13$ ; (2) the benchmark $\int_0^1 e^x\,\mathrm{d}x=1.71828183$ by trapezoidal ( $N=4\Rightarrow 1.72722190$ ), Simpson ( $N=4\Rightarrow 1.71831884$ ), and three-point Gauss ( $1.71828100$ ); (3) the Romberg corner $R_{3,3}=1.71828183$ ; (4) that two-point Gauss integrates $x^3$ exactly but errs on $x^4$ by $0.17778$ , matching the error bound with $f^{(4)}=24$ .

Prerequisites. Taylor's Theorem and Polynomial Interpolation · Orthogonal Polynomials · Order Symbols and Gauge Functions. Leads to. Laplace's Method for Integrals · Method of Stationary Phase · Bohr–Sommerfeld Quantization · WKB Approximation.

Ask the tutor: Ask me to walk through the Romberg table for any integral of your choice, or to show exactly where an equally spaced rule breaks down and a Gauss rule saves the day.

What is the "degree of exactness" and why does it predict accuracy on non-polynomials?

The degree of exactness is the largest polynomial degree a rule integrates with zero error. It matters far beyond polynomials because of Taylor's theorem: any smooth integrand looks locally like its Taylor polynomial, so a rule that nails polynomials to degree $d$ leaves an error controlled by the first neglected derivative, $f^{(d+1)}$ . That is why the trapezoidal rule (exact to degree $1$ ) errs in $f''$ , Simpson (degree $3$ ) in $f^{(4)}$ , and $n$ -point Gauss (degree $2n-1$ ) in $f^{(2n)}$ . Higher degree of exactness means the error rides on a higher derivative, a smaller effect for a smooth function — provided that derivative stays bounded.

When should I reach for Gauss quadrature instead of Simpson or Romberg?

Use a Gauss rule when the integrand is smooth over the whole interval and you cannot afford many evaluations — each is an expensive simulation, or you need a known smooth kernel to high precision. Its degree- $2n-1$ exactness extracts the most accuracy per sample. Simpson and Romberg shine when samples are cheap and you want a built-in convergence check by halving the step, or when your data already sit on an equally spaced grid (Gauss samples at its own irrational nodes). For a singular or infinite endpoint, prefer the matching Gauss variant — Laguerre, Hermite, or Chebyshev — which folds the awkward weight function into the rule.

Why does the trapezoidal rule become spectacularly accurate for periodic functions?

Because its entire error is a sum of boundary terms. The Euler–Maclaurin formula writes the trapezoidal error as a series in $h^{2k}$ whose coefficients are odd-derivative differences at the two endpoints, $f^{(2k-1)}(b)-f^{(2k-1)}(a)$ . For a smooth function periodic over the interval, the endpoints agree in every derivative, so every one of those differences is exactly zero. The whole power-law error series vanishes and what remains decays faster than any power of $h$ — spectral accuracy. This is why the trapezoidal rule, not Simpson, is the natural partner of Fourier analysis and the FFT.

Is Romberg integration really just Simpson in disguise?

Only the first extrapolation column is. Romberg's first column is trapezoidal estimates on successively halved grids; the first Richardson step, $(4T(h/2)-T(h))/3$ , cancels the leading $h^2$ error and is algebraically identical to composite Simpson. But it does not stop there — because the trapezoidal error contains only even powers of $h$ , each further column cancels the next even power ( $h^4$ , then $h^6$ , and so on), producing rules of ever higher order that Simpson never reaches. Romberg is a whole tower of extrapolations of which Simpson is merely the ground floor, and for smooth integrands the deeper columns converge dramatically faster.

Why must the Gauss weights be positive, and does it matter?

It matters a great deal for stability, and it is provable directly: the weight $w_i$ equals $\int (\ell_i(x))^2\,\mathrm{d}x$ for the Lagrange basis at the Gauss nodes — an integral of a square, hence positive. Positive weights make the rule a genuine weighted average of samples, so rounding errors in the $f(x_i)$ cannot be amplified: the sum of the absolute weights equals $b-a$ . High-order Newton–Cotes rules lose exactly this property — some weights turn negative and large, so cancellation magnifies floating-point noise. Positivity is the quiet reason Gauss rules can be pushed to very high order while equally spaced rules cannot.

How do adaptive quadrature routines decide where to add points?

They estimate the local error and refine where it is largest. The standard trick, used in QUADPACK and hence in most `quad` functions, is to evaluate two rules of different order on the same panel — a Gauss rule and a Gauss–Kronrod rule that reuses the Gauss nodes and adds a few more. Their difference is a cheap, reliable proxy for the panel's error: if it exceeds the tolerance the panel is bisected and the procedure recurses on each half, otherwise the panel is accepted. This concentrates evaluations exactly where the integrand misbehaves — near a spike or a near-singularity — and leaves smooth stretches coarsely sampled, which is why adaptive routines handle real-world integrands so robustly.