Numerical Integration

An integral is an area under a curve — exact formulas exist for some functions, but many real functions (or functions you only have as code, not a formula) need a numerical approximation. The trapezoidal rule approximates the area under the curve as a series of trapezoids instead of the true curved shape:

def trapezoidal_rule(f, a, b, n):
    h = (b - a) / n
    total = (f(a) + f(b)) / 2
    for i in range(1, n):
        total += f(a + i * h)
    return total * h

print(trapezoidal_rule(lambda x: x**2, 0, 1, 1000))   # very close to 1/3

Each trapezoid's area is h * (left height + right height) / 2 — summed across all n of them, with each interior point's height counted once (it's the right edge of one trapezoid and the left edge of the next).

Simpson's rule fits a parabola through each group of three points instead of a straight line — parabolas hug curves far more closely than straight edges do, so Simpson's rule reaches the same accuracy with far fewer intervals:

def simpsons_rule(f, a, b, n):   # n must be even
    h = (b - a) / n
    total = f(a) + f(b)
    for i in range(1, n):
        weight = 4 if i % 2 != 0 else 2
        total += weight * f(a + i * h)
    return total * h / 3

The alternating 4/2/4/2 weighting comes directly from the algebra of fitting parabolas through consecutive triples of points — it looks arbitrary but isn't. Simpson's rule is exact (zero error) for any polynomial up to degree 3, which is why it converges so much faster than trapezoidal on smooth, curved functions.