Trapezoidal Rule — $2M Options Pricing Failure

The Gaussian normal distribution has no closed-form integral — there is no formula for erf(x). Machine learning models using Gaussian distributions (virtually all of them) require numerical integration constantly. Physics simulations, financial option pricing (Black-Scholes), and engineering stress analysis all reduce to integrals that cannot be solved analytically.

Numerical integration is where theory meets real engineering constraints: function evaluation is often expensive, so minimising the number of evaluations while achieving target precision is the practical problem. Simpson's rule achieves O(h^4) error versus Trapezoidal's O(h^2) — for the same number of function evaluations, you get roughly 10^4 more accuracy. That difference matters when each evaluation takes a second to compute.

Why Numerical Integration Methods Are Not Optional

Numerical integration methods approximate the definite integral of a function when an analytic solution is impossible or impractical. Instead of finding a closed-form antiderivative, they sample the function at discrete points and combine those values using a weighted sum. The core mechanic is replacing a continuous curve with a finite set of geometric shapes — rectangles, trapezoids, or polynomials — whose area can be computed exactly.

Every method introduces error proportional to the step size raised to some power. The Trapezoidal Rule, for example, uses linear interpolation between sample points, yielding O(h²) error for smooth functions. Halving the step size quadruples accuracy — but only if the function's second derivative is bounded. That assumption fails spectacularly near singularities or discontinuities, producing silent garbage.

Use numerical integration when the integrand is expensive to evaluate, comes from experimental data, or has no closed form. It's the backbone of pricing engines, physics simulations, and signal processing pipelines. Choosing the wrong method or step size doesn't just waste cycles — it can misprice a portfolio by millions.

Trapezoidal Rule

Approximate the area as a sum of trapezoids: h/2 × (f(x0) + 2f(x1) + 2f(x2) + ... + 2f(x_{n-1}) + f(xn)). Error: O(h²) per interval, O(h²) total.

The trapezoidal rule is the simplest numerical integration method. It approximates the integral by dividing the region into n strips of equal width h = (b-a)/n and summing the area of each trapezoid. The error term for a single interval is - (h³/12) f''(ξ), meaning the method is exact for linear functions (second derivative zero). For smooth functions, doubling n reduces the error by a factor of 4.

def trapezoidal(f, a: float, b: float, n: int) -> float:
    """Integrate f from a to b using n trapezoids."""
    h = (b - a) / n
    total = f(a) + f(b)
    for i in range(1, n):
        total += 2 * f(a + i * h)
    return total * h / 2

import math

# Integrate sin(x) from 0 to pi = 2.0 (exact)
for n in [10, 100, 1000]:
    result = trapezoidal(math.sin, 0, math.pi, n)
    print(f'n={n:5}: {result:.10f}, error={abs(result-2):.2e}')

Simpson's Rule

Fit a parabola through each three consecutive points. Requires even n. h/3 × (f(x0) + 4f(x1) + 2f(x2) + 4f(x3) + ... + 4f(x_{n-1}) + f(xn)) Error: O(h⁴) — much better than Trapezoidal O(h²).

Simpson's rule uses quadratic interpolation over pairs of intervals. The pattern "4, 2, 4, 2, ..., 4" alternates coefficients. It is exact for polynomials up to degree 3 (cubic polynomials), because the error term involves the fourth derivative: - (h⁵/90) f^{(4)}(ξ). For smooth functions where the fourth derivative is bounded, Simpson's rule converges incredibly fast — notice with n=10, error is already 1e-7 for sin(x).

def simpsons(f, a: float, b: float, n: int) -> float:
    """Simpson's rule. n must be even."""
    if n % 2 != 0:
        n += 1
    h = (b - a) / n
    total = f(a) + f(b)
    for i in range(1, n):
        coeff = 4 if i % 2 != 0 else 2
        total += coeff * f(a + i * h)
    return total * h / 3

for n in [10, 100, 1000]:
    result = simpsons(math.sin, 0, math.pi, n)
    print(f'n={n:5}: {result:.10f}, error={abs(result-2):.2e}')

print('\nSimpson vs Trapezoid comparison (n=10):')
print(f'Trapezoid error: {abs(trapezoidal(math.sin, 0, math.pi, 10)-2):.2e}')
print(f'Simpson error:   {abs(simpsons(math.sin, 0, math.pi, 10)-2):.2e}')

Error Analysis

Simpson's with n=10 is more accurate than Trapezoidal with n=10,000 — same number of function evaluations, 1000x better.

