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.
C o m p . S c i .
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.
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)
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 {ax2+bx+c}
--i.e. area under quadratic
= [x(ax2/3 + bx/2 +c)]0..2w
= 2w(4aw2/3 + bw +c)
= 2w(8aw2 + 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: