Numerical Integration

⚡ Key Takeaways

Numerical integration (or quadrature) approximates a definite integral $\int_a^b f(x)\,dx$ by a finite weighted sum of function values $\sum_i w_i f(x_i)$ . It is the standard tool when an antiderivative is unavailable, unknown, or only sampled data exist. Core methods range from the trapezoidal and Simpson's rules to Gaussian quadrature and adaptive/Monte Carlo schemes. Developed from Newton and Cotes (1700s) through Gauss (1814) to modern adaptive algorithms.

💡 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.

DefinitionQuadrature Rule

A quadrature rule for approximating $I[f] = \int_a^b f(x)\,dx$ is a formula

1

\quad Q_n[f] = \sum_{i=0}^{n} w_i\, f(x_i) \approx \int_a^b f(x)\,dx

Derivation hint: integrate an interpolating polynomial $p(x)=\sum_i f(x_i)\ell_i(x)$ and set $w_i=\int_a^b \ell_i(x)\,dx$ , using linearity of the integral. 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]$ .

DefinitionNewton–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:

2

\quad \int_a^b f(x)\,dx \approx \frac{h}{2}\Big[f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n)\Big]\quad\text{(composite trapezoidal)}

Derivation hint: on each subinterval $[x_i,x_{i+1}]$ approximate $f$ by the straight line through its endpoints (area of a trapezoid $=\tfrac{h}{2}[f(x_i)+f(x_{i+1})]$ ), then sum.

3

\quad \int_a^b f(x)\,dx \approx \frac{h}{3}\Big[f(x_0) + 4\sum_{i\,\text{odd}} f(x_i) + 2\sum_{i\,\text{even}} f(x_i) + f(x_n)\Big]\quad\text{(composite Simpson)}

Derivation hint: fit a parabola through each consecutive triple $(x_{i-1},x_i,x_{i+1})$ , integrate exactly to get $\tfrac{h}{3}[f_{i-1}+4f_i+f_{i+1}]$ , and sum over pairs of panels (requires $n$ even).

TheoremError 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)$ ,

4

\quad E_{\text{trap}} = -\frac{(b-a)h^2}{12}\,f''(\xi) = O(h^2)

Derivation hint: expand $f$ in a Taylor series about each subinterval midpoint and integrate the remainder, or apply the Peano kernel theorem; sum the local errors and use the intermediate value theorem on $f''$ . If $f\in C^4[a,b]$ , the composite Simpson rule satisfies, for some $\eta\in(a,b)$ ,

5

\quad E_{\text{Simp}} = -\frac{(b-a)h^4}{180}\,f^{(4)}(\eta) = O(h^4)

Derivation hint: Simpson's local error per double-panel is $-\tfrac{h^5}{90}f^{(4)}(\xi_i)$ ; sum over $n/2$ panels and apply the intermediate value theorem.

TheoremGaussian 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

6

\quad \int_{-1}^{1} f(x)\,dx = \sum_{i=0}^{n} w_i\, f(x_i),\qquad w_i = \frac{2}{(1-x_i^2)\,[P_{n+1}'(x_i)]^2}

Derivation hint: require the rule to be exact for $1,x,\dots,x^{2n+1}$ ; the constraints force the nodes to be roots of the orthogonal polynomial $P_{n+1}$ , since any such $f$ decomposes as $f=qP_{n+1}+r$ with $\int qP_{n+1}=0$ by orthogonality. 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 = \tfrac{a+b}{2} + \tfrac{b-a}{2}t$ , which multiplies the result by $\tfrac{b-a}{2}$ .

ProofInterpolatory Quadrature and Its Weights

Goal: show that any rule built by integrating a polynomial interpolant has the form $\sum_i w_i f(x_i)$ with $f$ -independent weights, and is exact for low-degree polynomials. 1. Choose nodes. Fix distinct nodes $x_0,\dots,x_n\in[a,b]$ . 2. Interpolate. The unique degree- $\le n$ Lagrange interpolant of $f$ at these nodes is

p(x) = \sum_{i=0}^n f(x_i)\,\ell_i(x),\qquad \ell_i(x)=\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}.

3. Integrate the interpolant. By linearity,

\int_a^b p(x)\,dx = \sum_{i=0}^n f(x_i)\underbrace{\int_a^b \ell_i(x)\,dx}_{=:\,w_i}.

The weights $w_i = \int_a^b \ell_i(x)\,dx$ depend only on the nodes, establishing equation [1]. 4. 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$ . 5. Error from interpolation error. Since $f(x)-p(x) = \dfrac{f^{(n+1)}(\zeta_x)}{(n+1)!}\prod_{i}(x-x_i)$ , integrating gives

E_n[f] = \int_a^b \frac{f^{(n+1)}(\zeta_x)}{(n+1)!}\prod_i (x-x_i)\,dx,

which is small when $f^{(n+1)}$ is bounded and the node polynomial $\prod_i(x-x_i)$ has small integrated magnitude. $\blacksquare$

∎

ProofWhy 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]$ . 2. Take any polynomial $f$ of degree $\le 2n+1$ . Divide by $P_{n+1}$ : $f = q\,P_{n+1} + r$ with $\deg q\le n$ , $\deg r\le n$ . 3. 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 ( $\deg q\le n$ ). 4. 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)$ . 5. 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$ . 6. Combining: $\sum_i w_i f(x_i) = \int_{-1}^1 r = \int_{-1}^1 f$ . The rule is exact to degree $2n+1$ . $\blacksquare$

