Numerical Integration with Simpson's Rule

Approximate a definite integral by fitting parabolas over pairs of intervals, achieving fourth-order accuracy for smooth functions.

The idea

The trapezoidal rule fits straight lines between sample points. Simpson's rule fits a parabola through every three consecutive points, integrating that exactly. For smooth integrands it is dramatically more accurate for the same number of evaluations.

The formula

With an even number n of intervals of width h = (b-a)/n, the integral is approximately (h/3)[f0 + 4(f1+f3+...) + 2(f2+f4+...) + fn]: endpoints weight 1, odd interior points weight 4, even interior points weight 2.

import numpy as np
def simpson(f,a,b,n):
    if n%2: n+=1
    x=np.linspace(a,b,n+1); h=(b-a)/n; y=f(x)
    return h/3*(y[0]+y[-1]+4*y[1:-1:2].sum()+2*y[2:-1:2].sum())
print(round(simpson(np.sin,0,np.pi,10),8))  # exact = 2
print(round(simpson(lambda x:x**4,0,1,4),8)) # exact = 0.2

Error behaviour

Simpson's error scales as h^4, so doubling the sample count cuts error by sixteen. Notably it integrates cubics exactly - the parabola fit gets one bonus order for free by symmetry. Only from the fourth derivative onward does error appear.

Caveats

The accuracy assumes a smooth integrand. For functions with kinks, singularities, or sharp peaks, blindly refining Simpson's rule converges slowly; adaptive quadrature that subdivides where the integrand is difficult is the practical answer. Simpson's rule remains the go-to for smooth, well-behaved integrals.