Integration

Integration is the operation of finding the area under a curve, f(x), between given limits, from x=lo to x=hi. If the form of f(x) is suitable, and if we are smart enough, it may be possible to solve the integration analytically to get an exact answer. If this is not possible, f(x) can be integrated numerically to give an approximate answer. Such a computation is known as quadrature. The errors introduced in the process depend upon the properties of f(x) and on the method of integration. Common methods require that f(x) and its derivatives are continuous and that its high-order derivatives are "small".

Some simple methods of numerical integration are illustrated below.

Rectangle Rule

The range [lo,hi] is divided up into N equally sized intervals of width (hi-lo)/N. In a given interval, f(x) is approximated by its value at the centre of an interval, hence rectangles:

function rectangle(f, lo, hi, N)
 { var width = (hi-lo)/N;
   var sum = 0;

   var i;
   for(i=0; i < N; i++)
      sum += f( lo+(i+0.5)*width );
      // f() at centre of i-th interval

   return sum*width;
 }

Trapezoidal Rule

C
o
m
p
.
S
c
i
.

The range [lo,hi] is divided up into N intervals, as before. A straight-line approximation for f(x) is used in each interval, i.e. the area under f(x) is approximated by a series of trapeziums. Note that the area of a trapezium is its width multiplied by the average of the two parallel sides.

function trapezoidal(f, lo, hi, N)
 { var width = (hi-lo)/N;
   var sum = (f(lo)+f(hi))/2;

   var i;
   for(i=1; i < N; i++)
      sum += f( (lo*(N-i) + hi*i)/N );

   return sum*width;
 }

Simpson's Rule

A straight line can be fitted through any two points and a quadratic can be fitted through three points.

  |         |         |
  |      **q***       |
  |  ****   |  ****   |      * ax²+bx+c
  |*        |      ****r
 p*         |         |
  |         |         |
  |         |         |
  |         |         |
  |         |         |
  |         |         |
 -.---------------------------> x
  |         |         |
  0         w        2w

Fitting ax²+bx+c through (0,p), (w,q) & (2w,r) gives us some constraints:

         c = p         -- (0,p)
 aw²+ bw+c = q         -- (w,q)
4aw²+2bw+c = r         -- (2w,r)

We could solve the constraints to find a, b & c, but there is no need to do so directly...

The area under a quadratic is easy to work out:

∫_x=0..2w {ax²+bx+c} --i.e. area under quadratic

= [x(ax²/3 + bx/2 +c)]_0..2w

= 2w(4aw²/3 + bw +c)

= 2w(8aw² + 6bw + 6c)/6 --Note r+4q = 8aw+6bw+5c

= w*(r + 4q + c)/3

= w(r + 4q + p)/3

So, if the range [lo,hi] is divided into an even number of intervals, giving an odd number of boundaries, the area under f(x) can approximated by a sequence of quadratics. Code for the method can be optimised by realising that one needs a sum of f(x) evaluated at the boundaries of the intervals, weighted by four at odd-numbered internal boundaries and by two at even-numbered internal boundaries.

function Simpson(f, lo, hi, N)
// PRE: N is even
 { var width = (hi-lo)/N;
   var sum = f(lo)+f(hi);

   var i, odd=true;
   for(i=1; i < N; i++)
    { sum += f(lo+i*width) * (odd ? 4 : 2);
      odd = !odd;
    }

   return sum*width/3;
 }

Simpson's method is much more accurate than the rectangle and trapezoidal rules, for a given number of intervals.

Complexity

The rectangle rule evaluates f(x) N times. The trapezoidal rule and Simpson's method evaluate f(x) N+1 times.

Demonstration

The HTML FORM demonstrates the integration of a function, f(x), between x=lo and x=hi in N intervals. In this particular case f(x) is a polynomial so it is easy to do the integration analytically to find the correct answer for comparison with the numerical methods; this is not possible in general. Change f(x), lo, hi and the number of intervals, and experiment:

— 1999 L.A.