Numerical Integration

Numerical Integration

Definitions

Formulas

Trapezoidal Rule

Derivation

From the single-interval trapezoidal rule applied to each subinterval and summed:

• Total = Σᵢ₌₀ⁿ⁻¹ h·(f(xᵢ) + f(xᵢ₊₁))/2
• Interior points appear twice (once as right endpoint, once as left)
• Endpoints appear once
• Result: h/2·[f(x₀) + 2f(x₁) + ... + 2f(x₋₁) + f(xₙ)]

Simpson's Rule

Weight Pattern

For n = 6 (7 nodes): weights are [1, 4, 2, 4, 2, 4, 1] × h/3 Pattern: 1, 4, 2, 4, 2, 4, ..., 4, 1

Gaussian Quadrature

Legendre Polynomial Nodes and Weights

Transformation

To integrate over [a, b] instead of [−1, 1]: x = ((b−a)t + (b+a))/2, dx = (b−a)/2 dt

Examples

import numpy as np


def trapezoidal(f, a, b, n=1000):
    """Composite trapezoidal rule."""
    x = np.linspace(a, b, n + 1)
    y = f(x)
    h = (b - a) / n
    return h * (y[0] / 2 + np.sum(y[1:-1]) + y[-1] / 2)


def simpson(f, a, b, n=1000):
    """Composite Simpson's rule (n must be even)."""
    if n % 2 != 0:
        n += 1
    x = np.linspace(a, b, n + 1)
    y = f(x)
    h = (b - a) / n
    return h / 3 * (y[0] + 4 * np.sum(y[1::2]) + 2 * np.sum(y[2:-1:2]) + y[-1])


def gauss_legendre(f, a, b, n=5):
    """Gaussian quadrature using numpy's Gauss-Legendre nodes."""
    nodes, weights = np.polynomial.legendre.leggauss(n)
    # Transform from [-1, 1] to [a, b]
    t = 0.5 * (b - a) * nodes + 0.5 * (b + a)
    w = 0.5 * (b - a) * weights
    return np.sum(w * f(t))


# Test function
f = np.exp
exact = np.e - 1

print("Comparing Integration Methods for ∫₀¹ eˣ dx")
print(f"Exact value: {exact:.12f}")
print(f"{'Method':<20} {'n':<6} {'Result':<16} {'Error':<12}")
print("-" * 56)

for n in [4, 8, 16, 32, 64]:
    trap = trapezoidal(f, 0, 1, n)
    simp = simpson(f, 0, 1, n)
    gauss = gauss_legendre(f, 0, 1, n // 2 + 1)

    print(f"{'Trapezoidal':<20} {n:<6} {trap:<16.12f} {abs(trap - exact):<12.2e}")
    print(f"{'Simpson':<20} {n:<6} {simp:<16.12f} {abs(simp - exact):<12.2e}")
    print(f"{'Gauss-Legendre':<20} {n:<6} {gauss:<16.12f} {abs(gauss - exact):<12.2e}")
    print()

import numpy as np
from scipy.integrate import quad


def integrand(x):
    """A function with a sharp peak — challenges uniform methods."""
    return np.exp(-x**2) * np.sin(10 * x) ** 2


# Adaptive quadrature
result, error = quad(integrand, -5, 5)
print(f"Adaptive quadrature: {result:.12f}")
print(f"Estimated error:     {error:.2e}")
print(f"Actual evaluations:  likely ~50-200 (adaptive)")

# Compare with fixed methods
print(f"\nFixed n=100:")
from numpy import sin, exp

def trapezoidal(f, a, b, n):
    x = np.linspace(a, b, n + 1)
    y = f(x)
    h = (b - a) / n
    return h * (y[0]/2 + np.sum(y[1:-1]) + y[-1]/2)

trap_100 = trapezoidal(integrand, -5, 5, 100)
print(f"  Trapezoidal (n=100): {trap_100:.12f}")
print(f"  Error:               {abs(trap_100 - result):.2e}")

# Oscillatory integral
def oscillatory(x):
    return np.sin(100 * x) / (1 + x**2)

result_osc, error_osc = quad(oscillatory, 0, 10)
print(f"\nOscillatory integral:  {result_osc:.8f}")
print(f"Estimated error:      {error_osc:.2e}")

import numpy as np


def monte_carlo_circle(n_samples):
    """Estimate π using Monte Carlo integration."""
    x = np.random.uniform(-1, 1, n_samples)
    y = np.random.uniform(-1, 1, n_samples)

    # Count points inside the unit circle
    inside = x**2 + y**2 <= 1.0
    pi_estimate = 4.0 * np.sum(inside) / n_samples

    return pi_estimate


# Demonstrate convergence
print("Monte Carlo Estimation of π")
print(f"{'Samples':<12} {'Estimate':<14} {'Error':<12}")
print("-" * 38)

for n in [100, 1_000, 10_000, 100_000, 1_000_000]:
    np.random.seed(42)
    est = monte_carlo_circle(n)
    print(f"{n:<12,} {est:<14.8f} {abs(est - np.pi):<12.6f}")

# Verify O(1/√N) convergence
print(f"\nExpected error reduction for 10x samples: {1/np.sqrt(10):.4f}")

import numpy as np


# Simulated experimental data: temperature measurements over time
time = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])  # seconds
temp = np.array([20, 22, 25, 28, 31, 33, 34, 34, 33, 31, 28])  # °C

# Trapezoidal rule on discrete data
dt = time[1] - time[0]  # Uniform spacing
area_trap = dt * (temp[0]/2 + np.sum(temp[1:-1]) + temp[-1]/2)
print(f"Trapezoidal (manual):  {area_trap:.2f} °C·s")

# Using numpy
area_np = np.trapz(temp, time)
print(f"numpy.trapz:           {area_np:.2f} °C·s")

# Simpson's rule (n must be even — we have 11 points, 10 intervals)
n = len(time) - 1  # 10 intervals
h = (time[-1] - time[0]) / n
area_simp = h/3 * (temp[0] + 4*np.sum(temp[1::2]) + 2*np.sum(temp[2:-1:2]) + temp[-1])
print(f"Simpson's rule:        {area_simp:.2f} °C·s")

# Physical interpretation
print(f"\nMean temperature: {area_trap / (time[-1] - time[0]):.2f} °C")

Theorems

Python Implementation

import numpy as np
from typing import Callable, Tuple


class NumericalIntegrator:
    """Comprehensive numerical integration toolkit."""

    @staticmethod
    def trapezoidal(f: Callable, a: float, b: float,
                    n: int = 1000) -> Tuple[float, float]:
        """Composite trapezoidal rule. Returns (result, error_estimate)."""
        x = np.linspace(a, b, n + 1)
        y = f(x)
        h = (b - a) / n
        result = h * (y[0] / 2 + np.sum(y[1:-1]) + y[-1] / 2)

        # Error estimate via Richardson extrapolation
        x2 = np.linspace(a, b, 2 * n + 1)
        y2 = f(x2)
        h2 = h / 2
        result2 = h2 * (y2[0] / 2 + np.sum(y2[1:-1]) + y2[-1] / 2)

        error = abs(result2 - result) / 3  # Richardson estimate
        return result, error

    @staticmethod
    def simpson(f: Callable, a: float, b: float,
                n: int = 1000) -> Tuple[float, float]:
        """Composite Simpson's rule (n forced even). Returns (result, error)."""
        if n % 2 != 0:
            n += 1
        x = np.linspace(a, b, n + 1)
        y = f(x)
        h = (b - a) / n
        result = h / 3 * (y[0] + 4 * np.sum(y[1::2]) +
                          2 * np.sum(y[2:-1:2]) + y[-1])

        # Error estimate
        n2 = 2 * n
        x2 = np.linspace(a, b, n2 + 1)
        y2 = f(x2)
        h2 = (b - a) / n2
        result2 = h2 / 3 * (y2[0] + 4 * np.sum(y2[1::2]) +
                            2 * np.sum(y2[2:-1:2]) + y2[-1])

        error = abs(result2 - result) / 15  # Richardson estimate
        return result, error

    @staticmethod
    def gauss_legendre(f: Callable, a: float, b: float,
                       n: int = 5) -> float:
        """n-point Gauss-Legendre quadrature."""
        nodes, weights = np.polynomial.legendre.leggauss(n)
        t = 0.5 * (b - a) * nodes + 0.5 * (b + a)
        w = 0.5 * (b - a) * weights
        return np.sum(w * f(t))

    @staticmethod
    def romberg(f: Callable, a: float, b: float,
                max_level: int = 10) -> np.ndarray:
        """Romberg integration using Richardson extrapolation."""
        R = np.zeros((max_level + 1, max_level + 1))

        # Trapezoidal estimates at various levels
        n = 1
        h = b - a
        R[0, 0] = h / 2 * (f(a) + f(b))

        for i in range(1, max_level + 1):
            n = 2 ** i
            h = (b - a) / n
            # Add new midpoints
            x_new = np.arange(1, n, 2) * h + a
            R[i, 0] = 0.5 * R[i-1, 0] + h * np.sum(f(x_new))

            # Richardson extrapolation
            for j in range(1, i + 1):
                R[i, j] = R[i, j-1] + (R[i, j-1] - R[i-1, j-1]) / (4**j - 1)

        return R

    @staticmethod
    def monte_carlo(f: Callable, domain: list,
                    n_samples: int = 100000) -> Tuple[float, float]:
        """Monte Carlo integration over a rectangular domain."""
        a, b = domain[0], domain[1]
        volume = b - a

        x = np.random.uniform(a, b, n_samples)
        values = f(x)
        result = volume * np.mean(values)
        std_err = volume * np.std(values) / np.sqrt(n_samples)

        return result, std_err


