Adaptive Simpson's method

Adaptive Simpson's method, also called adaptive Simpson's rule, is a method of numerical integration proposed by William M. McKeeman in 1962.William M. McKeeman: Algorithm 145: Adaptive numerical integration by Simpson's rule. Commun. ACM 5(12): 604 (1962). doi|10.1145/355580.369102] It is perhaps the first adaptive algorithm for numerical integration to appear in print, although more modern adaptive methods based on Gauss–Kronrod quadrature and Clenshaw–Curtis quadrature are now generally preferred. Adaptive Simpson's method uses an estimate of 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 usually much 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 criterion for determining when to stop subdividing an interval is

:

where $[a,b]\; ,!$ is an interval with midpoint $c,!$, $S(a,b),!$, $S(a,c),!$, and $S(c,b),!$ are the estimates given by Simpson's rule on the corresponding intervals and $epsilon,!$ is the desired tolerance for the interval.Numerical Analysis. Mathematics of Scientific Computing Third Edition. by David Kincaid and Ward Cheney. Published by Brooks/Cole Third Edition, 2002]

Simpson's rule is a special case of Romberg's method. The theory of this method shows that Simpson's rule is exact when the integrand is a polynomial of degree three or lower. According to Romberg's method, the more accurate Simpson estimate $S(a,c)\; +\; S(c,b),$ for six function values is combined with the less accurate estimate $S(a,b),$ for three function values by applying the correction $[S(a,c)\; +\; S(c,b)\; -\; S(a,b)]\; /15\; ,$. Here the constant 15 is chosen to ensure that after applying the correction, an estimate is obtained that is exact for polynomials of degree five or less.

Here is an implementation of adaptive Simpson's method in Python. Note that this is explanatory code, without regard for efficiency. Every call to recursive_asr entails six function evaluations. For actual use, one will want to modify it so that the minimum of two function evaluations are performed.

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): "Recursive implementation of adaptive Simpson's rule." 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) + recursive_asr(f,c,b,eps/2,right)

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 sinprint adaptive_simpsons_rule(sin,0,1,.000000001)

Bibliography

External links

* [http://math.fullerton.edu/mathews/n2003/AdaptiveQuadMod.html Module for Adaptive Simpson's Rule]