Numerical Integration
February 9, 2010
The three simple methods are similar; w is the width of a slice, g is the approximation function, and lo, mid and hi are the x-coordinates at three points in the slice:
(define (int-rect f a b n)
(let* ((w (/ (- b a) n))
(mid (lambda (i) (+ a (/ w -2) (* i w))))
(g (lambda (i) (* w (f (mid i))))))
(sum-of (g i) (i range 1 (+ n 1)))))
(define (int-trap f a b n)
(let* ((w (/ (- b a) n))
(lo (lambda (i) (+ a (* (- i 1) w))))
(hi (lambda (i) (+ a (* i w))))
(g (lambda (i) (* w 1/2 (+ (f (lo i)) (f (hi i)))))))
(sum-of (g i) (i range 1 (+ n 1)))))
(define (int-simp f a b n)
(let* ((w (/ (- b a) n))
(lo (lambda (i) (+ a (* (- i 1) w))))
(mid (lambda (i) (+ a (/ w -2) (* i w))))
(hi (lambda (i) (+ a (* i w))))
(g (lambda (i) (* w 1/6 (+ (f (lo i)) (* 4 (f (mid i))) (f (hi i)))))))
(sum-of (g i) (i range 1 (+ n 1)))))
Adaptive quadrature integrates over five and ten slices and measures the difference between the two approximations. Recursion stops when the difference is within the desired tolerance; otherwise, the slice is split in two and the two sub-slices are integrated separately. Integrate allows you to select the simple method you wish:
(define (integrate method f a b . epsilon)
(let ((g5 (method f a b 5))
(g10 (method f a b 10))
(mid (/ (+ a b) 2))
(eps (if (null? epsilon) #e1e-7 (car epsilon))))
(if (< (abs (- g10 g5)) eps) g10
(+ (integrate method f a mid eps)
(integrate method f mid b eps)))))
The logarithmic integral is a well-known approximation to the prime-counting function pi:
(define (approx-pi n)
(inexact->exact
(round
(integrate int-simp
(lambda (x) (/ (log x)))
2 n 1e-7))))
We give examples below, including the definite integral of the cube function from zero to one, which is 1/4 by analysis:
> (define (cube x) (* x x x))
> (int-rect cube 0 1 10000.)
0.24999999875000073
> (int-trap cube 0 1 10000.)
0.2500000025
> (int-simp cube 0 1 10000.)
0.25
> (approx-pi 1e21)
21127269486616129536
The correct value of the logarithmic integral for 1021 is 21127269486616126181.3; we lost some accuracy after fifteen digits. The correct value of pi is 21127269486018731928, which is astonishingly close to the approximation.
We used sum-of from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/EeqvRWXv.
Programming Praxis 
[…] Praxis – Numerical Integration By Remco Niemeijer In today’s Programming Praxis exercise we have to do some numerical integration. Let’s get […]
My Haskell solution (see http://bonsaicode.wordpress.com/2010/02/09/programming-praxis-numerical-integration/ for a version with comments):
int combine f a b n = sum $ map g [0..n - 1] where w = (b - a) / n lo i = a + w * i g i = w * combine (f $ lo i) (f $ lo i + w/2) (f $ lo i + w) intRect = int (\_ m _ -> m) intTrap = int (\l _ h -> (l + h) / 2) intSimp = int (\l m h -> (l + 4 * m + h) / 6) intAdapt m f a b epsilon = if abs (g10 - g5) < e then g10 else intAdapt m f a mid (Just e) + intAdapt m f mid b (Just e) where g5 = m f a b 5 g10 = m f a b 10 mid = (a + b) / 2 e = maybe 1e-7 id epsilon approxPi n = round $ intAdapt intSimp (recip . log) 2 n NothingPython code.
def rectangular(f, lo, hi, nsteps=1): area = 0 spread = float(hi-lo) for n in range(nsteps): x = lo + spread*(2*n+1)/(2*nsteps) area += f(x) return area*spread/nsteps def trapezoidal(f, lo, hi, nsteps=1): area = -(f(lo) + f(hi)) spread = float(hi-lo) for n in range(nsteps+1): x = lo + n*spread/nsteps area += 2*f(x) return area*float(hi-lo)/(2*nsteps) def simpson(f, lo, hi, nsteps=1): area = -(f(lo) + f(hi)) spread = float(hi-lo) for n in range(2*nsteps+1): x = lo + n*spread/(2*nsteps) area += 4*f(x) if n&1 else 2*f(x) return area*spread/(6*nsteps) def adaptive(f, lo, hi, integrator=simpson, eps=1e-3): lo,hi = float(lo),float(hi) area = 0 todo = [(lo, hi,integrator(f, lo, hi))] while todo: a,b,slice_area = todo.pop() mid = (a+b)/2 lohalf = integrator(f, a, mid) hihalf = integrator(f, mid, b) if abs(slice_area - (lohalf + hihalf)) > eps: todo.append((a,mid,lohalf)) todo.append((mid,b,hihalf)) else: area += lohalf + hihalf return area # tests def cube(x): return x*x*x def invlog(x): from math import log return 1/log(x) print rectangular(cube, 0, 1, 10000) print trapezoidal(cube, 0, 1, 10000) print simpson(cube, 0, 1, 10000) print adaptive(invlog, 2, 1e21, eps=1e-7)output:
0.24999999875
0.2500000025
0.25
2.11272694866e+19
A bit late to the party, but my solution in available here.
Github account moved, so gist disappeared. Here’s my submission (in Python 2.x):
from __future__ import division from math import log def dx(a, b, n): return (b - a) / n def rectangle(f, a, b, n=1): x = lambda k: a + (k + 0.5) * dx(a, b, n) return dx(a, b, n) * sum(f(x(k)) for k in xrange(n)) def trapezoid(f, a, b, n=1): x = lambda k: a + k * dx(a, b, n) return dx(a, b, n) * (f(a) + f(b) + sum(0.5 * (f(x(k)) + f(x(k + 1))) for k in xrange(1, n))) def simpson(f, a, b, n=1): x = lambda k: a + k * dx(a, b, 2 * n) return dx(a, b, 2 * n) / 3 * (f(a) + f(b) + 2 * sum(f(x(k)) for k in xrange(2, 2 * n, 2)) + 4 * sum(f(x(k)) for k in xrange(1, 2 * n, 2))) def adaptive_quad(f, a, b, quad=simpson, eps=1e-7): int_with5 = quad(f, a, b, 5) int_with10 = quad(f, a, b, 10) m = (a + b) / 2 if abs(int_with5 - int_with10) < eps: return int_with10 else: return (adaptive_quad(f, a, m, quad, eps) + adaptive_quad(f, m, b, quad, eps)) def approx_pi(b): return adaptive_quad(lambda x: 1 / log(x), 2, b) if __name__ == "__main__": cube = lambda x: pow(x, 3) print rectangle(cube, 0, 1, 10000) print trapezoid(cube, 0, 1, 10000) print simpson(cube, 0, 1, 10000) print approx_pi(1e21)