Numerical Integration — Simpson's Rule and Trapezoidal Method
- Trapezoidal rule: O(h²) error — doubling n reduces error by 4x.
- Simpson's rule: O(h⁴) error — doubling n reduces error by 16x. Prefer over Trapezoidal.
- Simpson's with n=10 beats Trapezoidal with n=10,000 for smooth functions.
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.
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}')
n= 100: 1.9998355038, error=1.64e-04
n= 1000: 1.9999983550, error=1.64e-06
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²).
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}')
n= 100: 2.0000000000, error=1.06e-15
n= 1000: 2.0000000000, error=0.00e+00
Simpson vs Trapezoid comparison (n=10):
Trapezoid error: 1.65e-02
Simpson error: 1.07e-07
Error Analysis
| Method | Error Order | n=10 error (sin) | n=100 error | |---|---|---|---| | Trapezoidal | O(h²) | 1.65×10⁻² | 1.64×10⁻⁴ | | Simpson's | O(h⁴) | 1.07×10⁻⁷ | 1.06×10⁻¹⁵ | | Gaussian (n=5) | Exact for deg≤9 | ~0 | ~0 |
Simpsons's with n=10 is more accurate than Trapezoidal with n=10,000 — same number of function evaluations, 1000x better.
Adaptive Integration
Fixed n doesn't work for all functions — some regions need more resolution. Adaptive integration refines intervals where the error is large.
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}')
scipy.quad = 0.746824132812
🎯 Key Takeaways
- Trapezoidal rule: O(h²) error — doubling n reduces error by 4x.
- Simpson's rule: O(h⁴) error — doubling n reduces error by 16x. Prefer over Trapezoidal.
- Simpson's with n=10 beats Trapezoidal with n=10,000 for smooth functions.
- Adaptive integration automatically refines intervals where error is high.
- Use scipy.integrate.quad in production — it implements adaptive Gaussian quadrature.
Interview Questions on This Topic
- QWhy is Simpson's rule more accurate than the Trapezoidal rule for the same number of intervals?
- QHow many intervals does the Trapezoidal rule need to achieve 10^-6 error vs Simpson's rule?
- QWhat is adaptive integration and why is it useful?
- QDerive the Trapezoidal rule formula geometrically.
Frequently Asked Questions
When should I use Gaussian quadrature instead of Simpson's rule?
Gaussian quadrature (n points integrate polynomials of degree 2n-1 exactly) is superior when the function is smooth and you want maximum accuracy with minimum evaluations. Simpson's rule is simpler and sufficient for most engineering applications. scipy.integrate.quad uses adaptive Gauss-Kronrod quadrature.
Developer and founder of TheCodeForge. I built this site because I was tired of tutorials that explain what to type without explaining why it works. Every article here is written to make concepts actually click.