Numerical Integration and Quadrature

Numerically approximating definite integrals, a process known as quadrature, is a cornerstone of scientific computing. When an antiderivative is unavailable or a function is only known at discrete data points, these algorithms transform an intractable analytic problem into a solvable numerical one. Mastering the trade-offs between accuracy, computational cost, and stability is essential for effectively applying these tools to simulations, data analysis, and high-dimensional models.

The Foundation: Newton-Cotes Formulas

Newton-Cotes formulas provide the most intuitive entry point to numerical integration. They approximate a function $f(x)$ over an interval $[a,b]$ by replacing it with a simple polynomial interpolant at equally spaced nodes, then integrating that polynomial exactly. The resulting integral is a weighted sum of function values: ${\int }_a^{b}f(x)dx\approx {\sum }_{i=0}^n{w}_{i}f(x_i)$.

The simplest rules are the low-degree cases:

• The Trapezoidal Rule uses a linear interpolant ($n=1$). It approximates the area under the curve as a trapezoid: ${\int }_a^{b}f(x)dx\approx (b-a)\left[\frac{f(a)+f(b)}{2}\right]$.
• Simpson's Rule uses a quadratic interpolant ($n=2$), fitting a parabola through three points: ${\int }_a^{b}f(x)dx\approx \frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]$.

These basic rules are often applied in composite forms, where the interval $[a,b]$ is subdivided into $m$ smaller panels and the rule is applied to each. The error analysis for these methods reveals a crucial pattern. The error term for the composite Trapezoidal Rule is proportional to $(b-a{)}^3{f}^{\prime \prime }(\xi )/m^2$, meaning it has a convergence rate of $O(m^{-2})$. Simpson's Rule converges faster at $O(m^{-4})$, provided the function is sufficiently smooth. This illustrates a general principle: higher-degree Newton-Cotes formulas can offer higher order convergence, but they can also suffer from numerical instability (Runge's phenomenon) for very high degrees and are inefficient if the function is expensive to evaluate.

Optimal Sampling: Gaussian Quadrature

While Newton-Cotes formulas fix the node locations (as equally spaced points), Gaussian quadrature optimizes both the nodes $x_i$ and the weights $w_i$ to achieve the highest possible degree of precision. For an $n$-point rule, Gaussian quadrature is designed to be exact for all polynomials of degree up to $2n-1$, double what a Newton-Cotes rule using the same number of nodes can achieve. This dramatic gain in efficiency—achieving high accuracy with fewer function evaluations—makes it a preferred method for integrating smooth functions.

The nodes for Gaussian quadrature on the interval $[-1,1]$ are the roots of Legendre polynomials. For a general interval $[a,b]$, a linear change of variable is applied. The weights are then derived from the theory of orthogonal polynomials. For example, a 2-point Gaussian rule on $[-1,1]$ is: $${\int }_{-1}^{1}f(x)dx\approx f\left(-\frac{1}{\sqrt{3}}\right)+f\left(\frac{1}{\sqrt{3}}\right)$$ and is exact for cubics. The error bound for Gaussian quadrature is complex but relates to the best polynomial approximation of $f(x)$. The key takeaway is its superior convergence rate for smooth integrands compared to Newton-Cotes methods using the same computational budget.

Intelligent Refinement: Adaptive Integration

Both composite Newton-Cotes and Gaussian quadrature use a fixed, uniform mesh. Adaptive integration algorithms dynamically allocate computational effort by concentrating subdivisions where the integrand is most troublesome (e.g., rapid oscillations, sharp peaks). A common strategy is adaptive Simpson's quadrature.

The algorithm works recursively:

1. Estimate the integral $I$ over an interval using Simpson's rule.
2. Split the interval in half and estimate the integral on each subinterval, summing them to get $I_L+I_R$.
3. Compare the estimates. If $|I-(I_L+I_R)|$ is less than a specified error tolerance, accept the refined estimate.
4. If the error is too large, recursively apply the same procedure to each subinterval.

This approach ensures that the final result meets a user-specified global accuracy requirement without wasting evaluations on regions where the function is well-behaved. It directly manages the trade-off between error bounds and computational cost, making it the workhorse for many one-dimensional integration problems in software libraries.

Tackling Complexity: Multidimensional Integration

Integrating over domains in two or more dimensions, termed cubature, introduces exponential complexity. A "curse of dimensionality" emerges: applying a 1D rule with $n$ points to a $d$-dimensional hypercube requires $n^d$ function evaluations, which quickly becomes prohibitive.

• Product Rules: The direct approach is to form a product of one-dimensional quadrature rules. For a double integral, this looks like:

$${\int }_a^b{\int }_c^{d}f(x,y)dydx\approx \sum _{i=1}^{n_x}\sum _{j=1}^{n_y}{w}_i^x{w}_j^{y}f(x_i,y_j).$$ While straightforward and effective for low dimensions (2D, 3D), its cost grows exponentially.

• Monte Carlo Integration: For high-dimensional integrals, Monte Carlo integration becomes essential. It estimates the integral by sampling the integrand at random points. The estimate is $V\cdot ⟨f⟩$, where $V$ is the volume of the integration region and $⟨f⟩$ is the average of $f$ over $N$ random samples. Crucially, the statistical error converges as $O(1/\sqrt{N})$, independent of dimension. This breaks the curse of dimensionality, making it the only feasible method for integrals in dozens or hundreds of dimensions, despite its relatively slow convergence rate. Variance reduction techniques (e.g., importance sampling, stratified sampling) are used to improve its efficiency.

Common Pitfalls

1. Applying High-Order Rules to Non-Smooth Functions: Using Gaussian quadrature or a high-order Newton-Cotes formula on a function with a discontinuity or cusp will yield poor results. The high-order convergence rates depend on the function's smoothness. For such functions, a low-order composite rule or an adaptive method that can localize the discontinuity is better.
2. Ignoring Error Estimates: Treating a single quadrature result as exact is dangerous. Always use an error estimate, whether it's the theoretical asymptotic bound, the difference between two rule orders (e.g., trapezoid vs. Simpson), or the local error estimate in an adaptive routine. Without it, you have no gauge of the answer's reliability.
3. Misapplying 1D Thinking to Multiple Dimensions: Attempting to use a product rule for a high-dimensional integral will lead to unacceptably long computation times. Recognize the dimensionality of your problem and switch to Monte Carlo or sparse grid methods when needed.
4. Overlooking the Integration Domain: The standard rules apply to simple intervals or hypercubes. Integrating over a complex or implicitly defined region requires a change of variables to map it to a standard shape or the use of specialized methods like Markov Chain Monte Carlo (MCMC).

Summary

• Newton-Cotes formulas (Trapezoidal, Simpson's) are built on polynomial interpolation at equally spaced nodes and are the basis for many composite and adaptive schemes.
• Gaussian quadrature optimizes node placement and weighting to achieve exact integration for polynomials up to degree $2n-1$, offering superior efficiency for smooth integrands.
• Adaptive integration algorithms, like adaptive Simpson, dynamically refine the mesh to control error efficiently, making them robust for functions with difficult behavior.
• For multidimensional integration, product rules work for low dimensions, but Monte Carlo methods are necessary for high-dimensional problems due to their dimension-independent convergence rate, overcoming the curse of dimensionality.
• The choice of method always involves balancing the convergence rate, the computational efficiency (number of function evaluations), and the need for reliable error bounds or estimates.