def demonstrate_integration():
    """Full demonstration of all integration methods."""
    print("=" * 60)
    print("NUMERICAL INTEGRATION METHOD COMPARISON")
    print("=" * 60)

    # Test 1: Smooth function
    f1 = lambda x: np.exp(-x**2)
    exact1 = np.sqrt(np.pi)  # ∫₋∞^∞ e^(-x²) dx = √π
    # But we integrate from -5 to 5 (close enough)

    from scipy.integrate import quad
    exact1_quad, _ = quad(f1, -5, 5)

    print(f"\n∫₋₅⁵ e^(-x²) dx")
    print(f"scipy quad: {exact1_quad:.12f}")

    integrator = NumericalIntegrator()

    for n in [10, 50, 200]:
        trap, err = integrator.trapezoidal(f1, -5, 5, n)
        simp, err = integrator.simpson(f1, -5, 5, n)
        print(f"  n={n:>3}: Trapezoidal {trap:.10f} (err {err:.2e}) | "
              f"Simpson {simp:.10f} (err {err:.2e})")

    # Gauss-Legendre
    for n in [3, 5, 10]:
        gl = integrator.gauss_legendre(f1, -5, 5, n)
        print(f"  Gauss-Legendre (n={n}): {gl:.10f}")

    # Romberg
    R = integrator.romberg(f1, -5, 5, max_level=6)
    print(f"\nRomberg integration:")
    for i in range(min(7, R.shape[0])):
        print(f"  Level {i}: {R[i, i]:.12f}")

    # Test 2: Oscillatory function
    print(f"\n{'=' * 60}")
    f2 = lambda x: np.sin(100 * x) / (1 + x**2)
    mc_result, mc_err = integrator.monte_carlo(f2, [-10, 10], 1_000_000)
    quad_result, _ = quad(f2, -10, 10)
    print(f"∫₋₁₀¹⁰ sin(100x)/(1+x²) dx")
    print(f"scipy quad:  {quad_result:.10f}")
    print(f"Monte Carlo: {mc_result:.10f} (±{mc_err:.2e})")


demonstrate_integration()

Applications in AI/ML

Variational Inference

In Bayesian neural networks, the posterior P(θ|D) is intractable. Variational inference approximates it with a tractable distribution q(θ) and minimizes KL(q||p), which requires computing integrals like E_q[log p(D|θ)] using quadrature or Monte Carlo.

Monte Carlo in ML

Loss Function Integration

Some loss functions require integration: Expected Calibration Error (ECE) integrates over confidence bins, and proper scoring rules integrate over predicted distributions.

Common Mistakes

Interview Questions

Q1: Why is Gaussian quadrature more accurate than equally-spaced rules? A: Gaussian quadrature optimizes both node placement and weights to maximize the degree of exactness. With n nodes it integrates polynomials of degree 2n−1 exactly, while equally-spaced rules (Newton-Cotes) with n nodes achieve at most degree n (or n+1 for Simpson's). The trade-off is that Gaussian nodes are irrational, requiring more complex implementation.

Q2: When would you use Monte Carlo over deterministic quadrature? A: Monte Carlo is preferred for high-dimensional integrals (d > 4), where deterministic methods suffer the curse of dimensionality. It is also useful when the integrand is noisy, discontinuous, or only available through simulations. For low-dimensional smooth functions, deterministic methods are far more efficient.

Q3: How does adaptive quadrature decide where to subdivide? A: Adaptive quadrature estimates the local error on each subinterval. If the error exceeds the tolerance, the interval is split (usually in half). The process recurses until all subintervals meet the error tolerance. This concentrates evaluations where the integrand varies rapidly.

Q4: Explain the relationship between numerical integration and numerical differentiation. A: Numerical differentiation amplifies noise (it's an ill-conditioned problem), while numerical integration smooths it (it's well-conditioned). The error in numerical differentiation is O(ε/h) from roundoff plus O(h^p) from truncation, creating an optimal h ≈ √ε. Integration error simply decreases with h.

Q5: What is Richardson extrapolation and how is it used in Romberg integration? A: Richardson extrapolation combines numerical results at different step sizes to eliminate leading error terms. If R(h) = I + c₁h^p + ..., then combining R(h) and R(h/2) eliminates the h^p term, achieving higher accuracy. Romberg integration applies this recursively to trapezoidal rule estimates, achieving exponential convergence for smooth functions.

Q6: How do you handle integrable singularities? A: Strategies include: (1) variable substitution to remove the singularity (e.g., t = √x for x^(-1/2) singularity), (2) splitting the integral and handling the singular part analytically, (3) using Gauss-Jacobi or other specialized quadrature designed for singular weight functions, (4) adaptive quadrature which may concentrate points near the singularity.

Practice Problems

import numpy as np
from sympy import symbols, integrate, Poly

# Symbolic derivation
x = symbols('x')

# Fit parabola through (-h, y0), (0, y1), (h, y2)
# Lagrange basis:
# L0 = (x)(x-h) / ((-h)(-2h)) = x(x-h)/(2h²)
# L1 = (x+h)(x-h) / ((h)(-h)) = -(x²-h²)/h²
# L2 = (x)(x+h) / ((2h)(h)) = x(x+h)/(2h²)
# P(x) = y0*L0 + y1*L1 + y2*L2

h_sym = symbols('h', positive=True)
y0, y1, y2 = symbols('y0 y1 y2')

L0 = x * (x - h_sym) / (2 * h_sym**2)
L1 = -(x**2 - h_sym**2) / h_sym**2
L2 = x * (x + h_sym) / (2 * h_sym**2)

P = y0 * L0 + y1 * L1 + y2 * L2

# Integrate from -h to h
result = integrate(P, (x, -h_sym, h_sym))
print(f"∫ P(x) dx from -h to h = {result}")
print(f"Simplified: {result.simplify()}")
# Result: h/3 * (y0 + 4*y1 + y2) — Simpson's 1/3 rule

import numpy as np


def monte_carlo_2d(f, n_samples):
    """Standard Monte Carlo integration in 2D."""
    x = np.random.uniform(0, 1, n_samples)
    y = np.random.uniform(0, 1, n_samples)
    return np.mean(f(x, y))


def pseudo_random_2d(f, n_samples):
    """Quasi-Monte Carlo using scrambled Sobol-like sequence."""
    # Simple low-discrepancy: nested uniform grid with random shift
    n_side = int(np.ceil(np.sqrt(n_samples)))
    grid = np.linspace(0, 1, n_side, endpoint=False)
    xx, yy = np.meshgrid(grid, grid)
    x = (xx.ravel()[:n_samples] + np.random.uniform(0, 1/n_side)) % 1
    y = (yy.ravel()[:n_samples] + np.random.uniform(0, 1/n_side)) % 1
    return np.mean(f(x, y))


# Test function
f = lambda x, y: np.sin(np.pi * x) * np.cos(np.pi * y)
exact = (2/np.pi) * 0  # ∫₀¹ sin(πx)dx × ∫₀¹ cos(πy)dy = (2/π) × 0 = 0
# Actually: ∫₀¹ sin(πx) dx = 2/π, ∫₀¹ cos(πy) dy = 0
exact = 0  # So this is a bad test... let me use a different one

f = lambda x, y: x * y + np.sin(2 * np.pi * x)
exact = 0.25  # ∫₀¹∫₀¹ (xy + sin(2πx)) dxdy = 1/4 + 0 = 1/4

print("Monte Carlo vs Quasi-Monte Carlo Convergence")
print(f"{'N':<12} {'MC Error':<14} {'QMC Error':<14}")
print("-" * 40)

for n in [100, 1_000, 10_000, 100_000]:
    np.random.seed(42)
    mc_errors = [abs(monte_carlo_2d(f, n) - exact) for _ in range(5)]
    qmc_errors = [abs(pseudo_random_2d(f, n) - exact) for _ in range(5)]

    print(f"{n:<12} {np.mean(mc_errors):<14.6f} {np.mean(qmc_errors):<14.6f}")

Quick Reference

Cross-References

• 086-Floating-Point — Numerical integration error is bounded by floating point precision; long summations need Kahan summation
• 087-Root-Finding — Adaptive quadrature uses error estimation similar to root-finding convergence criteria
• 089-Interpolation — Quadrature rules are derived by integrating interpolating polynomials; Gaussian quadrature uses optimal node placement
• 090-Error Analysis — Truncation error analysis determines the order of convergence for integration methods
• Monte Carlo Methods — Integration is the foundation of Monte Carlo simulation in physics, finance, and ML
• ODE Solvers: Every ODE solver (Euler, Runge-Kutta) uses quadrature internally to advance the solution