Numerical Integration

Numerical integration is the branch of numerical analysis concerned with computing approximate values of definite integrals using only evaluations of the integrand. It matters because the overwhelming majority of integrals arising in physics, engineering, statistics, and applied science have no closed-form antiderivative expressible in elementary functions — even something as simple as ${\int }_0^1{e}^{-x^2}dx$. The central strategy is to replace the integrand by an interpolating polynomial (or another tractable model) over each subinterval and integrate that exactly. A surprising fact is that by choosing the sample points cleverly rather than uniformly — as Gauss did — an $n$-point rule can integrate polynomials of degree up to $2n-1$ exactly, doubling the accuracy obtainable from equally spaced nodes.

The Core Insight

You cannot integrate an arbitrary function symbolically, but you can always evaluate it at chosen points. Replace $f$ on each piece of $[a,b]$ by a polynomial that agrees with $f$ at the sample points and integrate the polynomial exactly. The integral of the polynomial is a linear combination ${\sum }_i{w}_{i}f(x_i)$ where the weights $w_i$ depend only on the nodes — not on $f$. All of quadrature is the art of choosing nodes $x_i$ and weights $w_i$ to make the error small.

Formal Statement

Quadrature Rule

A quadrature rule for approximating $I[f]={\int }_a^{b}f(x)dx$ is a formula

where the points $x_i\in [a,b]$ are the nodes (or abscissae), the numbers $w_i$ are the weights, and $n+1$ is the number of points. The degree of exactness (or order) of the rule is the largest integer $m$ such that $Q_n[p]=I[p]$ for every polynomial $p$ of degree $\le m$. The quadrature error is $E_n[f]=I[f]-Q_n[f]$.

Newton–Cotes Formulas

The Newton–Cotes formulas use $n+1$ equally spaced nodes $x_i=a+ih$, where $h=(b-a)/n$, and take the weights from the integral of the Lagrange interpolating polynomial. The two most important closed forms are:

Error of the Composite Trapezoidal and Simpson Rules

If $f\in C^2[a,b]$, the composite trapezoidal rule satisfies, for some $\xi \in (a,b)$,

If $f\in C^4[a,b]$, the composite Simpson rule satisfies, for some $\eta \in (a,b)$,

Gaussian Quadrature

For the weight function $w(x)\equiv 1$ on $[-1,1]$, there exist $n+1$ nodes (the roots of the degree-$(n+1)$ Legendre polynomial $P_{n+1}$) and weights $w_i$ such that

This Gauss–Legendre rule with $n+1$ nodes is exact for all polynomials of degree $\le 2n+1$. A general interval $[a,b]$ is reached by the affine change of variable $x=\frac{a+b}{2}+\frac{b-a}{2}t$, which multiplies the result by $\frac{b-a}{2}$.

Why It Holds / Derivation

Interpolatory Quadrature and Its Weights

1. Choose nodes. Fix distinct nodes $x_0,\dots ,x_n\in [a,b]$.

1. Interpolate. The unique degree-$\le n$ Lagrange interpolant of $f$ at these nodes is

1. Integrate the interpolant. By linearity,

1. Degree of exactness. If $f$ is itself a polynomial of degree $\le n$, then the interpolant reproduces it exactly: $p=f$. Hence $Q_n[f]={\int }_a^{b}p={\int }_a^{b}f=I[f]$. So every $(n+1)$-point interpolatory rule is exact to degree at least $n$.

1. Error from interpolation error. Since $f(x)-p(x)=\frac{f^{(n+1)}({\zeta }_x)}{(n+1)!}{\prod }_i(x-x_i)$, integrating gives

Why Gaussian Nodes Double the Accuracy

1. Let the nodes $x_0,\dots ,x_n$ be the roots of $P_{n+1}$, the degree-$(n+1)$ Legendre polynomial, orthogonal to all lower-degree polynomials on $[-1,1]$.

1. Take any polynomial $f$ of degree $\le 2n+1$. Divide by $P_{n+1}$: $f=q{P}_{n+1}+r$ with $\mathrm{deg}q\le n$, $\mathrm{deg}r\le n$.

1. Exact integral: ${\int }_{-1}^{1}f={\int }_{-1}^{1}q{P}_{n+1}+{\int }_{-1}^{1}r=0+{\int }_{-1}^{1}r$, the first term vanishing by orthogonality ($\mathrm{deg}q\le n$).

1. Quadrature: at each node $P_{n+1}(x_i)=0$, so $f(x_i)=r(x_i)$. Thus ${\sum }_i{w}_{i}f(x_i)={\sum }_i{w}_{i}r(x_i)$.

1. The rule is interpolatory of degree $\ge n$, so it is exact for $r$ (degree $\le n$): ${\sum }_i{w}_{i}r(x_i)={\int }_{-1}^{1}r$.

1. Combining: ${\sum }_i{w}_{i}f(x_i)={\int }_{-1}^{1}r={\int }_{-1}^{1}f$. The rule is exact to degree $2n+1$. $■$

Carl Friedrich Gauss | 1777–1855

Gauss presented his quadrature method in Methodus nova integralium valores per approximationem inveniendi (1814). His insight — that the placement of nodes is a free resource to be optimized, not fixed by convenience — was revolutionary. By demanding exactness for as many monomials as the $2(n+1)$ free parameters (nodes + weights) allow, he discovered that the optimal nodes are precisely the roots of orthogonal polynomials, linking quadrature to the theory of orthogonal polynomials he and Legendre were developing. The result remains the gold standard for smooth integrands and underlies modern spectral methods.

Building Intuition

Think of estimating the area under a curve from a few altitude readings. The trapezoidal rule connects adjacent readings with straight lines — it overestimates concave-down curves and underestimates concave-up ones, hence the $-f^{\prime \prime }$ in its error. Simpson's rule fits little parabolas, capturing curvature, so it nails any cubic exactly (a happy accident from symmetry) and only fails at the fourth derivative.

The leap to Gaussian quadrature is the realization that equally spaced samples waste information. If you are allowed to pick where to measure, you place nodes denser near the ends (where polynomial interpolation is shakiest) and let each node "carry" twice its weight in exactness. Two numbers per node — its location and its weight — give two equations satisfied, so $n+1$ nodes satisfy $2n+2$ conditions, i.e. exactness through degree $2n+1$.

The governing principle throughout: convergence rate is dictated by smoothness. Doubling the number of panels cuts trapezoidal error by 4 ($h^2$) but Simpson's by 16 ($h^4$); for analytic functions, Gaussian and spectral methods converge faster than any power of $h$.

Trapezoidal vs. Simpson approximation of a curve over one panel pair

Worked Examples

Trapezoidal and Simpson for $\int_0^1 e^{-x^2}\,dx$

Sample values:

• $f(0)=1.000000$
• $f(0.25)=0.939413$
• $f(0.50)=0.778801$
• $f(0.75)=0.569783$
• $f(1.00)=0.367879$

The true value is $0.7468241\dots$; Simpson's error is $\sim 3\times 10^{-5}$ versus the trapezoidal error $\sim 4\times 10^{-3}$ — a hundredfold improvement at the same five evaluations.

Two-point Gauss–Legendre for $\int_1^3 \ln x\,dx$

• $t_1=-1/\sqrt{3}\approx -0.577350\Rightarrow x_1=1.422650$
• $t_2=+1/\sqrt{3}\approx +0.577350\Rightarrow x_2=2.577350$

• $\ln (1.422650)=0.352573$
• $\ln (2.577350)=0.946785$

Only two evaluations achieve accuracy comparable to a 4-panel Simpson rule — the power of optimal node placement.

Common Misconceptions

Runge's function defeats high-order Newton–Cotes

Consider $f(x)=\frac{1}{1+25x^2}$ on $[-1,1]$ — perfectly smooth, even analytic. The single-panel Newton–Cotes rule of degree $n$ (one polynomial through $n+1$ equally spaced nodes) has error that grows without bound as $n\to \infty$, because the interpolating polynomial $p_n$ oscillates with amplitude diverging near $x=\pm 1$. The integrated error does not go to zero. By contrast, Gauss–Legendre and the composite trapezoidal rule both converge for this function. The lesson: smoothness of $f$ does not guarantee convergence of a fixed-node scheme — node distribution matters as much as integrand regularity.

Connections

