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

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***       |
  |  ****   |  ****   |      * ax2+bx+c
  |*        |      ****r
 p*         |         |
  |         |         |
  |         |         |
  |         |         |
  |         |         |
  |         |         |
 -.---------------------------> x
  |         |         |
  0         w        2w

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

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

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

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: