Numerical Integration — Theory Notebook

"Numerical integration is the art of turning hard integrals into easy sums."

Interactive derivations: Newton-Cotes rules, Gaussian quadrature, Richardson extrapolation, adaptive quadrature, Monte Carlo, and ML applications.

Companion: notes.md | exercises.ipynb

Code cell 2

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

try:
    import seaborn as sns
    sns.set_theme(style="whitegrid", palette="colorblind")
    HAS_SNS = True
except ImportError:
    plt.style.use("seaborn-v0_8-whitegrid")
    HAS_SNS = False

mpl.rcParams.update({
    "figure.figsize":    (10, 6),
    "figure.dpi":         120,
    "font.size":           13,
    "axes.titlesize":      15,
    "axes.labelsize":      13,
    "xtick.labelsize":     11,
    "ytick.labelsize":     11,
    "legend.fontsize":     11,
    "legend.framealpha":   0.85,
    "lines.linewidth":      2.0,
    "axes.spines.top":     False,
    "axes.spines.right":   False,
    "savefig.bbox":       "tight",
    "savefig.dpi":         150,
})
np.random.seed(42)
print("Plot setup complete.")

Code cell 3

import numpy as np
import scipy.linalg as la
from scipy import integrate
from scipy.fft import fft

try:
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    try:
        import seaborn as sns
        sns.set_theme(style='whitegrid', palette='colorblind')
        HAS_SNS = True
    except ImportError:
        plt.style.use('seaborn-v0_8-whitegrid')
        HAS_SNS = False
    mpl.rcParams.update({
        'figure.figsize': (10, 6), 'figure.dpi': 120,
        'font.size': 13, 'axes.titlesize': 15, 'axes.labelsize': 13,
        'xtick.labelsize': 11, 'ytick.labelsize': 11,
        'legend.fontsize': 11, 'lines.linewidth': 2.0,
        'axes.spines.top': False, 'axes.spines.right': False,
    })
    HAS_MPL = True
except ImportError:
    HAS_MPL = False

COLORS = {
    'primary':   '#0077BB',
    'secondary': '#EE7733',
    'tertiary':  '#009988',
    'error':     '#CC3311',
    'neutral':   '#555555',
    'highlight': '#EE3377',
}

np.set_printoptions(precision=8, suppress=True)
np.random.seed(42)
print('Setup complete.')

print("Chapter helper setup complete.")

1. Newton-Cotes Rules

Deriving and verifying Newton-Cotes rules: midpoint, trapezoidal, and Simpson's.

Code cell 5

# === 1.1 Newton-Cotes Rules: Basic Implementation ===

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

def composite_midpoint(f, a, b, m):
    """Composite midpoint rule with m intervals."""
    h = (b - a) / m
    x = a + h * (np.arange(m) + 0.5)  # Midpoints
    return h * f(x).sum()

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

# Test: integral of exp(x) from 0 to 1
f = np.exp
exact = np.e - 1  # = 1.718281828...