∎

mathematician · 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.

🧠 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''$ 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$ .

Figure

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

ExampleTrapezoidal and Simpson for

\int_0^1 e^{-x^2}\,dx

Given: $f(x)=e^{-x^2}$ on $[0,1]$ with $n=4$ subintervals, $h=0.25$ . 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$

Trapezoidal (eq. [2]):

\frac{0.25}{2}\big[1.000000 + 2(0.939413+0.778801+0.569783) + 0.367879\big]

= 0.125\,[1.367879 + 2(2.287997)] = 0.125 \times 5.943873 = 0.742984.

Simpson (eq. [3]), odd $i=1,3$ , even $i=2$ :

\frac{0.25}{3}\big[1.000000 + 4(0.939413+0.569783) + 2(0.778801) + 0.367879\big]

= 0.0833333\,[1.367879 + 4(1.509196) + 1.557602] = 0.0833333 \times 8.962265.

\boxed{I_{\text{trap}} \approx 0.74298,\quad I_{\text{Simp}} \approx 0.746855\ (\text{dimensionless})}

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

ExampleTwo-point Gauss–Legendre for

\int_1^3 \ln x\,dx

Given: $f(x)=\ln x$ , interval $[1,3]$ . Use the 2-point Gauss rule: on $[-1,1]$ nodes $t=\pm 1/\sqrt3$ , weights $w=1$ each.

Step 1 — map $[-1,1]\to[1,3]$ : midpoint $\tfrac{a+b}{2}=2$ , half-width $\tfrac{b-a}{2}=1$ , so $x=2+t$ .

$t_1=-1/\sqrt3 \approx -0.577350 \Rightarrow x_1 = 1.422650$
$t_2=+1/\sqrt3 \approx +0.577350 \Rightarrow x_2 = 2.577350$

Step 2 — evaluate:

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

Step 3 — apply rule (multiply by half-width 1):

I \approx 1\cdot\big[1\cdot 0.352573 + 1\cdot 0.946785\big] = 1.299358.

Exact: $\int_1^3\ln x\,dx = [x\ln x - x]_1^3 = (3\ln3 - 3)-(0-1) = 3\ln3 - 2 = 1.295837.$

\boxed{I_{\text{Gauss-2}} \approx 1.29936,\quad \text{error}\approx 3.5\times10^{-3}}

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

⚠️ Common mistake

More nodes is not always better with equal spacing. High-degree Newton–Cotes formulas on equally spaced points develop negative weights and suffer Runge's phenomenon — the interpolating polynomial oscillates wildly near the endpoints. Negative weights destroy numerical stability (cancellation amplifies rounding error) and can make the rule diverge as $n\to\infty$ for some smooth functions. The cure is either compositing (many low-order panels) or non-uniform nodes (Gauss, Clenshaw–Curtis).

CounterexampleRunge's function defeats high-order Newton–Cotes

Consider $f(x) = \dfrac{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=\pm1$ . 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.

⚠️ Common mistake

Simpson's rule requires an even number of subintervals. Each parabola spans two panels. With an odd $n$ , you cannot pair the panels; applying the weights $1,4,2,4,\dots$ blindly gives a wrong answer. Handle a leftover panel separately (e.g. Simpson's 3/8 rule on three panels) or choose $n$ even.

connection

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.

🔬 Where this appears

⏳ History

- 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.

open-problem

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.

Formal Statement

Why It Holds / Derivation

Building Intuition

Worked Examples

Common Misconceptions

Connections

Real-World Applications

Historical Development

Subtleties & Limits

See Also