Prev
Unit 8 · Numerical Differentiation and Integration
Gaussian Quadrature
55m100 XP
Next
∫Numerical Differentiation and Integration

Gaussian Quadrature

Free nodes, orthogonal polynomials, Gauss–Legendre rules, change of interval and finite elements

Göttingen, 1814
Gauss is computing perturbations of planetary orbits — integrals that must be evaluated again and again, each function value a long hand calculation. Every saved evaluation is hours of work. Newton–Cotes rules fix the nodes at equal spacing and choose only the weights: nn free parameters, so they can match polynomials of degree n−1n - 1 (one more by symmetry). Gauss sees 2n2n parameters — nn weights and nn nodes — and asks whether all of them can be used. In a paper presented to the Göttingen Society in 1814, Methodus nova integralium valores per approximationem inveniendi, he shows that they can: nn well-placed points integrate every polynomial of degree up to 2n−12n - 1 exactly. Twelve years later Jacobi explained why the magic points are what they are: the roots of orthogonal polynomials. Two centuries later, every finite-element program on Earth evaluates its element matrices at Gauss's points.
💡Gauss–Legendre quadrature
On [−1,1][-1, 1], the nn-point Gauss–Legendre rule ∫−11f(x) dx≈∑i=1nwi f(xi)\int_{-1}^{1}f(x)\,dx \approx \sum_{i=1}^{n}w_i\,f(x_i) takes as nodes the nn roots of the Legendre polynomial PnP_n — the polynomials orthogonal on [−1,1][-1, 1]: ∫PmPn dx=0\int P_mP_n\,dx = 0 for m≠nm \ne n. With the weights from integrating the interpolant, the rule is exact for every polynomial of degree ≤2n−1\le 2n - 1, and no nn-point rule can do better. The weights are all positive and sum to 2, so the rule is stable for any nn — unlike high-order Newton–Cotes.
Derivation
The two-point rule by hand, and the first Legendre polynomials
Two points, four conditions. Look for ∫−11f dx≈w1f(x1)+w2f(x2)\int_{-1}^1 f\,dx \approx w_1f(x_1) + w_2f(x_2) exact for 1,x,x2,x31, x, x^2, x^3. Symmetry suggests x1=−ax_1 = -a, x2=ax_2 = a, w1=w2=ww_1 = w_2 = w; then the odd powers xx and x3x^3 give 0=00 = 0 automatically, and
1: 2w=2  ⇒  w=1,x2: 2a2=∫−11x2 dx=23  ⇒  a=13=0.577350.1:\ 2w = 2 \;\Rightarrow\; w = 1, \qquad x^2:\ 2a^2 = \int_{-1}^1x^2\,dx = \tfrac23 \;\Rightarrow\; a = \tfrac{1}{\sqrt3} = 0.577350.
Two evaluations, exact up to degree 3 — where the trapezoid rule's two evaluations stop at degree 1.
The Legendre polynomials. Build polynomials of rising degree, each orthogonal on [−1,1][-1, 1] to all lower ones (Gram–Schmidt on 1,x,x2,…1, x, x^2, \dots, lesson 5-2, with ∫−11pq dx\int_{-1}^1 pq\,dx as the inner product), scaled so Pn(1)=1P_n(1) = 1:
P0=1,P1=x,P2=12(3x2−1),P3=12(5x3−3x).P_0 = 1, \quad P_1 = x, \quad P_2 = \tfrac12(3x^2 - 1), \quad P_3 = \tfrac12(5x^3 - 3x).
Check one: ∫−11P2 dx=12(2−2)=0\int_{-1}^1 P_2\,dx = \tfrac12(2 - 2) = 0, and ∫−11xP2 dx=0\int_{-1}^1 xP_2\,dx = 0 because the integrand is odd. The roots of P2P_2 are ±1/3\pm1/\sqrt3 — exactly the nodes found by hand. The roots of P3P_3 are 00 and ±3/5=±0.774597\pm\sqrt{3/5} = \pm0.774597, the three-point nodes in the table.
Why the weights are positive. Let ℓi\ell_i be the Lagrange basis polynomial of node ii (degree n−1n - 1). Then ℓi2\ell_i^2 has degree 2n−2≤2n−12n - 2 \le 2n - 1, so the rule integrates it exactly, and ℓi2(xj)\ell_i^2(x_j) is 1 at j=ij = i and 0 elsewhere:
0<∫−11ℓi2 dx=∑jwj ℓi(xj)2=wi.0 < \int_{-1}^1\ell_i^2\,dx = \sum_j w_j\,\ell_i(x_j)^2 = w_i.
nnodesweightsexact up to degree
1021
2±0.577350 (= ±1/√3)1, 13
30, ±0.7745970.888889, 0.5555565
50, ±0.538469, ±0.9061800.568889, 0.478629, 0.2369279

