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.

About these ads

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

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

    Python 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

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 634 other followers

%d bloggers like this: