Numerical Integration

Assignment :

Read Chpt 15, Chpt. 18.8, and Web Notes on Numerical Differentiation. Finish Hw13

New Fortran:


Basic Numerical Integration

You could spend most of a semester learning all of the methods to evaluate integrals (single and multiple) on a computer. I'm just going to cover two simple methods for evaluation of definite integrals.

Computers don't know anything about Calculus in their CPUs. They simply do basic arithmetic. One way to deal with this is to go back to thinking of Integrals as sums of areas in small intervals.

Another approach is to obtain exact values of the function f(x) at some number of points, then fit curves to these points for which we know the value of the integral. The trapezoidal method and Simpson's rule can be regarded as variations on this curve fitting approach. Keep in mind that these involve exact fits to the points not "best" fits (although that can be made to work also.

Trapezoidal Rule.

The basic approach hear is to convert the problem of finding the area under some curve, to one of finding the combined areas of a set of trapezoids. The upper side is obtained by approximating the actual curve with a straight line, fit to two points on the curve.

We start with some x-y curve where y= f(x). Next we select a set of n values of x ( xi, I=1,n) from x1 = a through xn = b (bounds of the definite integral) and evaluate the function to obtain corresponding values of y (yi, I=1,n). The integral (area under the curve) is approximated by the sum of the areas of n-1 trapezoids bounded by lines connecting the points (xi,0.), (xi, yi), (xi+1, yi+1), and (xi+1,0.). The area of one such trapezoid is:

Ai = 0.5 ( yi+1 + yi ) ( xi+1 - xi )

so the integral can be approximated as

if the x sample points are spaced evenly (x2-x1=x3-x2=... ) then this simplifies to

How much error have we introduced with this approximation? Within a given interval we can think of our function as being approximated by a Taylor series:

The Trapezoidal rule only covers the right hand side through the first power of x-xi. The largest errors involve first integrating the quadratic term between two points giving a local error proportional to the cube of the x spacing and summing all of these errors gives a total error proportional to the square of the x spacing . If you halve the value of x spacing (double the number of points at which you evaluate f(x)), then you drop the error in your approximation to the integral by a factor of 4.

You can look at some examples of the Trapezoidal rule for integration in trapezoid.f, trapz1.f, and trapz2.f. trapz1.f needs the files trapcom.h and trapcom1.h to compile properly. You can get all of these files by typing:

cp ~jhm/201/integrate/* .

These examples have been programmed to illustrate passing information between program units (main program, subroutines, and functions) in three different ways. You should by now be used to argument lists and modules. The third is an old Fortran statement called a COMMON block. It is similar to a MODULE but more difficult to use. More on COMMON blocks in the next section.

Simpson's Rule

Simpson's rule starts by approximating a segment of the function f(x) as quadratic curve. The basic method operates on groups of 3 evenly spaced points xi-1, xi, and xi+1. We want to find coefficients a0, a1, and a2 for the quadratic

such that:

After a little algebra, the following values are obtained for the coefficients:

We are now in a position to approximate the integral of f(x) between xi-1 and xi+1.


or in terms of y=f(x):

Notice 3 things about this formula:

  1. You need to evaluate the function for at least 3 values of x;
  2. To integrate with many points (more accuracy) you must add 2 points for each additional integral of a quadratic fit (third point comes from the end of the last set);
  3. This means that the application of Simpson's rule requires evaluation of the function at an odd number of points (> 3).

The general summation formula for Simpson's rule becomes (for n and odd integer):

If you store the function values in an array "y" such that y(1) is f(a) and y(n) is f(b) then a simple Fortran 90 expression for the above equation is:

simpint = .3333333*(y(1) + y(n) + 4*sum(y(2:n-1:2)+2*sum(3:n-2:2))

Back to the Table of Contents / Home

Written and Maintained by John Mahaffy :