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.
A criterion for determining when to stop subdividing an interval, suggested by J.N. Lyness, is
where is an interval with midpoint,,, and are the estimates given by Simpson's rule on the corresponding intervals and is the desired tolerance for the interval.
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. The thus obtained estimate is exact for polynomials of degree five or less.

Numerical consideration

Some inputs will fail to converge in adaptive Simpson's method quickly, resulting in the tolerance underflowing and producing an infinite loop. Simple methods of guarding against this problem include adding a depth limitation, verifying that in floating-point arithmetics, or both. The interval size may also approach the local machine epsilon, giving.
Lychee's 1969 paper includes a "Modification 4" that addresses this problem in a more concrete way:
The epsilon-raising maneuver allows the routine to be used in a "best effort" mode: given a zero initial tolerance, the routine will try to get the most precise answer and return an actual error level.

Sample implementations

A common implementation technique shown below is passing down along with the interval. These values, used for evaluating at the parent level, will again be used for the subintervals. Doing so cuts down the cost of each recursive call from 6 to 2 evaluations of the input function. The size of the stack space used stays linear to the layer of recursions.

Python

Here is an implementation of adaptive Simpson's method in Python.

from __future__ import division # python 2 compat
  1. "structured" adaptive version, translated from Racket
def _quad_simpsons_mem:
"""Evaluates the Simpson's Rule, also returning m and f to reuse"""
m = / 2
fm = f
return / 6 * )
def _quad_asr:
"""
Efficient recursive implementation of adaptive Simpson's rule.
Function values at the start, middle, end of the intervals are retained.
"""
lm, flm, left = _quad_simpsons_mem
rm, frm, right = _quad_simpsons_mem
delta = left + right - whole
if abs <= 15 * eps:
return left + right + delta / 15
return _quad_asr +\
_quad_asr
def quad_asr:
"""Integrate f from a to b using Adaptive Simpson's Rule with max error of eps."""
fa, fb = f, f
m, fm, whole = _quad_simpsons_mem
return _quad_asr
from math import sin
print

C

Here is an implementation of the adaptive Simpson's method in C99 that avoids redundant evaluations of f and quadrature computations. It includes all three "simple" defenses against numerical problems.

  1. include // include file for fabs and sin
  2. include // include file for printf and perror
  3. include
/** Adaptive Simpson's Rule, Recursive Core */
float adaptiveSimpsonsAux, float a, float b, float eps,
float whole, float fa, float fb, float fm, int rec)
/** Adaptive Simpson's Rule Wrapper
* */
float adaptiveSimpsons, // function ptr to integrate
float a, float b, // interval
float epsilon, // error tolerance
int maxRecDepth)
/** usage example */
  1. include // for the hostile example
static int callcnt = 0;
static float sinfc
static float frand48c
int main

This implementation has been incorporated into a C++ ray tracer intended for X-Ray Laser simulation at Oak Ridge National Laboratory, among other projects. The ORNL version has been enhanced with a call counter, templates for different datatypes, and wrappers for integrating over multiple dimensions.

Racket

Here is an implementation of the adaptive Simpson method in Racket with a behavioral software contract. The exported function computes the indeterminate integral for some given function f.

;; -----------------------------------------------------------------------------
;; interface, with contract
)

)])
;; -----------------------------------------------------------------------------
;; implementation


)
)
;; the "efficient" implementation
)
)

) ]
))
;; evaluate half an intervel


))

Related algorithms