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 
Programming Praxis 
[…] Praxis – Numerical Integration By Remco Niemeijer In today’s Programming Praxis exercise we have to do some numerical integration. Let’s get […]
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 NothingPython 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
A bit late to the party, but my solution in available here.
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)#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 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 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));
}
int main() {
float (invlog_ptr)(float) = &invlog;
float (cube_ptr)(float) = &cube;
float limit = pow(10, 21);
}