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.

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

Frequently Asked Questions