Numerical Integration

February 9, 2010

Computing the definite integral of a function is fundamental to many computations. In many cases, the integral of a function can be determined symbolically and the definite integral computed directly. Another approach, which we shall examine in today’s exercise, is to compute the definite integral numerically, an operation that numerical analysts call quadrature. We will consider a function f for which we are to find the definite integral over the range a to b, with a < b.

Recall from your study of calculus that the definite integral is the area under a curve. If we slice the curve into hundreds or thousands or millions of vertical strips, we can approximate the area under the curve by summing the areas of the individual strips, assuming that each is a rectangle with height equal to the value of the curve at the midpoint of the strip; this is, obviously, the rectangular method of quadrature. The trapezoidal method is similar, but uses the average of the value of the curve at the two ends of the strip; the strip then forms a trapezoid with parallel sides, a flat bottom on the horizontal axis, and a slanted top that follows the curve. A final method, called Simpson’s method, fits a parabola to the top of the strip by taking the height of the curve as its value at the left side of the strip plus its value at the right side of the strip plus four times its value at the midpoint of the strip. Accuracy improves as we move from the rectangular method (one point) to the trapezoidal method (two points) to Simpson’s method (three points), because the approximation more closely fits the curve, and also as the number of slices increases. The trapezoid rule is illustrated above right.

Many curves are characterized by large portions that are nearly flat, where the approximations work well, with a few small portions that change rapidly, where the approximations fail. If we increase the number of slices, there is little effect on the flat parts of the curve, so much of the work is wasted. Adaptive quadrature works by measuring the approximations at various slices and narrowing the slices where the curve is changing rapidly, allowing us to control the approximation precisely. Adaptive quadrature is a recursive process; if approximations with n and n/2 slices are close enough, recursion stops, otherwise the range is split in two and adaptive quadrature is applied to each of the two halves separately.

Your task is to write functions that perform numerical integration by the rectangular, trapezoidal and Simpson’s methods, as well as a function that performs adaptive quadrature. Use your functions to calculate the logarithmic integral \mathrm{li}(x) = \int_2^x \frac{dt}{\log{t}} for x = 1021, where the logarithmic integral is an approximation of the number of primes less than its argument. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.


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


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

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

Facebook photo

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

Connecting to %s

%d bloggers like this: