Fundamental Numerical Integration Methods

Why This Matters

Numerical integration is what you turn to whenever a closed-form antiderivative doesn't exist, which is most of the time in real applications. In Numerical Analysis II, you're expected to go beyond applying formulas: you need to understand error behavior, convergence rates, and when to choose one method over another. These concepts tie directly into polynomial interpolation, Taylor series error bounds, and approximation theory more broadly.

What separates strong exam performance from mediocre recall is understanding the why behind each method. Why does Simpson's Rule outperform the Trapezoidal Rule? Why might you abandon classical methods entirely for Monte Carlo in high dimensions? For each method, know what order of accuracy it achieves, what assumptions it requires, and what trade-offs it makes between function evaluations and precision.

---

Polynomial-Based Interpolation Methods

These foundational methods approximate the integrand using polynomials of increasing degree over subintervals. Higher-degree polynomial interpolation generally yields faster error convergence, but there are important caveats once the degree gets large enough.

Rectangle Rule (Midpoint Rule)

The Midpoint Rule approximates the integrand as a constant (degree-0) polynomial on each subinterval, using the function value at the midpoint as the rectangle's height.

• The composite error is $O(h^2)$, where $h$ is the subinterval width. This matches the Trapezoidal Rule despite the simpler construction, because evaluating at the midpoint causes the first-order error term to cancel by symmetry.
• Left-endpoint and right-endpoint rectangle rules only achieve $O(h)$, so the midpoint choice is doing real work here.

Trapezoidal Rule

The Trapezoidal Rule uses linear (degree-1) interpolation, connecting function values at subinterval endpoints to form trapezoids.

• The composite error is $O(h^2)$ with leading error term $-\frac{(b-a)}{12}h^2 f''(\xi)$ for some $\xi$ in the interval. The second derivative controls accuracy: low curvature means small error.
• It's exact for linear functions (degree $\leq 1$) and performs well on smooth, non-oscillatory integrands.

Simpson's Rule

Simpson's Rule fits a quadratic (degree-2) polynomial through three equally spaced points (both endpoints and the midpoint of each pair of subintervals). This requires an even number of subintervals.

• The composite error is $O(h^4)$, two orders better than the Trapezoidal Rule. For smooth functions, this improvement is dramatic.
• It's exact for polynomials up to degree 3, not just degree 2. This "bonus" degree of exactness comes from the symmetry of the nodes and weights: the error term involving $f'''$ vanishes because the odd-powered error contribution cancels over symmetric intervals.
• The weighting pattern across nodes follows the familiar (1, 4, 2, 4, 2, ..., 4, 1) scheme in the composite version.

Newton-Cotes Formulas

The Newton-Cotes family is the general framework using equally spaced nodes. The Rectangle, Trapezoidal, and Simpson's Rules are all special cases (closed formulas using degree-0, 1, and 2 interpolation, respectively).

• Higher-degree Newton-Cotes formulas exist (e.g., Simpson's 3/8 rule at degree 3, Boole's rule at degree 4), but beyond moderate degree they suffer from Runge's phenomenon: oscillatory behavior in the interpolating polynomial can cause errors to grow with more nodes.
• Additionally, high-order closed Newton-Cotes formulas develop negative weights, which leads to numerical instability through cancellation.
• Open vs. closed variants: closed formulas include the interval endpoints as nodes; open formulas (like the Midpoint Rule) exclude them. Open formulas are useful when the integrand is undefined or singular at an endpoint.

---

Extrapolation and Refinement Techniques

These methods build on basic rules by systematically combining approximations to cancel error terms. The core principle: if you know the structure of the error expansion, you can eliminate leading terms through careful linear combinations.

Romberg Integration

Romberg Integration applies Richardson extrapolation to the composite Trapezoidal Rule. It exploits the fact that the Trapezoidal Rule's error has an expansion in even powers of $h$:

$T(h) = I + c_1 h^2 + c_2 h^4 + c_3 h^6 + \cdots$

By computing $T(h)$ for successively halved step sizes ($h, h/2, h/4, \ldots$), you build a triangular table where each new column eliminates another power of $h^2$ from the error.

• The first extrapolation column reproduces Simpson's Rule.
• Subsequent columns achieve even higher orders of accuracy.
• For smooth functions, the diagonal entries converge rapidly. But the method assumes the error expansion in even powers of $h$ is valid, which requires sufficient smoothness of the integrand.

Composite Integration Methods

