Adaptive Simpson's method

Adaptive Simpson's method, also called adaptive Simpson's rule, is a method of numerical integration proposed by G. F. Kuncir in 1962. It is probably the first recursive 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 cubic function. Simpson's rule is an interpolatory quadrature rule which is exact when the integrand is a polynomial of degree three or lower. Using Richardson extrapolation, the more accurate Simpson estimate for six function values is combined with the less accurate estimate for three function values by applying the correction So, the obtained estimate is exact for polynomials of degree five or less.

A criterion for determining when to stop subdividing an interval, suggested by J. N. Lyness, is

where is an interval with midpoint , while , , and given by Simpson's rule are the estimates of , , and respectively, and is the desired maximum error tolerance for the interval. Note, .

To perform adaptive Simpson's method, do the following: if , add and to the sum of Simpson's rules which are used to approximate the integral, otherwise, perform the same operation with and instead of .

Henriksson (1961) is a non recursive variant of Simpson's Rule. It "adapts" by integrating from left to right and adjusting the interval width as needed. The 1962 algorithm, found to be over cautious, uses for termination, so a 1963 improvement uses instead. Lyness (1969) is almost the current integrator. Created as a set of four modifications of McKeeman 1962, it replaces trisection with bisection to lower computational costs (Modifications 1+2, coinciding with the Kuncir integrator) and improves McKeeman's 1962/63 error estimates to the fifth order (Modification 3), in a way related to Boole's rule and Romberg's method. Modification 4, not implemented here, contains provisions for roundoff error that allows for raising the ε to the minimum allowed by current precision and returning the new error.