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 • "structured" adaptive version, translated from Racket def _quad_simpsons_mem(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 _quad_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 = _quad_simpsons_mem(f, a, fa, m, fm) rm, frm, right = _quad_simpsons_mem(f, m, fm, b, fb) delta = left + right - whole if abs(delta)
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. • include // include file for fabs and sin • include // include file for printf and perror • include /** Adaptive Simpson's Rule, Recursive Core */ float adaptiveSimpsonsAux(float (*f)(float), float a, float b, float eps, 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; // serious numerical trouble: it won't converge if ((eps/2 == eps) || (a == lm)) { errno = EDOM; return whole; } 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 // for the hostile example (rand function) static int callcnt = 0; static float sinfc(float x) { callcnt++; return sinf(x); } static float frand48c(float x) { callcnt++; return drand48(); } int main() { // Let I be the integral of sin(x) from 0 to 2 float I = adaptiveSimpsons(sinfc, 0, 2, 1e-5, 3); printf("integrate(sinf, 0, 2) = %lf\n", I); // print the result perror("adaptiveSimpsons"); // Was it successful? (depth=1 is too shallow) printf("(%d evaluations)\n", callcnt); callcnt = 0; srand48(0); I = adaptiveSimpsons(frand48c, 0, 0.25, 1e-5, 25); // a random function printf("integrate(frand48, 0, 0.25) = %lf\n", I); perror("adaptiveSimpsons"); // won't converge printf("(%d evaluations)\n", callcnt); return 0; } 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. ;; ----------------------------------------------------------------------------- ;; 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 [( ==Related algorithms==