print('Convergence of Newton-Cotes rules:')
print(f'  {'m':>6}  {'Trap error':>12}  {'Midpt error':>12}  {'Simp error':>12}')
for m in [10, 100, 1000]:
    et = abs(composite_trap(f, 0, 1, m) - exact)
    em = abs(composite_midpoint(f, 0, 1, m) - exact)
    es = abs(composite_simpson(f, 0, 1, m//2*2) - exact)
    print(f'  {m:>6}  {et:>12.4e}  {em:>12.4e}  {es:>12.4e}')

print('\nMidpoint constant / Trapezoidal constant (should be ~0.5):')
print(f'  {abs(composite_midpoint(f,0,1,10)-exact) / abs(composite_trap(f,0,1,10)-exact):.4f}')

Code cell 6

# === 1.2 Error Scaling: O(h^2) vs O(h^4) ===

m_vals = [10, 20, 40, 80, 160, 320]
errs_trap = [abs(composite_trap(np.exp, 0, 1, m) - (np.e-1)) for m in m_vals]
errs_simp = [abs(composite_simpson(np.exp, 0, 1, m) - (np.e-1)) for m in m_vals]
hs = [1.0/m for m in m_vals]

print('Error ratios when m doubles (h halves):')
print(f'  {'m':>5}  {'h':>8}  {'Trap err':>12}  {'ratio':>6}  {'Simp err':>12}  {'ratio':>6}')
for i, (m, h, et, es) in enumerate(zip(m_vals, hs, errs_trap, errs_simp)):
    rt = errs_trap[i-1]/et if i > 0 else float('nan')
    rs = errs_simp[i-1]/es if i > 0 else float('nan')
    print(f'  {m:>5}  {h:>8.5f}  {et:>12.4e}  {rt:>6.2f}  {es:>12.4e}  {rs:>6.2f}')

print('\nExpected: Trapezoidal ratio → 4 (O(h^2)), Simpson → 16 (O(h^4))')

if HAS_MPL:
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.loglog(hs, errs_trap, 'o-', color=COLORS['error'], label='Trapezoidal')
    ax.loglog(hs, errs_simp, 's-', color=COLORS['primary'], label="Simpson's")
    h_ref = np.array(hs)
    ax.loglog(h_ref, h_ref**2 * errs_trap[0]/hs[0]**2, '--',
              color=COLORS['neutral'], alpha=0.7, label='O(h²) reference')
    ax.loglog(h_ref, h_ref**4 * errs_simp[0]/hs[0]**4, ':',
              color=COLORS['neutral'], alpha=0.7, label='O(h⁴) reference')
    ax.set_title('Newton-Cotes Convergence')
    ax.set_xlabel('Step size h'); ax.set_ylabel('Error')
    ax.legend(); plt.show()

2. Gaussian Quadrature

Deriving GL nodes via the Golub-Welsch algorithm and verifying the degree of exactness.

Code cell 8

# === 2.1 Gauss-Legendre via Golub-Welsch ===

def gauss_legendre(n):
    """
    Compute n-point Gauss-Legendre nodes and weights.
    Uses the Golub-Welsch algorithm (eigenvalue of Jacobi matrix).
    """
    k = np.arange(1, n)
    beta = k / np.sqrt(4*k**2 - 1)  # Off-diagonal elements
    J = np.diag(beta, -1) + np.diag(beta, 1)  # Symmetric tridiagonal
    eigvals, eigvecs = np.linalg.eigh(J)
    idx = np.argsort(eigvals)
    nodes = eigvals[idx]
    weights = 2 * eigvecs[0, idx]**2  # First row of Q
    return nodes, weights

# Verify GL nodes for n=3
n = 3
x, w = gauss_legendre(n)
print(f'GL-{n} nodes:   {x.round(8)}')
print(f'GL-{n} weights: {w.round(8)}')
print(f'Sum of weights: {w.sum():.8f} (should be 2 = integral of 1 over [-1,1])')

# Compare with known values
# n=3: nodes at 0 and ±sqrt(3/5), weights 8/9, 5/9, 5/9
known_nodes = np.array([-np.sqrt(3/5), 0, np.sqrt(3/5)])
known_weights = np.array([5/9, 8/9, 5/9])
print(f'\nExpected nodes:   {known_nodes.round(8)}')
print(f'Expected weights: {known_weights.round(8)}')
ok = np.allclose(np.sort(x), np.sort(known_nodes), atol=1e-12)
print(f'PASS: GL-3 matches known values' if ok else 'FAIL')

Code cell 9

# === 2.2 Degree of Exactness Verification ===

def gl_integrate(f, a, b, n):
    """Gauss-Legendre quadrature on [a,b] with n nodes."""
    x, w = gauss_legendre(n)
    # Map from [-1,1] to [a,b]
    x_mapped = (a+b)/2 + (b-a)/2 * x
    return (b-a)/2 * np.dot(w, f(x_mapped))

# GL-n should integrate polynomials of degree <= 2n-1 exactly
n = 5  # Test GL-5
print(f'GL-{n} degree of exactness verification:')
print(f'  Degree  Exact integral    GL estimate      Error')
for deg in range(0, 12):
    f_poly = lambda x, d=deg: x**d
    exact_val = (1 - (-1)**(deg+1)) / (deg+1)  # integral of x^d on [-1,1]
    gl_val = gl_integrate(f_poly, -1, 1, n)
    err = abs(gl_val - exact_val)
    status = 'EXACT' if err < 1e-10 else ''
    print(f'  {deg:>6}  {exact_val:>15.8f}  {gl_val:>15.8f}  {err:.2e} {status}')

print(f'\nGL-{n} is exact for degrees 0 to {2*n-1} (as predicted by theory)')

Code cell 10

# === 2.3 GL vs Composite Trapezoidal Efficiency ===

f_test = np.exp
exact = np.e - 1

print('Comparison: GL-n vs Composite Trapezoidal for same work:')
print(f'  {'Method':>30}  {'Evals':>6}  {'Error':>12}')

for n_gl in [2, 3, 5, 8, 10]:
    err = abs(gl_integrate(f_test, 0, 1, n_gl) - exact)
    print(f'  {'GL-'+str(n_gl):>30}  {n_gl:>6}  {err:>12.4e}')

for m in [10, 100, 1000, 10000]:
    err = abs(composite_trap(f_test, 0, 1, m) - exact)
    print(f'  {'Trapezoidal m='+str(m):>30}  {m+1:>6}  {err:>12.4e}')

print('\nGL-10 (10 evals) vs Trap-10000 (10001 evals):')
print(f'  GL-10 needs 10 evals; Trap needs 10001 evals for comparable accuracy')

3. Richardson Extrapolation and Romberg Integration

Code cell 12

# === 3.1 Romberg Integration ===

def romberg(f, a, b, max_k=8, tol=1e-12):
    """Romberg integration via repeated Richardson extrapolation."""
    R = np.zeros((max_k, max_k))
    h = b - a
    R[0, 0] = h/2 * (f(a) + f(b))

    for k in range(1, max_k):
        h /= 2
        n_new = 2**(k-1)
        x_new = a + h * (2*np.arange(1, n_new+1) - 1)
        R[k, 0] = R[k-1, 0]/2 + h * f(x_new).sum()
        for j in range(1, k+1):
            R[k, j] = (4**j * R[k, j-1] - R[k-1, j-1]) / (4**j - 1)
        if k > 0 and abs(R[k, k] - R[k-1, k-1]) < tol:
            return R[:k+1, :k+1], k
    return R, max_k - 1

# Test on e^x
exact = np.e - 1
R_table, final_k = romberg(np.exp, 0, 1, max_k=8)

print('Romberg table for integral of e^x from 0 to 1:')
print('Diagonal entries (main estimate at each level):')
for k in range(final_k + 1):
    evals = 2**k + 1
    err = abs(R_table[k, k] - exact)
    print(f'  k={k}: R[{k},{k}] = {R_table[k,k]:.12f}, error = {err:.2e}, '
          f'evals = {evals}')

print(f'\nExact: {exact:.12f}')
print(f'Romberg converged at k={final_k} with error {abs(R_table[final_k,final_k]-exact):.2e}')

Code cell 13

# === 3.2 Euler-Maclaurin: Periodic vs Non-Periodic ===

# Periodic function: e^{sin(2*pi*x)} on [0,1]
f_periodic = lambda x: np.exp(np.sin(2*np.pi*x))

# Non-periodic: e^x on [0,1]
f_nonperiodic = np.exp

# Exact integrals (via scipy)
exact_p, _ = integrate.quad(f_periodic, 0, 1)
exact_np = np.e - 1

print('Trapezoidal rule: periodic vs non-periodic functions')
print(f'  {'m':>5}  {'Periodic error':>16}  {'Non-periodic error':>20}')
for m in [5, 10, 20, 40]:
    ep = abs(composite_trap(f_periodic, 0, 1, m) - exact_p)
    enp = abs(composite_trap(f_nonperiodic, 0, 1, m) - exact_np)
    print(f'  {m:>5}  {ep:>16.4e}  {enp:>20.4e}')

print('\nPeriodic: spectral accuracy (error decreases exponentially)')
print('Non-periodic: O(h^2) convergence')

4. Monte Carlo Integration

Code cell 15

# === 4.1 Basic Monte Carlo ===

np.random.seed(42)

# Estimate pi/4 = Prob(X^2 + Y^2 <= 1) via MC
def estimate_pi(n):
    x, y = np.random.rand(n), np.random.rand(n)
    inside = (x**2 + y**2 <= 1)
    pi_est = 4 * inside.mean()
    std_err = 4 * np.sqrt(inside.mean() * (1-inside.mean()) / n)
    return pi_est, std_err

print('Monte Carlo estimation of pi:')
print(f'  {'n':>9}  {'pi estimate':>12}  {'std error':>10}  {'true error':>12}')
for n in [100, 1000, 10000, 100000, 1000000]:
    pi_est, se = estimate_pi(n)
    print(f'  {n:>9}  {pi_est:>12.6f}  {se:>10.6f}  {abs(pi_est-np.pi):>12.6f}')

print(f'\nTrue pi = {np.pi:.6f}')
print(f'Error ~ O(1/sqrt(n)): halving error needs 4x more samples')

Code cell 16

# === 4.2 Variance Reduction: Antithetic Variates ===

np.random.seed(42)

def mc_plain(f, n):
    u = np.random.rand(n)
    return f(u).mean(), f(u).std() / np.sqrt(n)

def mc_antithetic(f, n):
    u = np.random.rand(n // 2)
    vals = (f(u) + f(1 - u)) / 2
    return vals.mean(), vals.std() / np.sqrt(n // 2)

# f(u) = sqrt(u), exact = 2/3
f = np.sqrt
exact = 2/3

print('Antithetic variates for integral of sqrt(u) on [0,1]:')
print(f'  {'n':>7}  {'Plain SE':>10}  {'Antithetic SE':>14}  {'Speedup':>8}')
for n in [100, 1000, 10000]:
    _, se_p = mc_plain(f, n)
    _, se_a = mc_antithetic(f, n)
    print(f'  {n:>7}  {se_p:>10.6f}  {se_a:>14.6f}  {se_p/se_a:>8.2f}x')

Code cell 17

# === 4.3 Monte Carlo vs Quasi-MC (Sobol) ===

try:
    from scipy.stats import qmc
    HAS_QMC = True
except ImportError:
    HAS_QMC = False

# Integrate f(x) = sin(pi*x1)*sin(pi*x2) on [0,1]^2
# Exact: (2/pi)^2 = 0.4053...
f_2d = lambda X: np.prod(np.sin(np.pi * X), axis=1)
exact_2d = (2/np.pi)**2

print(f'Integrating prod(sin(pi*xi)) on [0,1]^d:')
print(f'Exact value: {exact_2d:.6f}')

np.random.seed(42)
errors_mc, errors_qmc = [], []
n_vals = [64, 128, 256, 512, 1024, 2048]

for n in n_vals:
    # Plain MC
    X_mc = np.random.rand(n, 2)
    est_mc = f_2d(X_mc).mean()
    errors_mc.append(abs(est_mc - exact_2d))

    # Quasi-MC (Sobol)
    if HAS_QMC:
        sampler = qmc.Sobol(d=2, scramble=True, seed=42)
        X_qmc = sampler.random(n)
        est_qmc = f_2d(X_qmc).mean()
        errors_qmc.append(abs(est_qmc - exact_2d))
    else:
        errors_qmc.append(float('nan'))

print(f'  {'n':>6}  {'MC error':>12}  {'QMC error':>12}')
for n, em, eq in zip(n_vals, errors_mc, errors_qmc):
    print(f'  {n:>6}  {em:>12.4e}  {eq:>12.4e}')

if not HAS_QMC:
    print('Note: scipy.stats.qmc not available; install scipy >= 1.7 for QMC support')

5. Gaussian Quadrature for ML: Gauss-Hermite and ELBO

Code cell 19

# === 5.1 Gauss-Hermite Quadrature ===

def gauss_hermite_expect(f, mu, sigma, n=20):
    """
    Compute E_{x~N(mu,sigma^2)}[f(x)] via n-point GH quadrature.
    Uses the transform x = sqrt(2)*sigma*t + mu.
    """
    t, w = np.polynomial.hermite.hermgauss(n)
    x_mapped = np.sqrt(2) * sigma * t + mu
    return np.dot(w, f(x_mapped)) / np.sqrt(np.pi)

# Test: E[x^2] = mu^2 + sigma^2 (known analytically)
mu, sigma = 2.0, 0.5
f_sq = lambda x: x**2
exact_e2 = mu**2 + sigma**2  # = 4.25

for n in [2, 5, 10, 20]:
    gh_est = gauss_hermite_expect(f_sq, mu, sigma, n)
    print(f'GH-{n:>2}: E[x^2] = {gh_est:.10f}, error = {abs(gh_est-exact_e2):.2e}')

print(f'Exact E[x^2] = {exact_e2:.10f}')

# ML application: E[log(1+exp(x))] for x ~ N(0,1) (softplus expectation)
f_softplus = lambda x: np.log(1 + np.exp(x))
exact_sp, _ = integrate.quad(lambda x: f_softplus(x) * np.exp(-x**2/2) / np.sqrt(2*np.pi),
                              -10, 10)

print(f'\nE[softplus(x)] for x~N(0,1):')
print(f'  Exact (quad): {exact_sp:.10f}')
for n in [3, 5, 10, 20]:
    est = gauss_hermite_expect(f_softplus, 0, 1, n)
    print(f'  GH-{n:>2}:       {est:.10f}, error = {abs(est-exact_sp):.2e}')

# Monte Carlo comparison
mc_est = f_softplus(np.random.randn(100000)).mean()
print(f'  MC-1e5:       {mc_est:.6f}, error = {abs(mc_est-exact_sp):.4f}')

Code cell 20

# === 5.2 Reparameterization Trick vs Score Function ===

# Compute d/dmu E_{x~N(mu,1)}[-x^2/2] = E[d/dmu (-x^2/2)]
# = E[-x] = -mu  (via reparameterization: x = mu + eps, eps~N(0,1))
# Also = -mu via score function: E[(-x^2/2) * d/dmu log N(x;mu,1)]
#                              = E[(-x^2/2) * (x-mu)]

np.random.seed(42)
mu_val = 1.5
n_samples = 10000

# Exact gradient
exact_grad = -mu_val  # d/dmu E[-x^2/2] = -mu

# Reparameterization: x = mu + eps, eps~N(0,1), d/dmu f(x) = d/dmu(-x^2/2) * dx/dmu = -x * 1
eps = np.random.randn(n_samples)
x_reparam = mu_val + eps
grad_reparam = -x_reparam  # d/dmu (-x^2/2) = -x, dx/dmu = 1
mean_reparam = grad_reparam.mean()
std_reparam = grad_reparam.std() / np.sqrt(n_samples)

# Score function: E[f(x) * d/dmu log p(x;mu)] = E[(-x^2/2) * (x-mu)]
x_score = np.random.normal(mu_val, 1, n_samples)
f_score = -x_score**2/2
score = (x_score - mu_val)  # d/dmu log N(x;mu,1)
grad_score = (f_score * score).mean()
std_score = (f_score * score).std() / np.sqrt(n_samples)

print(f'Exact gradient: {exact_grad}')
print(f'Reparameterization: {mean_reparam:.4f} ± {std_reparam:.4f} (SE)')
print(f'Score function:     {grad_score:.4f} ± {std_score:.4f} (SE)')
print(f'\nVariance ratio (score/reparam): {std_score**2/std_reparam**2:.2f}x higher variance')

6. Summary Verification Suite

Code cell 22

# === 6. Summary Verification ===

print('=' * 60)
print('SUMMARY: Numerical Integration Verification')
print('=' * 60)

results = []
exact = np.e - 1

# Test 1: Trapezoidal O(h^2)
e1 = abs(composite_trap(np.exp, 0, 1, 100) - exact)
e2 = abs(composite_trap(np.exp, 0, 1, 200) - exact)
ok1 = abs(e1/e2 - 4.0) < 0.5
results.append(('Trapezoidal O(h^2) convergence', ok1))

# Test 2: Simpson's O(h^4)
e1 = abs(composite_simpson(np.exp, 0, 1, 10) - exact)
e2 = abs(composite_simpson(np.exp, 0, 1, 20) - exact)
ok2 = abs(e1/e2 - 16.0) < 2.0
results.append(("Simpson's O(h^4) convergence", ok2))

# Test 3: GL-n degree of exactness
n_gl = 5
ok3 = abs(gl_integrate(lambda x: x**(2*n_gl-1), -1, 1, n_gl) - 0) < 1e-10
results.append((f'GL-{n_gl} exact for degree {2*n_gl-1}', ok3))

# Test 4: Romberg super-convergence
R_tbl, _ = romberg(np.exp, 0, 1, max_k=6)
ok4 = abs(R_tbl[-1,-1] - exact) < 1e-14
results.append(('Romberg near machine precision', ok4))

# Test 5: Gauss-Hermite E[x^2] = mu^2 + sigma^2
gh_test = gauss_hermite_expect(lambda x: x**2, mu=2.0, sigma=0.5, n=5)
ok5 = abs(gh_test - (4.0 + 0.25)) < 1e-12
results.append(('Gauss-Hermite E[x^2] exact', ok5))

# Test 6: GL weights sum to 2
_, w_test = gauss_legendre(7)
ok6 = abs(w_test.sum() - 2.0) < 1e-12
results.append(('GL-7 weights sum to 2', ok6))

print()
all_pass = True
for name, ok in results:
    print(f'  {'PASS' if ok else 'FAIL'}  {name}')
    if not ok: all_pass = False

print()
print('=' * 60)
print(f'  {'ALL TESTS PASS' if all_pass else 'SOME TESTS FAILED'}')
print('=' * 60)