Simpson's Rule

Simpson's rule fits parabolas through triples of points, achieving fourth-order accuracy that integrates cubics exactly.

Parabolas instead of lines

Simpson's rule improves on the trapezoidal rule by fitting a parabola through three equally spaced points and integrating it exactly. On a panel of width 2h the rule is (h/3)(f_0 + 4 f_1 + f_2). The composite version pairs up panels across the interval, requiring an even number of subintervals.

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

Fourth-order accuracy

The composite Simpson's rule has error proportional to h^4, so halving the step reduces error roughly sixteenfold. Although built from parabolas, it integrates cubic polynomials exactly, an extra degree of exactness that comes free from symmetry. Its error coefficient involves the fourth derivative of the integrand.

Variants

Simpson's 3/8 rule uses cubics through four points and handles panel counts not divisible by two in the 1/3 scheme. Boole's rule and higher Newton-Cotes rules raise the order further, but very high-order Newton-Cotes rules develop negative weights and become unstable, so practice rarely goes beyond Simpson unless using Gaussian quadrature.

When it shines and when it does not

Simpson's rule is an excellent general-purpose choice for smooth integrands, offering high accuracy for little code. It struggles with singularities, discontinuities, and sharp peaks, where adaptive subdivision or specialized rules are needed. For those cases adaptive Simpson recursively refines only the difficult subintervals.

Simpson-type quadrature appears throughout scientific computing, from energy-integrated cross sections to spatial moments in the breeder Hyperion models, wherever a smooth one-dimensional integrand needs efficient evaluation.