Numerical integration sits at the crossroads of several fields. The weights and exactness conditions come from Polynomial Interpolation and Lagrange basis polynomials. Gaussian quadrature is an application of Orthogonal Polynomials (Legendre, Chebyshev, Hermite, Laguerre — each suited to a weight function and interval). The error analysis relies on Taylor's Theorem and the Peano kernel theorem. Richardson Extrapolation applied to the trapezoidal rule produces Romberg integration, systematically cancelling the leading $h^2,h^4,\dots$ error terms. For high-dimensional integrals, where grid methods suffer the curse of dimensionality, one turns to Monte Carlo Integration whose $O(N^{-1/2})$ error is independent of dimension. Finally, quadrature is the engine inside Numerical Solution of ODEs, since integrating a differential equation step is exactly an integration problem.

Real-World Applications

• Computational physics: evaluating partition functions, scattering cross-sections, and matrix elements $⟨{\psi }_i|\hat{H}|{\psi }_j⟩$ in quantum chemistry uses Gauss–Hermite and Gauss–Legendre quadrature on every integral.
• Finite element analysis (FEA): structural and thermal simulations integrate stiffness matrices element-by-element with Gauss points; engineers literally speak of "2×2 Gauss integration" inside each element.
• Computer graphics and rendering: the rendering equation's light-transport integrals are evaluated by Monte Carlo and quasi-Monte Carlo quadrature (path tracing).
• Statistics and Bayesian inference: marginal likelihoods and posterior expectations $\int f(\theta )p(\theta )d\theta$ use Gauss–Hermite quadrature (for Gaussian-like posteriors) or adaptive Monte Carlo.
• Signal processing: computing energies, correlations, and Fourier coefficients of sampled signals via trapezoidal/Simpson integration of discrete data.
• Astrodynamics and orbital mechanics: propagating spacecraft trajectories and computing perturbation integrals with high-order Gauss and Romberg schemes.
• Financial engineering: pricing options requires integrating payoff functions against probability densities (e.g. Heston model characteristic-function integrals via adaptive Gauss–Kronrod).

Historical Development

• c. 250 BCE — Archimedes: computes areas by the method of exhaustion, the conceptual ancestor of quadrature; his parabola result anticipates Simpson's rule.
• 1670s — Newton & Gregory: interpolation formulas on equally spaced points lay groundwork for Newton–Cotes rules.
• 1722 — Roger Cotes: publishes the family of equal-interval quadrature weights now called Newton–Cotes formulas (posthumously, after his 1716 death).
• 1743 — Thomas Simpson: popularizes the parabolic rule (known earlier to Cavalieri and Gregory) in his textbook, lending it his name.
• 1814 — Carl Friedrich Gauss: introduces optimal-node quadrature, achieving degree $2n+1$ exactness; ties nodes to roots of Legendre polynomials.
• 1826 — Carl Jacobi: generalizes to weighted integrals and the connection with orthogonal polynomials, broadening Gauss's method.
• 1900 — Walter Romberg's precursors (Richardson): L. F. Richardson develops extrapolation (1911) that later yields Romberg integration.
• 1955 — Werner Romberg: publishes the systematic $h^2$-extrapolation tableau combining trapezoidal estimates into high-order results.
• 1960s — Clenshaw, Curtis, Kronrod: Clenshaw–Curtis quadrature (Chebyshev nodes, 1960) and Gauss–Kronrod extensions (1964) bring efficient, error-estimating adaptive quadrature into wide computational use.

Subtleties & Limits

Several issues remain practically and theoretically delicate. (1) Singular and oscillatory integrands: standard rules degrade catastrophically for ${\int }_0^1{x}^{-1/2}f(x)dx$ or highly oscillatory $\int f(x)e^{i\omega x}dx$ with large $\omega$; specialized methods (Levin, Filon, modified Gauss weights, double-exponential "tanh-sinh" transforms) are needed, and choosing among them automatically is unsolved in full generality. (2) High-dimensional quadrature: the curse of dimensionality makes tensor-product grids ($n^d$ points) infeasible; sparse grids (Smolyak) and quasi-Monte Carlo lattice rules push the boundary, but optimal constructions for general integrands remain an active research area. (3) A posteriori error estimation: rigorous, tight, automatically computable error bounds without knowing $f^{(k)}$ are hard; practical codes use embedded pairs (Gauss–Kronrod) heuristically, which can be fooled. (4) Verified/interval quadrature producing guaranteed enclosures is an ongoing field at the intersection of numerical analysis and computer-assisted proof.

See Also