# Programming Praxis

A collection of etudes, updated weekly, for the education and enjoyment of the savvy programmer

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

### Like this:

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
5 Comments »

### 5 Responses to “Numerical Integration”

1. 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)
```
%d bloggers like this:
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 1, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "github"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "english_confidence": 0.9074325561523438, "perplexity_flag": "middle"}