Gauss Quadrature

The idea — both nodes and weights are free parameters

The standard Newton–Cotes quadrature rules (trapezoidal, Simpson) pin the nodes at equally spaced points and choose the weights to maximise the degree of precision — the largest degree of polynomial the rule integrates exactly. With n equally spaced nodes, the best Newton–Cotes rule integrates polynomials up to degree n − 1 or n exactly.

Gauss's insight (1814): let the nodes be free parameters too. With n nodes and n weights, you have 2n unknowns. Insist that the rule integrates every monomial 1, x, x², …, x^{2n−1} exactly. That's 2n linear constraints — 2n equations, 2n unknowns. Generically solvable. The resulting rule is exact for any polynomial of degree ≤ 2n − 1 — twice the polynomial degree of the equally spaced rule.

∫_a^b f(x) w(x) dx ≈ ∑_{i=1}^n w_i · f(x_i)
nodes x_i and weights w_i chosen so that the rule is exact for x⁰, x¹, ..., x^(2n − 1)

The miracle is what the nodes turn out to be: the n roots of the degree-n orthogonal polynomial p_n(x) associated with the weight w on [a, b]. Gauss's discovery was that the optimal locations are exactly the roots of orthogonal polynomials.

Why the nodes are orthogonal-polynomial roots

Take any polynomial f of degree ≤ 2n − 1. Divide it by p_n(x):

f(x) = q(x) · p_n(x) + r(x),    deg q ≤ n − 1,  deg r ≤ n − 1

Now compute the integral:

∫ f w dx = ∫ q · p_n · w dx + ∫ r · w dx
        = 0 + ∫ r · w dx          [because q has degree < n and p_n is orthogonal to all such q]

And evaluate the quadrature sum at nodes x_i = roots of p_n:

∑ w_i · f(x_i) = ∑ w_i · (q(x_i) p_n(x_i) + r(x_i))
              = ∑ w_i · (q(x_i) · 0 + r(x_i))         [p_n(x_i) = 0 by construction]
              = ∑ w_i · r(x_i)

So matching ∫ f w dx to Σ w_i f(x_i) reduces to matching ∫ r w dx to Σ w_i r(x_i) — but r has degree ≤ n − 1, which only needs n parameters to integrate exactly. So choosing the n weights to be exact for {1, x, …, x^{n−1}} suffices. The full quadrature rule is then exact for polynomials up to degree 2n − 1, as advertised.

Worked example — Gauss-Legendre with n = 3

For ∫_{−1}^1 f(x) dx (weight w = 1), the relevant orthogonal polynomial is P_3(x) = (5x³ − 3x)/2. Its three roots are:

x_1 = −√(3/5) ≈ −0.7746
x_2 = 0
x_3 = +√(3/5) ≈ +0.7746

The corresponding weights (computed from Lagrange interpolation at the nodes) are:

w_1 = w_3 = 5/9 ≈ 0.5556
w_2 = 8/9 ≈ 0.8889

So the 3-point Gauss-Legendre rule is:

∫_{−1}^1 f(x) dx ≈ (5/9) f(−√(3/5)) + (8/9) f(0) + (5/9) f(√(3/5))

This rule is exact for any polynomial of degree ≤ 5. Compare Simpson's rule with the same 3 evaluation points (at x = −1, 0, +1), which is exact only up to degree 3. With 3 function evaluations, Gauss reaches degree of precision 5 vs Simpson's 3 — a doubling.

Try it on ∫_{−1}^1 (x⁴ − x² + 1) dx, which has exact value 2 · (1/5 − 1/3 + 1) = 2 · (13/15) = 26/15 ≈ 1.7333:

f(±√(0.6)) = 0.6² − 0.6 + 1 = 0.36 − 0.6 + 1 = 0.76
f(0) = 1
Gauss-3 ≈ (5/9)(0.76) + (8/9)(1) + (5/9)(0.76) = (10/9)(0.76) + 8/9
       = 7.6/9 + 8/9 = 15.6/9 ≈ 1.7333 ✓

Exact, as the theory promises (the integrand is degree 4 ≤ 5).

Common Gauss-Legendre rules (interval [−1, 1], weight 1)

n	Nodes x_i	Weights w_i	Degree of precision
1	0	2	1 (midpoint rule)
2	±1/√3 ≈ ±0.57735	1, 1	3
3	0, ±√(3/5) ≈ ±0.77460	8/9, 5/9, 5/9	5
4	±0.33998, ±0.86114	0.65215, 0.34785	7
5	0, ±0.53847, ±0.90618	0.56889, 0.47863, 0.23693	9
10	see tables	see tables	19

Nodes and weights are symmetric about 0 for symmetric weight functions. They are tabulated in standard references (Abramowitz–Stegun, NIST DLMF) to high precision.

Gauss-Hermite and Gauss-Laguerre — integrating over unbounded intervals

For integrals of the form ∫_{−∞}^∞ g(x) e^{−x²} dx — common in statistics and quantum mechanics — Gauss-Hermite is the rule:

∫_{−∞}^∞ g(x) e^{−x²} dx ≈ ∑ w_i g(x_i),  x_i = roots of H_n(x)

