## 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 10^{21} 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

[…] 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):

Python code.

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