Adaptive Simpson's method

Adaptive Simpson's Method, also called Adaptive Simpson's Rule, is a method of numerical integration proposed by W.H.McKeeman in 1962.. It is perhaps the first adaptive algorithm for numerical integration to appear in print. Adaptive Simpson's Method uses a heuristic to estimate the error we get from calculating a definite integral using Simpson's Rule. If the error exceeds a user-specified tolerance, the algorithm calls for subdividing the interval of integration in two and applying Adaptive Simpson's Method to each subinterval in a recursive manner. The technique is theoretically more efficient than Composite simpson's rule since it uses fewer function evaluations in places where the function is well-approximated by a quadratic function.

The heuristic for determining when to stop subdividing an interval is

$|S(a,c)+S(c,b)-S(a,b)|/15<\epsilon$

where $[a,b]$ is an interval with midpoint $c$, $S(a,b)$ is the estimate given by Simpson's Rule on the interval $[a,b]$ and $\epsilon$ is the desired tolerance for the interval.

Here is an implementation of Adaptive Simpson's Rule in Python.

def simpsons_rule(f,a,b):
    c = (a+b) / 2.0
    h3 = abs(b-a) / 6.0
    return h3 * (f(a) + 4.0*f(c) + f(b))

def recursive_asr(f,a,b,eps,sum,n=0):
    "Recursive implementation of adaptive Simpson's rule."
    print "a: ", a, " b: ", b, " n: ", n, "\n"
    c = (a+b) / 2.0
    left = simpsons_rule(f,a,c)
    right = simpsons_rule(f,c,b)
    if ( abs(left + right - sum) <= 15 * eps ):
        return (left + right +(left + right - sum)/15)
    return (recursive_asr(f,a,c,eps/2,left,n+1) + recursive_asr(f,c,b,eps/2,right,n+1))

def adaptive_simpsons_rule(f,a,b,eps):
    "Calculate integral of f from a to b with max error of eps."
    return (recursive_asr(f,a,b,eps,simpsons_rule(f,a,b)))

from math import sin
print adaptive_simpsons_rule(sin,0,1,.000000001)