<- previous index next ->
Numerical integration could be simple summation, but
in order to get reasonable step size and good accuracy
we need better techniques.
integral f(x) dx from x=a to x=b with step size h=(b-a)/n
using simple summation would be computed:
area = (b-a)/n * ( sum i=0..n-1 f(a+i*h) )
"h" has become the approximation to "dx" using "n" steps.
A larger value of "n" gives a smaller value of "h" and
up to where roundoff error grows, a better approximation.
This simple calculation did not use the value f(b).
The Trapezoidal rule approximates the area between x and x+h
as h * (f(x)+f(x+h))/2 , base times average height.
Summing the areas we note that f(a) and f(b) are used once
while the other intermediate f(x) are used twice, thus:
area = (b-a)/n * ( (f(a)+f(b))/2 + sum i=1..n-1 f(a+i*h) )
note: i=0 is f(a) i=n is f(b) thus sum index range 1..n-1
h = (b-a)/n
cutting h in half generally cuts the error in half for large n
In order to get better accuracy with fewer evaluations of
the function, f(x), we have found a way to choose the values
of x and to apply a weight, w(x) to each evaluation of f(x).
The integral is evaluated using:
area = sum i=1..n w(i)*f(x(i))
The w(i) and x(i) are computed using the Legendre polynomials
covered in the previous lecture. The numerical analysis shows
that using n evaluations we obtain accuracy about equal to
fitting the f(x(i)) with an nth order polynomial and accurately
computing the integral using that nth order polynomial.
Some values of weights w(i) and ordinates x(i) -1 < x < 1 are:
x[1]= 0.0000000000000E+00, w[1]= 2.0000000000000E+00
x[1]= -5.7735026918963E-01, w[1]= 1.0000000000000E-00
x[2]= 5.7735026918963E-01, w[2]= 1.0000000000000E-00
x[1]= -7.7459666924148E-01, w[1]= 5.5555555555555E-01
x[2]= 0.0000000000000E+00, w[2]= 8.8888888888889E-01
x[3]= 7.7459666924148E-01, w[3]= 5.5555555555555E-01
x[1]= -8.6113631159405E-01, w[1]= 3.4785484513745E-01
x[2]= -3.3998104358486E-01, w[2]= 6.5214515486255E-01
x[3]= 3.3998104358486E-01, w[3]= 6.5214515486255E-01
x[4]= 8.6113631159405E-01, w[4]= 3.4785484513745E-01
x[1]= -9.0617984593866E-01, w[1]= 2.3692688505618E-01
x[2]= -5.3846931010568E-01, w[2]= 4.7862867049937E-01
x[3]= 0.0000000000000E+00, w[3]= 5.6888888888889E-01
x[4]= 5.3846931010568E-01, w[4]= 4.7862867049937E-01
x[5]= 9.0617984593866E-01, w[5]= 2.3692688505618E-01
x[1]= -9.3246951420315E-01, w[1]= 1.7132449237916E-01
x[2]= -6.6120938646626E-01, w[2]= 3.6076157304814E-01
x[3]= -2.3861918608320E-01, w[3]= 4.6791393457269E-01
x[4]= 2.3861918608320E-01, w[4]= 4.6791393457269E-01
x[5]= 6.6120938646626E-01, w[5]= 3.6076157304814E-01
x[6]= 9.3246951420315E-01, w[6]= 1.7132449237916E-01
Using the function gaulegf a typical integration could be:
double x[9], w[9] for 8 points
a = 0.5; integrate sin(x) from a to b
b = 1.0;
n = 8;
gaulegf(a, b, x, w, n); x's adjusted for a and b
area = 0.0;
for(j=1; j<=n; j++)
area = area + w[j]*sin(x[j]);
The following programs, gauleg for Gauss Legendre integration,
computes the x(i) and w(i) for various values of n. The
integration is tested on a constant f(x)=1.0 and then on
integral sin(x) from x=0.5 to x=1.0
integral exp(x) from x=0.5 to x=5.0
integral ((x^x)^x)*(x*(log(x)+1)+x*log(x)) from x=0.5 to x=5.0
Note how the accuracy increases with increasing values of n,
then, no further accuracy is accomplished with larger n.
Also note, the n where best numerical accuracy is achieved
is a far smaller n than where roundoff error would be
significant.
Choose your favorite language and study the .out file
then look over the source code.
gauleg.c
gauleg_c.out
gauleg.f90
gauleg_f90.out
gauleg.for
gauleg_for.out
gauleg.java
gauleg_java.out
gauleg.adb
gauleg_ada.out
Multidimensional integration extends by using
volume = volume + wx[i]*wy[j]*f(xx[i],yy[j])
as shown in:
int2d.c
int2d_c.out
Additional reference books include:
"Handbook of Mathematical Functions" by Abramowitz and Stegun
A classic reference book on numerical integration and differentiation.
"Numerical Recipes in Fortran" by Press, Teukolsky, Vetterling and Flannery
There is also a "Numerical recipes in C" yet in my copy, there have
been a number of errors in the C code due to poor translation from Fortran.
Many very good numerical codes for general purpose and special purposes.
<- previous index next ->