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

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
for n in range(nsteps):
area += f(x)

def trapezoidal(f, lo, hi, nsteps=1):
area = -(f(lo) + f(hi))
for n in range(nsteps+1):
area += 2*f(x)

return area*float(hi-lo)/(2*nsteps)

def simpson(f, lo, hi, nsteps=1):
area = -(f(lo) + f(hi))
for n in range(2*nsteps+1):
area += 4*f(x) if n&1 else 2*f(x)

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)



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

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: