Gaussian Quadrature

Every engineer, physicist, and data scientist eventually needs to compute an integral they cannot solve by hand. The obvious approach is to sample the function at many evenly spaced points and add up the pieces — the trapezoidal rule, or its smoother cousin Simpson's rule. Double the number of samples and you roughly halve the error. Eventually you converge, but slowly.

Carl Friedrich Gauss asked a sharper question in 1814: if I am free to choose where I sample a function, not just how many times, where should I look? The answer he found was remarkable. For polynomials, you can integrate exactly — zero error — using far fewer evaluations than you might think, simply by placing the sample points at the roots of a special family of polynomials called the Legendre polynomials.

The technique he discovered, now called Gaussian quadrature, is one of the great efficiency wins in all of numerical analysis. An $n$-point Gaussian rule integrates every polynomial of degree up to $2n - 1$ exactly. Three points suffice for any degree-5 polynomial. Five points handle degree 9. The number of free parameters — $n$ points plus $n$ weights, totalling $2n$ — is used to the last degree of freedom.

This is not a story about approximation getting lucky. It is about an optimal allocation of effort, and the mathematics behind it connects to polynomial approximation and the deep structure of orthogonal polynomials.

Choose a polynomial and see 3-point Gaussian quadrature integrate it exactly — then compare it with the trapezoidal rule using the same three panels.

<!-- {{c_demo_title}} -->
<div class="controls">
  <label for="poly-select">{{label_poly}}</label>
  <select id="poly-select">
    <option value="cubic">{{opt_cubic}}</option>
    <option value="degree5" selected>{{opt_deg5}}</option>
    <option value="degree3b">{{opt_deg3b}}</option>
    <option value="degree5b">{{opt_deg5b}}</option>
  </select>
</div>
<canvas id="chart" width="480" height="220"></canvas>
<div id="results" class="results"></div>
<div class="btns">
  <button id="btn-gauss" type="button">{{btn_gauss}}</button>
  <button id="btn-trap" type="button" class="ghost">{{btn_trap}}</button>
  <button id="btn-both" type="button" class="ghost">{{btn_both}}</button>
</div>

/* {{c_layout}} */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; padding: .4rem .5rem; }
.controls { display: flex; align-items: center; gap: .5rem; margin-bottom: .5rem; flex-wrap: wrap; }
label { font-size: .85rem; font-weight: 600; color: #444; }
select { font-size: .85rem; padding: .25rem .4rem; border: 1px solid #cdd9e3; border-radius: 6px; }
canvas { display: block; width: 100%; max-width: 480px; border: 1px solid #e0e8ef; border-radius: 8px; background: #f9fbfd; }
.results { font-size: .85rem; margin: .5rem 0; line-height: 1.6; min-height: 3.5rem; }
.results table { border-collapse: collapse; width: 100%; }
.results td, .results th { padding: .18rem .4rem; border-bottom: 1px solid #e8eef3; }
.results th { font-weight: 700; color: #1d3557; text-align: left; }
.gauss-val { color: #0a7d33; font-weight: 700; }
.trap-val  { color: #c92f3c; font-weight: 600; }
.exact-val { color: #555; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; }
button { font: 600 13px system-ui, sans-serif; padding: .4rem .8rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button.active { background: #0a7d33; border-color: #0a7d33; color: #fff; }

// Code not found

Notice what happens: the 3-point Gauss rule (nodes at $x \approx -0.775,\, 0,\, 0.775$ with carefully tuned weights) returns the exact answer for any polynomial up to degree 5. The trapezoidal rule with the same number of function evaluations can be wildly off. The difference is not effort — it is where the effort is placed.

Why does the trick work? The answer lives in orthogonal polynomials.

An $n$-point quadrature rule has $2n$ free parameters: the $n$ nodes $x_1, \dots, x_n$ and the $n$ weights $w_1, \dots, w_n$. A polynomial of degree $d$ is pinned down by $d + 1$ coefficients, so $2n$ free parameters can, in principle, enforce exactness for degrees $0$ through $2n - 1$ — that is $2n$ conditions. The question is whether a consistent solution always exists, and whether it is unique.

Gauss's key insight: set the nodes to be the roots of the degree-$n$ Legendre polynomial $P_n(x)$. These polynomials form an orthogonal family on $[-1, 1]$:

$$\int_{-1}^{1} P_m(x)\, P_n(x)\, dx = 0 \quad \text{when } m \neq n.$$

This orthogonality guarantees that the weights $w_i$, which are then uniquely determined by the Vandermonde-like system, are all positive — so the rule is numerically stable. Moreover, any polynomial $f$ of degree at most $2n - 1$ can be divided by $P_n$:

$$f(x) = q(x)\, P_n(x) + r(x)$$

where $q$ and $r$ each have degree at most $n - 1$. Integrating: the first term vanishes by orthogonality (since $q$ has degree $< n$), and the second term $r$ is integrated exactly by the $n$-point rule. Zero error.

The three Gauss-Legendre nodes for $n = 3$ on $[-1, 1]$ are $0$ and $\pm\sqrt{3/5}$, with weights $8/9$ and $5/9$ respectively. Every polynomial of degree $\leq 5$ is integrated exactly — provably, not approximately.

For non-polynomial functions the rule still performs well: the error shrinks as $f$ becomes smoother, and for analytic functions the convergence with $n$ is exponential — each extra node buys you far more accuracy than any fixed-step rule can match.

Gaussian quadrature is not a theoretical curiosity. It is the silent engine behind much of modern scientific computing:

• Finite-element analysis: every element in a structural or fluid simulation requires dozens of accurate integrals per timestep. Gauss points are the standard choice in commercial FEA codes.
• Computer graphics: rendering integrals that compute how light bounces between surfaces use Gaussian quadrature (and its generalizations) to achieve smooth, accurate results.
• Probability and statistics: computing expectations under smooth distributions — the bread and butter of Bayesian inference and Monte Carlo methods — benefits enormously from adaptive Gaussian rules.
• Special functions: the values of Bessel functions, elliptic integrals, and error functions used throughout physics and engineering are tabulated or computed with variants of Gaussian quadrature.
• Machine learning kernels: many kernel functions in support-vector machines and Gaussian processes require numerical integration; the accuracy of Gauss rules keeps the approximations tight.

The idea also generalises: Gauss-Chebyshev, Gauss-Laguerre, and Gauss-Hermite rules target different weight functions and domains, each derived by matching nodes to the roots of the corresponding orthogonal-polynomial family. One principle, an entire zoo of tools.

Gaussian quadrature makes a statement that goes beyond numerical analysis: where you measure matters as much as how often you measure. By placing sample points at the roots of Legendre polynomials and deriving the matching weights, Gauss extracted the maximum possible information from each function evaluation — hitting the theoretical ceiling of $2n - 1$ exactness with exactly $n$ points.

The lesson generalises. In signal processing, in statistical design of experiments, in machine learning, the insight that optimal placement of queries beats brute-force density keeps reappearing. Gaussian quadrature is perhaps its clearest expression: a closed-form recipe that provably cannot be beaten, derived in 1814, and still running inside the finite-element solver on your laptop.