Controls
Approximating definite integrals with quadrature rules
A quadrature rule approximates a definite integral by a weighted sum of function values at a finite number of sample points ("nodes"): $$\int_a^b f(x)\,dx \;\approx\; \sum_{i=1}^{n} w_i \, f(x_i).$$ All rules below are stated on the reference interval $u \in [-1,1]$. An integral over an arbitrary interval $[a,b]$ is mapped onto it with the linear substitution $$\int_a^b f(x)\,dx = \frac{b-a}{2} \int_{-1}^{1} f\!\left(\frac{b-a}{2}u + \frac{a+b}{2}\right) du,$$ so any rule for $[-1,1]$ immediately gives a rule for $[a,b]$.
Standard (Newton–Cotes) rules
All Newton–Cotes rules follow the same recipe: fix the nodes to be equally spaced in advance, replace $f$ on each subinterval by the polynomial that interpolates it at those nodes, and integrate that polynomial exactly — its integral is a fixed linear combination of the sampled values, which is exactly where the weights $w_i$ come from. Because the nodes are fixed rather than optimized, this only pins down $n$ degrees of freedom (the weights), not the $2n$ that Gauss–Legendre gets to use below. With $n$ subintervals of width $h = 2/n$ on $[-1,1]$:
- Midpoint rule replaces $f$ on each subinterval by the constant $f(x_i)$ sampled at its centre, so $w_i = h$ for all $i$. Concretely, with centres $x_i = -1 + h\big(i+\frac12\big)$, $i=0,\dots,n-1$, the integral is approximated by $$\int_{-1}^{1}\! f(u)\,du \;\approx\; h\sum_{i=0}^{n-1} f(x_i) \;=\; h\big[f(x_0)+f(x_1)+\cdots+f(x_{n-1})\big].$$ Left plot: every rectangle has the same width $h$ and height $f(x_i)$, and this sum is exactly the total shaded area. A Taylor expansion of $f$ around the centre $x_i$ gives the local (single-panel) error $$\int_{x_i-h/2}^{x_i+h/2}\! f(x)\,dx \;-\; h\,f(x_i) \;=\; \frac{h^3}{24}f''(\xi_i)$$ for some $\xi_i$ in the panel; summing the $n$ panels leaves a global error of order $O(h^2)$, explicitly bounded by $\frac{b-a}{24}h^2 \max|f''|$.
- Trapezoidal rule replaces $f$ on each subinterval by the line through its two endpoint samples, so consecutive panels share a node: $w_i = h$ except at the two ends of the composite rule, where $w_0 = w_n = h/2$. With boundary samples $x_i = -1+ih$, $i=0,\dots,n$, this gives $$\int_{-1}^{1}\! f(u)\,du \;\approx\; h\left[\frac{f(x_0)}{2} + f(x_1) + \cdots + f(x_{n-1}) + \frac{f(x_n)}{2}\right],$$ i.e. every interior sample counts once, and the two endpoint samples — shared by only one panel each — count half as much. Its local error $$\int_{x_i}^{x_i+h}\! f(x)\,dx \;-\; \frac{h}{2}\big(f(x_i)+f(x_i{+}h)\big) \;=\; -\frac{h^3}{12}f''(\xi_i)$$ has the same order $O(h^2)$ as the midpoint rule — globally $\frac{b-a}{12}h^2\max|f''|$ — but with twice the error constant and the opposite sign. (That the two errors are related this simply is exactly what Romberg-style extrapolation methods, not covered here, exploit to cancel the leading term and jump to a higher order for free.)
- Simpson's rule replaces $f$ on each pair of subintervals ($n$ even) by the parabola through its 3 samples. With the same boundary samples $x_i = -1+ih$, $i=0,\dots,n$, as the trapezoidal rule, integrating each parabola and summing gives the alternating-weight formula $$\int_{-1}^{1}\! f(u)\,du \;\approx\; \frac{h}{3}\Big[f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + \cdots + 4f(x_{n-1}) + f(x_n)\Big],$$ i.e. weights $\frac{h}{3},\frac{4h}{3},\frac{2h}{3},\frac{4h}{3},\dots,\frac{h}{3}$: endpoints count once, odd-indexed (parabola-centre) samples count $4\times$, even-indexed (shared parabola-boundary) samples count $2\times$. A parabola has one more free coefficient than a line, yet still only has to match its 3 sample points exactly — the same "one interpolation condition per coefficient" balance as the midpoint and trapezoidal rules — but this happens to buy two full polynomial degrees of exactness at once (Simpson's rule integrates cubics exactly, not just parabolas, because the leading cubic term is odd and cancels by symmetry around the midpoint). The local error over one pair of panels is $$-\frac{h^5}{90}f^{(4)}(\xi_i),$$ giving a global error of order $O(h^4)$, i.e. $\frac{b-a}{180}h^4\max|f^{(4)}|$ — markedly better than the other two for the same $h$, but the nodes are still equally spaced and fixed a priori, unlike Gauss–Legendre below.
All three formulas above assume $f$ is smooth enough for the relevant derivative ($f''$ or $f^{(4)}$) to exist and stay bounded; the "Non-smooth: $|x|$" test function violates this at $u=0$; try it with each rule and note how the error columns in the results table drop to the common, much slower, first-order rate regardless of which rule is selected.
Gauss–Legendre quadrature
Gauss–Legendre quadrature keeps the same basic shape as the rules above — a weighted sum of $n$ samples, $$\int_{-1}^{1}\! f(u)\,du \;\approx\; \sum_{a=1}^{n} w_a\, f(s_a) \;=\; w_1 f(s_1) + w_2 f(s_2) + \cdots + w_n f(s_n),$$ but no longer fixes the nodes $s_a$ in advance. Instead it treats both the $n$ nodes $s_a$ and the $n$ weights $w_a$ as free parameters — $2n$ degrees of freedom in total — and chooses them so that the rule integrates every polynomial of degree $\le 2n-1$ exactly. This is twice the polynomial degree a Newton–Cotes rule can reach with the same number of samples. The optimal nodes turn out to be the roots of the $n$-th Legendre polynomial $P_n$, and the corresponding weights follow directly from them: $$P_n(s_a) = 0, \qquad w_a = \frac{2}{\left(1-s_a^2\right)\left[P_n'(s_a)\right]^2}.$$ The Legendre polynomials satisfy the recurrence $P_0=1$, $P_1(u)=u$, $k P_k(u) = (2k-1)\,u\,P_{k-1}(u) - (k-1)P_{k-2}(u)$, which is used together with Newton's method (starting from $s_a^{(0)} = \cos\!\big(\pi\,\frac{a-0.25}{n+0.5}\big)$) to compute the $s_a$ in the code below — the same approach as the classic gauleg routine [PTVF07].
Because the weights are no longer all equal, the right plot draws the $a$-th term $w_a f(s_a)$ of the sum as a rectangle of height $\psi(s_a)$ centred at the node $s_a$, whose width is the weight $w_a$ itself (widths sum to $2$, the length of $[-1,1]$, just like the subinterval widths on the left). A wide rectangle marks a node that "represents" a larger share of the interval; nodes cluster more densely near the ends of $[-1,1]$, where the Legendre polynomial's roots naturally accumulate, which is what keeps their individual weights roughly balanced.
Why Gauss–Legendre usually wins
For a smooth (analytic) integrand, Gauss–Legendre quadrature converges spectrally — the error shrinks faster than any polynomial power of $1/n$ — while the Newton–Cotes rules above are stuck at their fixed algebraic order ($O(h^2)$ or $O(h^4)$). Increase $n$ for the smooth test functions below and watch the reported error columns: the Gauss–Legendre error collapses to machine precision within a handful of points, long before the other rules catch up. This advantage disappears for the non-smooth test function $f(x)=|x|$, whose derivative jumps at $x=0$: there, Gauss–Legendre still converges, but only at the same modest algebraic rate as the other rules, because the argument that guarantees spectral convergence relies on $f$ being (arbitrarily often) differentiable.
References
- [PTVF07] W. Press, S. Teukolsky, W. Vetterling, B. Flannery. Numerical Recipes: The Art of Scientific Computing, 3rd ed., Cambridge University Press, 2007. (Sect. 4.6, routine gauleg)
- [GW69] G. Golub, J. Welsch. Calculation of Gauss Quadrature Rules. Mathematics of Computation, 1969, 23, 221–230.
- [SB02] J. Stoer, R. Bulirsch. Introduction to Numerical Analysis, 3rd ed., Springer, 2002.
- [AS64] M. Abramowitz, I. Stegun (eds.). Handbook of Mathematical Functions. National Bureau of Standards, 1964. (Tables of Gauss–Legendre nodes/weights)