For n = 5: x_i ≈ 0, ±0.95857, ±2.02018; w_i ≈ 0.94531, 0.39362, 0.01996. Plugging in a polynomial g(x) of degree ≤ 9 gives the exact integral. To compute ∫ g(x) e^{−x²/2} dx (probabilists' Gaussian) just rescale x.

For ∫_0^∞ g(x) e^{−x} dx, Gauss-Laguerre uses roots of L_n(x). Useful for half-line integrals in radial quantum problems, lifetimes, and queueing theory.

JavaScript — Gauss-Legendre quadrature

// Hand-tabulated Gauss-Legendre nodes and weights on [-1, 1]
const GL_TABLES = {
  2: { x: [-0.5773502691896257, 0.5773502691896257],
       w: [1.0, 1.0] },
  3: { x: [-0.7745966692414834, 0, 0.7745966692414834],
       w: [5/9, 8/9, 5/9] },
  5: { x: [-0.9061798459386640, -0.5384693101056831, 0,
            0.5384693101056831, 0.9061798459386640],
       w: [0.2369268850561891, 0.4786286704993665,
           0.5688888888888889, 0.4786286704993665,
           0.2369268850561891] },
};

function gaussLegendre(f, a, b, n = 5) {
  const tab = GL_TABLES[n];
  if (!tab) throw new Error("Unsupported n; use 2, 3, or 5");
  // Map [-1, 1] -> [a, b]
  const half = (b - a) / 2;
  const mid  = (a + b) / 2;
  let sum = 0;
  for (let i = 0; i < n; i++) {
    const x = half * tab.x[i] + mid;
    sum += tab.w[i] * f(x);
  }
  return half * sum;
}

// Test: ∫_0^1 e^{-x²} dx ≈ 0.7468241...
const gauss = x => Math.exp(-x * x);
console.log(gaussLegendre(gauss, 0, 1, 5));  // ≈ 0.7468242
// 5-point Simpson would need ~15 evaluations for similar accuracy

// Test: ∫_{-1}^1 (x^4 - x^2 + 1) dx = 26/15 ≈ 1.7333
console.log(gaussLegendre(x => x**4 - x**2 + 1, -1, 1, 3));  // 1.7333... (exact)

Gauss vs Newton–Cotes — head to head

	Trapezoid (n nodes)	Simpson (n = 3)	Gauss-Legendre (n nodes)
Nodes	Equally spaced	Equally spaced	Orthogonal-polynomial roots
Degree of precision	n − 1	3	2n − 1
Per-evaluation accuracy	Low	Medium	Best possible
Nestable for adaptive use	Yes	Yes (Simpson-1/3)	Only Gauss-Kronrod variants
Best for	Data on a grid	Smooth functions, equally spaced	Smooth functions, freely chosen nodes
Practical use	Tabular data	Common in engineering	FEM, spectral methods, quantum chem

The reason adaptive quadrature codes (SciPy quad, GSL, Quadpack) use Gauss-Kronrod is that Kronrod adds 2n + 1 strategically chosen nodes to a Gauss rule to allow re-estimating the error without re-running the full computation. Pure Gauss is optimal for smooth integrands with fixed n; Kronrod handles adaptation.

Where Gauss quadrature lives

• Finite element methods. Element stiffness matrices ∫ B^T D B dx are integrated using Gauss-Legendre nodes on each element. 2-point rule for quadratic elements, 3-point for cubic; the FEM convergence theory hinges on Gauss exactness up to 2n − 1.
• Spectral methods. Pseudospectral collocation places the equation residual to zero at Gauss-Lobatto (or Gauss-Legendre) nodes. Exponentially convergent for smooth solutions.
• Computational physics and chemistry. Electron correlation, perturbation theory, scattering — every multi-electron integral uses Gauss quadrature in some basis.
• Bayesian inference. Marginalising likelihoods over Gaussian priors via Gauss-Hermite; alternatives to Monte Carlo when dimension is low.
• Computer graphics — Monte Carlo light transport. Quasi-Gauss schemes evaluate hemisphere integrals for BRDF and global illumination.
• Pricing and economic models. Numerical integration over option payoffs and probability distributions in quantitative finance.

Common mistakes

• Picking the wrong weight family for the integral. ∫ f(x) e^{−x²} dx is a Gauss-Hermite integral, not Gauss-Legendre. If you use the wrong weight, you compute ∫ f(x) e^{−x²} dx incorrectly because the weight is baked into the rule.
• Applying high-n Gauss to non-smooth integrands. High-degree exactness is wasted on functions with jumps or singularities. Subdivide the domain at the singularity or switch to special-purpose rules.
• Forgetting to map [a, b] → [−1, 1]. Gauss-Legendre tables are on [−1, 1]. For ∫_a^b f, substitute x = ((b − a)t + a + b)/2 and multiply by (b − a)/2 in front of the sum.
• Trying to add nodes to refine. Gauss nodes are not nested — refining n = 3 to n = 5 throws away all your previous evaluations. Use Gauss-Kronrod or Clenshaw-Curtis if nested rules are needed.
• Assuming Gauss is always best. For periodic integrands on a periodic interval, the trapezoidal rule converges spectrally (exponentially) — better than Gauss-Legendre. Gauss wins for non-periodic smooth integrands.
• Computing very high n nodes/weights from scratch each time. For n < 100, precomputed tables are faster than the Golub-Welsch algorithm. NumPy and SciPy cache nodes/weights for repeated calls.