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

Posted by programmingpraxis
Filed in Exercises
6 Comments »

6 Responses to “Numerical Integration”

  1. Programming Praxis – Numerical Integration « Bonsai Code said
    February 9, 2010 at 10:28 AM

    […] 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
    February 9, 2010 at 10:28 AM

    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
    April 15, 2010 at 8:40 PM

    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
    April 13, 2011 at 12:40 AM

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

  5. Graham said
    September 8, 2011 at 4:41 PM

    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)
  6. Eduardo Jesus said
    January 26, 2024 at 2:19 AM

    #include <stdio.h>
    #include <math.h>

    float invlog(float x) { return 1.0 / log(x); }

    float cube(float x) { return pow(x, 3); }

    float trapezoid(float (*f)(float), float a, float b, int steps) {
    if (steps <= 0) { printf(“Invalid number of intervals\n”); return 0; }

    float delta = fabs(b - a) / (float)steps;
float integral = delta * (f(a) + f(b)) / 2;

for (int i = 1; i < steps; ++i) {
    float x_n = a + i * delta;
    integral += delta * f(x_n);
}

return integral;

    }

    float simpson(float (*f)(float), float a, float b, int steps) {
    if (steps <= 0 || steps % 2 != 0) { printf(“Invalid steps. Must be even and positive\n”); }

    float delta = fabs(b - a) / (float)steps;
float integral = (delta / 3) * (f(a) + f(b));

for (int i = 1; i < steps; ++i) {
    float x_n = a + i * delta;
    integral += (delta / 3) * ((i % 2 != 0) * 4 * f(x_n) + (i % 2 == 0) * 2 * f(x_n));
}

return integral;

    }

    float quadrature(float (*f)(float), float a, float b) {
    float eps = 1e-7;
    float integral = 0;
    float midpoint = fabs(b – a) / 2.0;
    float error_margin = fabs(simpson(f, a, midpoint, 2) + simpson(f, midpoint, b, 2) – simpson(f, a, b, 2));

    if (error_margin > eps) {
    integral += quadrature(f, a, midpoint);
    integral += quadrature(f, midpoint, b);
} else {
    integral += simpson(f, a, b, 2);
}
return integral;

    }

    int main() {
    float (invlog_ptr)(float) = &invlog;
    float (    cube_ptr)(float) = &cube;
    float limit = pow(10, 21);

    float integral = quadrature(invlog_ptr, 2, limit);
float simp = simpson(cube_ptr, 0, 1, 10000);
float traps = trapezoid(cube_ptr, 0, 1, 10000);

printf("%f\n", integral);
printf("%f\n", simp);
printf("%f\n", traps);

invlog_ptr = NULL;
cube_ptr = NULL;

return 0;

    }

Leave a comment