The error analysis tells you more than just convergence rates. The constants matter: for trapezoidal, error ≈ (b-a)h²/12 max|f''|. For Simpson's, error ≈ (b-a)h⁴/180 max|f^{(4)}|. If your function has large second derivative but small fourth derivative, Simpson's wins dramatically. But if your function has huge fourth derivative (e.g., a sharp spike), Simpson's may actually be worse than a well-chosen trapezoidal with refined step near the spike.

Adaptive Integration

Fixed n doesn't work for all functions — some regions need more resolution. Adaptive integration refines intervals where the error is large.

Adaptive Simpson's method: compute the integral over [a,b] using Simpson's rule. Then split into two halves and apply Simpson's on each. The difference between the whole and the sum of halves gives an estimate of the error. If the error is larger than tolerance, recursively split each half. This concentrates evaluations where they're needed.

def adaptive_simpsons(f, a, b, tol=1e-10, depth=0, max_depth=50):
    """Adaptive Simpson's method."""
    c = (a + b) / 2
    fa, fb, fc = f(a), f(b), f(c)
    whole = (b-a)/6 * (fa + 4*fc + fb)
    left  = (c-a)/6 * (fa + 4*f((a+c)/2) + fc)
    right = (b-c)/6 * (fc + 4*f((c+b)/2) + fb)
    delta = left + right - whole
    if depth >= max_depth or abs(delta) < 15*tol:
        return left + right + delta/15
    return (adaptive_simpsons(f, a, c, tol/2, depth+1, max_depth) +
            adaptive_simpsons(f, c, b, tol/2, depth+1, max_depth))

# Integrate e^(-x^2) from 0 to 1 (no closed form)
result = adaptive_simpsons(lambda x: math.exp(-x**2), 0, 1)
print(f'∫e^(-x²)dx [0,1] = {result:.12f}')  # ≈ 0.746824132812
# Compare with scipy
from scipy import integrate
quad_result, _ = integrate.quad(lambda x: math.exp(-x**2), 0, 1)
print(f'scipy.quad         = {quad_result:.12f}')

Romberg Integration

Romberg integration uses Richardson extrapolation on the trapezoidal rule. By computing the trapezoidal rule for a sequence of step sizes (h, h/2, h/4, ...) and combining them, you can eliminate lower-order error terms and achieve very high accuracy with relatively few evaluations.

The method builds a triangular array called Romberg's table. The first column is the trapezoidal rule with decreasing step sizes. Each subsequent column applies Richardson extrapolation to remove error terms of increasing order. For sufficiently smooth functions, Romberg integration can achieve machine epsilon with fewer function evaluations than adaptive Simpson's.

def romberg(f, a, b, n=5):
    """Romberg integration with n steps of extrapolation."""
    h = b - a
    R = [[0.5 * h * (f(a) + f(b))]]
    for i in range(1, n+1):
        h /= 2
        # trapezoidal with step h
        total = 0.0
        x = a + h
        while x < b:
            total += f(x)
            x += h
        R_i0 = 0.5 * R[i-1][0] + total * h
        row = [R_i0]
        for j in range(1, i+1):
            extrap = row[j-1] + (row[j-1] - R[i-1][j-1]) / ((1 << (2*j)) - 1)
            row.append(extrap)
        R.append(row)
    return R[n][n]

import math
result = romberg(math.sin, 0, math.pi, n=5)
print(f'Romberg n=5: {result:.15f}, error={abs(result-2):.2e}')
# Compare with exact
print(f'Exact:        {2:.15f}')

When to Use Which Method

Choosing the right method depends on the integrand's behaviour, accuracy requirements, and computational budget.

High-level guidance: - Trapezoidal rule: use for quick estimates, when function evaluations are very cheap, or when the function is only piecewise smooth and you can manually split intervals. - Simpson's rule: use for smooth, well-behaved functions where you need moderate accuracy (e.g., engineering approximations). It's the workhorse of numerical integration. - Adaptive Simpson's: use when you don't know the function's behaviour a priori — it's robust and widely used in production systems. - Romberg integration: use for extremely smooth functions where you need very high accuracy with minimal evaluations (e.g., reference solutions for benchmarks). - Gaussian quadrature: use when the function is smooth and you can afford to precompute nodes and weights. Gauss-Legendre quadrature with n points integrates exactly polynomials of degree 2n-1. For many applications, n=5 gives machine epsilon.

Production recommendation: Default to scipy.integrate.quad (adaptive Gauss-Kronrod). It handles discontinuities, infinite limits, and oscillatory integrands better than any simple method.

import numpy as np
from scipy import integrate

def decide_method(f, a, b, max_evals=1000):
    """Return recommended method based on simple tests."""
    # Quick check for smoothness via finite difference
    xs = np.linspace(a, b, 10)
    try:
        deriv4 = np.gradient(np.gradient(np.gradient(np.gradient(f(xs)))))
        if np.max(np.abs(deriv4)) < 10:
            return "Romberg or Simpson's (fixed n)"
        else:
            # Check for rapid oscillation
            zero_crossings = np.sum(np.diff(np.sign(f(xs))) != 0)
            if zero_crossings > 3:
                return "Quad (adaptive Gauss-Kronrod)"
            else:
                return "Adaptive Simpson's"
    except:
        return "Quad (adaptive Gauss-Kronrod)"

# Example
f = lambda x: np.exp(-x**2)
print(decide_method(f, 0, 1))

Gauss-Legendre Quadrature: When You Need Bloody Accuracy Without Bloody Intervals

Trapezoidal and Simpson's use evenly spaced points. That's fine if you're teaching undergrads. In production you want maximal accuracy per function evaluation. Gauss-Legendre quadrature picks the optimal abscissas — the roots of Legendre polynomials — and weights them to integrate polynomials of degree 2n-1 exactly with only n points.

Why does that matter? Because every function call to your simulation costs CPU time. Your 3D fluid solver calls this a million times per timestep. Chebyshev nodes help fight Runge's phenomenon, but Gauss-Legendre gives you spectral accuracy for smooth integrands with half the evaluations.

You still need to transform your integration limits [a, b] into [-1, 1]. That's a linear mapping, trivial code. The abscissas and weights are precomputed for given n — grab them from a table or compute on startup. For n=5 you integrate a 9th-degree polynomial exactly. Compare that to Trapezoidal's first-order convergence and weep.

// io.thecodeforge — dsa tutorial

public class GaussLegendreQuad {
    // Precomputed abscissas and weights for n=5 (degree 9 exact)
    private static final double[] NODES = {
        -0.9061798459386640, -0.5384693101056831,
         0.0, 0.5384693101056831, 0.9061798459386640
    };
    private static final double[] WEIGHTS = {
        0.2369268850561891, 0.4786286704993665,
        0.5688888888888889, 0.4786286704993665,
        0.2369268850561891
    };

    public static double integrate(IntegrationFunction f, double a, double b) {
        double half = (b - a) * 0.5;
        double mid  = (b + a) * 0.5;
        double sum = 0.0;
        for (int i = 0; i < NODES.length; i++) {
            double x = mid + half * NODES[i];
            sum += WEIGHTS[i] * f.evaluate(x);
        }
        return half * sum;
    }
}

Clenshaw-Curtis: Chebyshev Nodes for the Practical Engineer

Gauss-Legendre is the gold standard, but its nodes are irrational. Clenshaw-Curtis quadrature uses Chebyshev extrema — x_k = cos(kπ/n). That's cheaper to compute and gives you nested grids for adaptive refinement. The weights come from a cosine transform, which you can slam through an FFT.

Why use it? When your integrand has edge singularities or rapid oscillations, Chebyshev clustering at boundaries fights Runge's phenomenon better than Gauss-Legendre's open nodes. Engineers in signal processing and filter design reach for this first because they already have FFT libraries.

You lose a factor of two in polynomial exactness compared to Gauss-Legendre for the same n, but you gain nesting and FFT speed. For n=1024, Clenshaw-Curtis with an O(n log n) FFT beats Gauss-Legendre's O(n^2) node computation cold. The convergence is still exponential for analytic functions — just slightly slower.

Rule of thumb: if you're doing fixed-order quadrature and your function calls are the bottleneck, go Gauss-Legendre. If you're adapting order or already have an FFT, Clenshaw-Curtis wins.

// io.thecodeforge — dsa tutorial

import java.util.function.DoubleUnaryOperator;

public class ClenshawCurtisDemo {
    // Compute Clenshaw-Curtis weights via FFT (simplified)
    public static double integrate(DoubleUnaryOperator f, int n) {
        // n+1 Chebyshev nodes: x_j = cos(j*pi/n)
        double sum = 0.5 * (f.applyAsDouble(1.0) + f.applyAsDouble(-1.0));
        for (int j = 1; j < n; j++) {
            double x = Math.cos(j * Math.PI / n);
            double fxj = f.applyAsDouble(x);
            // weight via Clenshaw-Curtis recurrence (not FFT for brevity)
            sum += fxj * computeWeight(j, n);
        }
        return sum * Math.PI / n;
    }

    private static double computeWeight(int j, int n) {
        // Full implementation uses FFT — this is placeholder
        return 1.0; 
    }
}