Composite methods subdivide the full interval $[a, b]$ and apply a basic rule to each subinterval (or pair of subintervals for Simpson's).

• Composite Trapezoidal has global error $O(h^2)$; Composite Simpson's has global error $O(h^4)$.
• Single-interval formulas rarely provide sufficient accuracy on their own. Composition is how these methods are actually used in practice.
• The trade-off is straightforward: more subintervals means more function evaluations but smaller error, with the rate of improvement governed by the method's order.

---

Optimal Node Selection Methods

Rather than using equally spaced points, these methods choose node locations strategically to maximize accuracy per function evaluation. Optimal node placement can achieve exponential convergence for analytic functions, far outperforming fixed-spacing methods.

Gaussian Quadrature

Gaussian quadrature selects nodes as roots of orthogonal polynomials (Legendre polynomials for standard intervals, Chebyshev, Laguerre, or Hermite for other weight functions) along with corresponding optimal weights.

• An $n$-point Gaussian quadrature rule is exact for polynomials of degree $2n - 1$. That's twice the degree you'd expect from $n$ points, because both the node locations and the weights are free parameters ($2n$ parameters total).
• For analytic (infinitely differentiable) functions, Gaussian quadrature achieves exponential convergence as $n$ increases, making it the gold standard when function evaluations are expensive.
• The main limitation: nodes are on a fixed reference interval (typically $[-1, 1]$), so you need a change of variables for other domains. Also, the nodes are not nested between different $n$-values, so you can't easily reuse function evaluations when increasing $n$.

Adaptive Quadrature

Adaptive quadrature dynamically refines subintervals based on local error estimates, concentrating computational effort where the integrand is hardest to approximate.

The typical procedure:

1. Apply a quadrature rule (e.g., Simpson's) on a subinterval to get an estimate $S_1$.
2. Subdivide that interval in half and apply the same rule to each half, getting a refined estimate $S_2$.
3. Compare $|S_2 - S_1|$ to a local error tolerance. If the difference is small enough, accept $S_2$. If not, recursively subdivide each half.

This approach is essential for functions with localized complexity (sharp peaks, rapid oscillations in a small region, near-singularities) where uniform spacing would waste evaluations on smooth regions.

---

Probabilistic and High-Dimensional Methods

When deterministic methods become impractical, especially in high dimensions, randomized approaches offer a fundamentally different strategy. Monte Carlo convergence depends on sample count, not dimension, which breaks the curse of dimensionality.

Monte Carlo Integration

Monte Carlo integration estimates integrals via random sampling:

$\int_\Omega f(\mathbf{x}) \, dV \approx V(\Omega) \cdot \frac{1}{N} \sum_{i=1}^{N} f(\mathbf{x}_i)$

where $\mathbf{x}_i$ are points sampled uniformly at random from the domain $\Omega$ and $V(\Omega)$ is the volume of that domain.

• The convergence rate is $O(N^{-1/2})$ regardless of dimension. In 1D, this is terrible compared to Simpson's $O(N^{-4})$ or Gaussian quadrature's exponential rate. But in high dimensions, deterministic methods require a number of nodes that grows exponentially with dimension (the curse of dimensionality), while Monte Carlo's rate stays fixed.
• Variance reduction techniques can significantly improve practical performance:
  • Importance sampling: sample more densely where $|f|$ is large
  • Stratified sampling: divide the domain into strata and sample each
  • Control variates: subtract a function with a known integral to reduce variance

---

Error Analysis and Convergence Theory

Understanding error behavior is how you choose methods, set parameters, and verify results. Error expansions reveal both the rate of convergence and the conditions under which methods succeed or fail.

Error Analysis and Convergence Rates

• Order of accuracy $O(h^p)$ means the error scales as the $p$-th power of the step size. Halving $h$ (doubling the number of subintervals) reduces the error by a factor of $2^p$. This is a practical tool: if you compute results at two different step sizes, you can estimate $p$ empirically and check whether your implementation matches the theoretical order.
• Error bounds depend on derivative magnitudes: the Trapezoidal error bound involves $\max|f''|$, Simpson's involves $\max|f^{(4)}|$. This means higher-order methods only deliver their promised convergence rate if the integrand has enough smoothness (continuous derivatives of sufficiently high order).
• Asymptotic vs. pre-asymptotic behavior: theoretical rates like $O(h^4)$ apply only for sufficiently small $h$. In practice, if $h$ isn't small enough relative to the integrand's features, you may not yet see the expected convergence rate. Conversely, once $h$ is very small, round-off error from floating-point arithmetic can start to dominate, and further refinement actually increases the total error.

---

Quick Reference Table

---

Self-Check Questions

1. Both the Midpoint Rule and Trapezoidal Rule have $O(h^2)$ error. What geometric or symmetry argument explains why the Midpoint Rule achieves the same order as the Trapezoidal Rule despite using degree-0 (rather than degree-1) interpolation?

2. Simpson's Rule is exact for polynomials up to degree 3, not just degree 2. What property of the method causes this "bonus" degree of exactness?

3. Compare Romberg Integration and Adaptive Quadrature: under what conditions would you prefer each, and what assumption does Romberg make that Adaptive does not?

4. A colleague suggests using 10-point Gaussian quadrature for a 15-dimensional integral. Explain why this is impractical and what method you would recommend instead.

5. You're integrating a function with a sharp spike near $x = 0.3$ but smooth behavior elsewhere. Would Composite Simpson's with uniform subintervals or Adaptive Quadrature be more efficient? Justify your answer in terms of how each method distributes its error budget.