The first Gauss–Legendre rules on [−1, 1]

Why the roots of PnP_n? Take any polynomial pp of degree ≤2n−1\le 2n - 1 and divide it by PnP_n: p=qPn+rp = qP_n + r, with qq and rr of degree ≤n−1\le n - 1. Then ∫qPn=0\int qP_n = 0 by orthogonality, and at the roots of PnP_n, p(xi)=r(xi)p(x_i) = r(x_i). So integrating pp reduces to integrating rr, which any nn-point interpolatory rule does exactly. Choosing the nodes at the roots makes the rule exact for the qPnqP_n part for free.
Any interval. Map [a,b][a, b] to [−1,1][-1, 1]: ∫abf dx=b−a2∑wif(b−a2xi+a+b2)\int_a^b f\,dx = \tfrac{b-a}{2}\sum w_i f\big(\tfrac{b-a}{2}x_i + \tfrac{a+b}{2}\big). Nodes and weights come from numpy.polynomial.legendre.leggauss(n), computed by the Golub–Welsch algorithm (1969) as eigenvalues of a small symmetric tridiagonal matrix — Unit 6 again.
📐Gauss versus Simpson
For ∫02excos⁡2x dx\int_0^2 e^x\cos 2x\,dx, with the same number of function evaluations:
  • 3 evaluations: Gauss error 1.6×10−21.6 \times 10^{-2}, Simpson 0.620.62;
  • 5: 5.4×10−65.4 \times 10^{-6} against 2.3×10−22.3 \times 10^{-2};
  • 9: 10−1410^{-14} against 1.2×10−31.2 \times 10^{-3}.
For a smooth integrand, Gauss's error falls geometrically, like Chebyshev interpolation (lesson 7-2). One exception: for a periodic integrand over a full period, the plain trapezoid rule is already spectrally accurate — 32 points give machine precision for ∫02π(1+12cos⁡θ)−1dθ\int_0^{2\pi}(1 + \tfrac12\cos\theta)^{-1}d\theta — which is why orbit averages are taken with equally spaced samples in the anomaly.
Mathematics is the queen of the sciences.
– Carl Friedrich Gauss, as quoted in Gauss zum Gedächtniss by W. Sartorius von Waltershausen (1856)
🛰Gauss points in every finite element
A finite-element stiffness matrix is an integral over each element: Ke=∫BTDB dVK_e = \int B^{\mathsf T}DB\,dV. Codes evaluate it with Gauss quadrature in the element's reference coordinates. For a two-node bar whose area varies linearly, the integrand is linear and one Gauss point gives the exact stiffness (30.0000 MN/m in the lab); for a quadratic taper one point is 5% off and two are exact. Choosing too few points ('reduced integration') makes elements too flexible and can admit spurious zero-energy 'hourglass' modes; too many wastes time and can lock thin elements. Choosing the number of Gauss points is a design decision every FE analyst makes, usually without knowing it.
💡Where you will use this
Structures, Geotechnics and Mechanical design: element matrices and stress recovery at Gauss points. Optimal Guidance: Gauss and Legendre pseudospectral methods for trajectory optimisation use these nodes as collocation points. Orbital Mechanics: Gauss–Legendre averaging of perturbations; Gauss's own method of orbit determination. Probability: Gauss–Hermite rules compute expectations over normal distributions with a handful of points — the basis of the unscented Kalman filter's sigma points (Unit 12).
📚Further reading
  • Gene Golub and John Welsch, 'Calculation of Gauss quadrature rules', Mathematics of Computation 23(106), 221–230 (1969) — nodes and weights as an eigenvalue problem.
  • Lloyd N. Trefethen, 'Is Gauss quadrature better than Clenshaw–Curtis?', SIAM Review 50(1), 67–87 (2008) — a readable comparison with experiments.
  • Philip Davis and Philip Rabinowitz, Methods of Numerical Integration (2nd ed., Academic Press, 1984) — Gauss rules of every kind, with tables.
↓But there is a problem
Every rule in this unit samples along a line. In two dimensions we take products — n2n^2 points; in three, n3n^3. A reliability engineer propagating twenty uncertain inputs through a structural model needs a twenty-dimensional integral: even three Gauss points per dimension means 320≈3.53^{20} \approx 3.5 billion model runs. In 1946 Stanisław Ulam, recovering from an illness in Los Angeles, passed the time playing solitaire, and wondered what fraction of games could be won. The combinatorics were hopeless. It would be far easier, he realised, to deal out a hundred games and simply count.