The free lunch hiding in the sample points
Every quadrature rule you have met so far approximates an integral as a weighted sum of function values: the integral of f over [a, b] is replaced by w_1 f(x_1) + w_2 f(x_2) + ... + w_n f(x_n). In the Newton-Cotes rules of the previous guide — the trapezoidal rule, Simpson's rule and their cousins — the points x_i were fixed in advance, spaced evenly across the interval, and only the weights w_i were yours to choose. With n evenly spaced points you can match every polynomial up to about degree n-1 (a touch more for the symmetric Simpson family). That felt like the whole game. It was not.
Here is the idea that started a whole subject. Why should the points be fixed? Treat the locations x_1, ..., x_n as unknowns too. Now you have 2n knobs to turn — n weights and n positions — instead of just n. A rule has 2n free parameters, so by a simple counting argument you ought to be able to make it exact for a family of functions with 2n free coefficients: the polynomials of degree up to 2n-1, which have exactly 2n coefficients (1, x, x^2, ..., x^{2n-1}). That is the promise of Gaussian quadrature: n points, chosen optimally, integrate every polynomial of degree 2n-1 or less with zero error. Two points reach degree 3; three points reach degree 5; five points reach degree 9.
Why the optimal points are the roots of a special polynomial
The secret is one beautiful fact: the optimal points are not arbitrary — they are the roots of an orthogonal polynomial. Work on the standard interval [-1, 1] (any [a, b] maps onto it by a linear change of variable). The relevant family is the Legendre polynomials P_0, P_1, P_2, ..., where P_n has degree n, and the family is orthogonal: the integral over [-1, 1] of P_m times P_n is zero whenever m is not equal to n. Place your n sample points exactly at the n roots of P_n — which all lie strictly inside (-1, 1) — and choose the matching weights, and the magic of degree 2n-1 falls out automatically. This is Gauss-Legendre quadrature.
Why orthogonality does the trick is worth seeing, because it is the whole engine. Take any polynomial f of degree up to 2n-1. Divide it by P_n: f = q P_n + r, where the quotient q and remainder r each have degree at most n-1. Integrate. The term involving q P_n vanishes — that is exactly the orthogonality statement, since q (degree < n) is a combination of P_0, ..., P_{n-1}, each orthogonal to P_n. So the true integral of f equals the integral of r alone. Meanwhile, at every sample point x_i we have P_n(x_i) = 0 (the points are its roots!), so f(x_i) = r(x_i): the rule never even sees the q P_n part. Both the exact integral and the rule reduce to the same degree-(n-1) remainder r, which an n-point rule handles exactly. The error is therefore zero. The two facts — orthogonality kills the quotient in the true integral, and the roots kill it in the sum — lock together perfectly.
Gauss-Legendre nodes/weights on [-1, 1]:
n = 1: x = 0 w = 2
n = 2: x = +-0.5773502692 (=+-1/sqrt3)
w = 1, 1
n = 3: x = 0, +-0.7745966692
w = 8/9, 5/9, 5/9
exact for all polynomials up to degree 2n-1.
map [-1,1] -> [a,b]: t = (a+b)/2 + (b-a)/2 * x
integral_a^b f dt ~= (b-a)/2 * sum_i w_i f(t_i)
weights are always positive and sum to 2 (the length of [-1,1]).What it buys you, and what it costs
The payoff is dramatic. For the same number of function evaluations, an n-point Gauss rule typically annihilates the error of an n-point Newton-Cotes rule. A tiny worked example: integrate f(x) = x^4 from -1 to 1, whose exact value is 2/5 = 0.4. The two-point Gauss rule samples at +-1/sqrt(3) with weights 1 and 1, giving 1*(1/sqrt3)^4 + 1*(-1/sqrt3)^4 = 2*(1/9) = 2/9... that is not 0.4, because degree 4 exceeds 2*2-1 = 3, so two points cannot reach it. Bump to three points (exact to degree 5) and the rule returns 0.4 to the last digit. Newton-Cotes would need many more evenly spaced points to match a smooth integrand to the same accuracy.
Two more practical virtues. First, the weights are all positive (unlike high-order Newton-Cotes, whose weights turn negative and amplify round-off catastrophically). Positive weights mean the rule cannot magnify the floating-point errors in your f(x_i) values — it stays numerically well-behaved even at high order. Second, Gaussian rules dodge Runge's phenomenon: because their nodes cluster toward the ends of the interval (rather than spreading evenly), they behave like interpolation on the kind of clustered nodes that tame the wild oscillations that ruin high-degree equispaced rules. High order is actually usable here.
Where Gauss fits among the rules, and its limits
Gaussian quadrature is not always the right reach, and it pays to know its rivals. One catch: the nodes for n points share nothing with the nodes for n+1 points, so if a 5-point estimate is not accurate enough, your 7-point retry must throw away all five old evaluations and start fresh. That is wasteful when each f is expensive. The fix is a nested family. Clenshaw-Curtis quadrature uses Chebyshev points, reuses old samples as you refine, and for many smooth functions runs neck-and-neck with Gauss; Gauss-Kronrod pairs add extra points to a Gauss rule precisely so the cheap-to-compute difference between the two estimates serves as a built-in error gauge.
That built-in error estimate is the seed of the next guide's theme — well, of the next rung's: an adaptive quadrature routine applies a Gauss-Kronrod pair on an interval, reads off the error estimate, and if it is too large, splits the interval in two and recurses — spending sample points densely where f wiggles and sparsely where it is calm. This is the same spirit as the Richardson extrapolation and Romberg ideas from earlier in this rung: don't pick a resolution blindly, let the method measure its own error and refine itself. A modern library integrator (the QUADPACK routines behind most scientific software) is exactly an adaptive Gauss-Kronrod engine.
Now the hard wall, and the rung's recurring lesson. Every number these rules return is an approximation, computed in floating point where 0.1 has no exact binary form and a long sum is never associative — so even the theoretically exact two-point Gauss answer for x^4 above arrives with a few units of round-off in its last digits. More importantly, all of this lives in one dimension. Gaussian quadrature's superpower is degree, and degree buys you almost nothing against the curse of dimensionality: an n-point 1D rule becomes an n^d-point grid in d dimensions, so 20 points per axis is 20^10, about ten trillion evaluations, in ten dimensions. Past a handful of dimensions you must abandon grids entirely and turn to Monte Carlo integration, whose O(1/sqrt(N)) error is painfully slow but — crucially — does not care how many dimensions there are. Gauss owns the smooth low-dimensional integral; randomness owns the high-dimensional one.