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

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

where $[a,b]\,\!$ is an interval with midpoint $m\,\!$, $S(a,b)\,\!$, $S(a,m)\,\!$, and $S(m,b)\,\!$ are the estimates given by Simpson's rule on the corresponding intervals and $\epsilon \,\!$ 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 $S(a,m)+S(m,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,m)+S(m,b)-S(m,b)]/15\,$. The thus obtained estimate is exact for polynomials of degree five or less.

Sample implementations

A common implementation technique shown below is passing down f(a), f(b), f(m) along with the interval [a, b]. These values, used for evaluating S(a, b) 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
# "structured" adaptive version, translated from Racket
def simpsons_1call(f, a, fa, b, fb):
    """Evaluates the Simpson's Rule, also returning m and f(m) to reuse"""
    m = (a+b) / 2
    fm = f(m)
    return (m, fm, abs(b-a) / 6 * (fa + 4 * fm + fb))

def recursive_asr(f, a, fa, b, fb, eps, whole, m, fm):
    """
    Efficient recursive implementation of adaptive Simpson's rule.
    Function values at the start, middle, end of the intervals are retained.
    """
    lm, flm, left  = simpsons_1call(f, a, fa, m, fm)
    rm, frm, right = simpsons_1call(f, m, fm, b, fb)
    delta = left + right - whole
    if abs(delta) <= 15 * eps:
        return left + right + delta / 15
    return recursive_asr(f, a, fa, m, fm, eps/2, left , lm, flm) +\
           recursive_asr(f, m, fm, b, fb, eps/2, right, rm, frm)

def adaptive_simpsons_rule(f, a, b, eps):
    """Calculate integral of f from a to b with max error of eps."""
    fa, fb = f(a), f(b)
    m, fm, whole = simpsons_1call(f, a, fa, b, fb)
    return recursive_asr(f, a, fa, b, fb, eps, whole, m, fm)

from math import sin
print(adaptive_simpsons_rule(sin, 0, 1, 1e-09))

C

Here is an implementation of the adaptive Simpson's method in C99 that avoids redundant evaluations of f and quadrature computations.

#include <math.h>  // include file for fabs and sin
#include <stdio.h> // include file for printf and perror
#include <errno.h>

/** Adaptive Simpson's Rule Wrapper
 *  (fills in cached function evaluations) */
float adaptiveSimpsons(float (*f)(float),     // function ptr to integrate
                       float a, float b,      // interval [a,b]
                       float epsilon,         // error tolerance
                       int maxRecDepth) {     // recursion cap
    errno = 0;
    float h = b - a;
    float fa = f(a), fb = f(b), fm = f((a + b)/2);
    float S = (h/6)*(fa + 4*fm + fb);
    return adaptiveSimpsonsAux(f, a, b, epsilon, S, fa, fb, fm, maxRecDepth);
}

/** Adaptive Simpson's Rule, Recursive Core */
float adaptiveSimpsonsAux(float (*f)(float), float a, float b, float epsilon,
                          float whole, float fa, float fb, float fm, int rec) {
    float m   = (a + b)/2,  h   = (b - a)/2;
    float lm  = (a + m)/2,  rm  = (m + b)/2;
    float flm = f(lm),      frm = f(rm);
    float left  = (h/6) * (fa + 4*flm + fm);
    float right = (h/6) * (fm + 4*frm + fb);
    float delta = left + right - whole;

    if (rec <= 0) errno = ERANGE; // Leave a note if we bailed out early
    // Lyness 1969 + Richardson extrapolation; see article
    if (rec <= 0 || fabs(delta) <= 15*epsilon)
        return left + right + (delta)/15;
    return adaptiveSimpsonsAux(f, a, m, epsilon/2, left,  fa, fm, flm, rec-1) +
           adaptiveSimpsonsAux(f, m, b, epsilon/2, right, fm, fb, frm, rec-1);
}

/** usage example */
int main() {
    // Let I be the integral of sin(x) from 0 to 2
    float I = adaptiveSimpsons(sin, 0, 2, 0.001, 100);
    printf("integrate(sin, 0, 2) = %lf\n", I);   // print the result
    perror("adaptiveSimpsons");                  // Was it successful?
    return 0;
}

This implementation has been incorprated 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 
(provide/contract
 [adaptive-simpson 
  (->i ((f (-> real? real?)) (L real?) (R  (L) (and/c real? (>/c L))))
       (#:epsilon (ε real?))
       (r real?))])


;; -----------------------------------------------------------------------------
;; implementation 
(define (adaptive-simpson f L R #:epsilon [ε .000000001])
  (define f@L (f L))
  (define f@R (f R))
  (define-values (M f@M whole) (simpson-1call-to-f f L f@L R f@R))
  (asr f L f@L R f@R ε whole M f@M))

;; the "efficient" implementation
(define (asr f L f@L R f@R ε whole M f@M)
  (define-values (leftM  f@leftM  left*)  (simpson-1call-to-f f L f@L M f@M))
  (define-values (rightM f@rightM right*) (simpson-1call-to-f f M f@M R f@R))
  (define delta* (- (+ left* right*) whole))
  (cond
    [(<= (abs delta*) (* 15 ε)) (+ left* right* (/ delta* 15))]
    [else (define epsilon1 (/ ε 2))
          (+ (asr f L f@L M f@M epsilon1 left*  leftM  f@leftM) 
             (asr f M f@M R f@R epsilon1 right* rightM f@rightM))]))

;; evaluate half an intervel (1 func eval)
(define (simpson-1call-to-f f L f@L R f@R)
  (define M (mid L R))
  (define f@M (f M))
  (values M f@M (* (/ (abs (- R L)) 6) (+ f@L (* 4 f@M) f@R))))

(define (mid L R) (/ (+ L R) 2.))

The code is an excerpt of a "#lang racket" module and that includes a (require rackunit) line.

Related algorithms

• 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.
• Kuncir's Algorithm 103 (1962) is the original recursive, bisecting, adaptive integrator. Algorithm 103 consists of a larger routine with a nested subroutine (loop AA), made recursive by the use of the goto statement. It guards against the underflowing of interval widths (loop BB), and aborts as soon as the user-specified eps is exceeded. The termination criteria is $|S^{(2)}(a,b)-S(a,b)|<2^{-n}\epsilon \,$, where n is the current level of recursion and S⁽²⁾ is the more accurate estimate.
• McKeeman's Algorithm 145 (1962) is a similarly recursive integrator that splits the interval into three instead of two parts. The recursion is written in a more familiar manner. The 1962 algorithm, found to be over-cautious, uses $|S^{(3)}(a,b)-S(a,b)|<3^{-n}\epsilon \,$ for termination, so a 1963 improvement uses ${\sqrt {3}}^{-n}\epsilon$ 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). Modification 4, not described here, contains provisions for roundoff error that override the user-specified precision to match the best achievable. The new return value with extrapolation is related to Boole's rule and Romberg's method.

External links

• Module for Adaptive Simpson's Rule