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

Pages: 1 2

### 5 Responses to “Numerical Integration”

1. […] Praxis – Numerical Integration By Remco Niemeijer In today’s Programming Praxis exercise we have to do some numerical integration. Let’s get […]

2. Remco Niemeijer said

```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 Nothing
```
3. Mike said

Python code.

```def rectangular(f, lo, hi, nsteps=1):
area = 0
for n in range(nsteps):
area += f(x)

def trapezoidal(f, lo, hi, nsteps=1):
area = -(f(lo) + f(hi))
for n in range(nsteps+1):
area += 2*f(x)

return area*float(hi-lo)/(2*nsteps)

def simpson(f, lo, hi, nsteps=1):
area = -(f(lo) + f(hi))
for n in range(2*nsteps+1):
area += 4*f(x) if n&1 else 2*f(x)

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)

```

output:
0.24999999875
0.2500000025
0.25
2.11272694866e+19

4. Graham said

A bit late to the party, but my solution in available here.

5. Graham said

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)))

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: