Adaptive Simpson's method, also called adaptive Simpson's rule, is a method of numerical integration

Numerical integration

In numerical analysis, numerical integration constitutes a broad family of algorithms for calculating the numerical value of a definite integral, and by extension, the term is also sometimes used to describe the numerical solution of differential equations. This article focuses on calculation of...

 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

Simpson's rule

In numerical analysis, Simpson's rule is a method for numerical integration, the numerical approximation of definite integrals. Specifically, it is the following approximation:...

. 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

Quadratic function

A quadratic function, in mathematics, is a polynomial function of the formf=ax^2+bx+c,\quad a \ne 0.The graph of a quadratic function is a parabola whose axis of symmetry is parallel to the y-axis....

.

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

Simpson's rule

In numerical analysis, Simpson's rule is a method for numerical integration, the numerical approximation of definite integrals. Specifically, it is the following approximation:...

 is an interpolatory quadrature rule which is exact when the integrand is a polynomial of degree three or lower. Using Richardson extrapolation

Richardson extrapolation

In numerical analysis, Richardson extrapolation is a sequence acceleration method, used to improve the rate of convergence of a sequence. It is named after Lewis Fry Richardson, who introduced the technique in the early 20th century. In the words of Birkhoff and Rota, ".....

, 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.

Python

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

Python (programming language)

Python is a general-purpose, high-level programming language whose design philosophy emphasizes code readability. Python claims to "[combine] remarkable power with very clear syntax", and its standard library is large and comprehensive...

. 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,whole):
"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 - whole) <= 15*eps:
return left + right + (left + right - whole)/15.0
return recursive_asr(f,a,c,eps/2.0,left) + recursive_asr(f,c,b,eps/2.0,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 sin
print adaptive_simpsons_rule(sin,0,1,.000000001)

C

Here is an implementation of the adaptive Simpson's method in C99 that avoids redundant evaluations of f and quadrature computations.
The amount of memory used is O(D) where D is the maximum recursion depth. Each stack frame caches computed values
that may be needed in subsequent calls.

1. include // include file for fabs and sin
2. include // include file for printf

//
// Recursive auxiliary function for adaptiveSimpsons function below
//
double adaptiveSimpsonsAux(double (*f)(double), double a, double b, double epsilon,
double S, double fa, double fb, double fc, int bottom) {
double c = (a + b)/2, h = b - a;
double d = (a + c)/2, e = (c + b)/2;
double fd = f(d), fe = f(e);
double Sleft = (h/12)*(fa + 4*fd + fc);
double Sright = (h/12)*(fc + 4*fe + fb);
double S2 = Sleft + Sright;
if (bottom <= 0 || fabs(S2 - S) <= 15*epsilon)
return S2 + (S2 - S)/15;
return adaptiveSimpsonsAux(f, a, c, epsilon/2, Sleft, fa, fc, fd, bottom-1) +
adaptiveSimpsonsAux(f, c, b, epsilon/2, Sright, fc, fb, fe, bottom-1);
}

//
// Adaptive Simpson's Rule
//
double adaptiveSimpsons(double (*f)(double), // ptr to function
double a, double b, // interval [a,b]
double epsilon, // error tolerance
int maxRecursionDepth) { // recursion cap
double c = (a + b)/2, h = b - a;
double fa = f(a), fb = f(b), fc = f(c);
double S = (h/6)*(fa + 4*fc + fb);
return adaptiveSimpsonsAux(f, a, b, epsilon, S, fa, fb, fc, maxRecursionDepth);
}
int main{
double I = adaptiveSimpsons(sin, 0, 1, 0.000000001, 10); // compute integral of sin(x)
// from 0 to 1 and store it in
// the new variable I
printf("I = %lf\n",I); // print the result
return 0;
}

External links

• Module for Adaptive Simpson's Rule