7.2 The Trapezoidal Algorithm
7.2 The Trapezoidal Algorithm As we can see in Figure 7-1, there is error at the top of each rectangle, since the rectangles don't fit exactly under the curve. These errors are smaller for narrower rectangles, but with a greater number of rectangles, we'll have a greater number of roundoff errors. Figure 7-2 shows how we can decrease the errors by using trapezoids instead of rectangles. Figure 7-2. Using trapezoids to estimate the area under a curve defined by a function f ( x ). Trapezoids are handy because there is a simple formula for computing their areas. If the heights of the two sides of a trapezoid are y 1 and y 2 , and its width is h, then the formula for the area A is Or, in other words, it is the width times the average of the two heights. [1]
That's the basic idea behind the trapezoidal integration algorithm. It simply fits a number of equal-width trapezoids under the curve of the integrand f ( x ). The heights of a trapezoid's sides are given by the integrand?aif the sides are at x 1 and x 2 , then the heights are f ( x 1 ) and f ( x 2 ). If we want to use n trapezoids to estimate the total area under the curve in the interval [ a, b ], then their common width is In theory, increasing n will increase the accuracy of the area estimate. However, in reality, a very large n (and hence a very small h ) will cause so many roundoff errors to accumulate that the accuracy will actually decrease. We will see this occurring in Program 7 §C1. Listing 7-1a shows the class TrapezoidalIntegrator in package numbercruncher.mathutils that implements the trapezoidal integration algorithm. Listing 7-1a The class that implements the trapezoidal integration algorithm. package numbercruncher.mathutils; /** * Function integrator that implements the trapezoidal algorithm. */ public class TrapezoidalIntegrator implements Integrator { /** the function to integrate */ private Evaluatable integrand; /** * Constructor. * @param integrand the function to integrate */ public TrapezoidalIntegrator(Evaluatable integrand) { this.integrand = integrand; } /** * Integrate the function from a to b using the trapezoidal * algorithm, and return an approximation to the area. * (Integrator implementation.) * @param a the lower limit * @param b the upper limit * @param intervals the number of equal-width intervals * @return an approximation to the area */ public float integrate(float a, float b, int intervals) { if (b <= a) return 0; float h = (b - a)/intervals; // interval width float totalArea = 0; // Compute the area using the current number of intervals. for (int i = 0; i < intervals; ++i) { float x1 = a + i*h; totalArea += areaOf(x1, h); } return totalArea; } /** * Compute the area of the ith trapezoidal region. * @param x1 the left bound of the region * @param h the interval width * @return the area of the region */ private float areaOf(float x1, float h) { float x2 = x1 + h; // right bound of the region float y1 = integrand.at(x1); // value at left bound float y2 = integrand.at(x2); // value at right bound float area = h*(y1 + y2)/2; // area of the region return area; } } The constructor takes an Evaluator argument, and so the argument can be a Function, InterpolationPolynomial , or RegressionLine object. Its second argument is the number of intervals. Method integrate() performs the integration by invoking method areaOf() for each region, summing the areas, and returning an approximation to the total area. Method areaOf() implements the formula for the area of a trapezoid. Class TrapezoidalIntegrator implements interface Integrator , which is also defined in package numbercruncher.mathutils . See Listing 7-1b. Listing 7-1b Interface integrator. package numbercruncher.mathutils; /** * Interface implemented by integrator classes. */ public interface Integrator { /** * Integrate the function from a to b, * and return an approximation to the area. * @param a the lower limit * @param b the upper limit * @param intervals the number of equal-width intervals * @return an approximation to the area */ float integrate(float a, float b, int intervals); } If we give Program 7 §C1 the command-line argument "trapezoidal," it instantiates a TrapezoidalIntegrator object. (We'll use the same program later when we examine Simpson's integration algorithm.) The noninteractive version of the program uses the object to compute an approximation to p ( since it integrates over the upper right quadrant of the circle. See Listing 7-1c. Listing 7-1c The noninteractive version of Program 7 §C1 computes an approximation to the value of p with the trapezoidal integration algorithm when given the command-line argument "trapezoidal." package numbercruncher.program7_1; import numbercruncher.mathutils.Function; import numbercruncher.mathutils.Integrator; import numbercruncher.mathutils.TrapezoidalIntegrator; import numbercruncher.mathutils.SimpsonsIntegrator; import numbercruncher.mathutils.Epsilon; import numbercruncher.mathutils.AlignRight; /** * PROGRAM 7-1: Integration * * Demonstrate numerical integration algorithms. */ public class Integration { private static final int MAX_INTERVALS = Integer.MAX_VALUE/2; private static final float TOLERANCE = 100*Epsilon.floatValue(); private static final float FROM_LIMIT = 0; private static final float TO_LIMIT = 1; private static final int TRAPEZOIDAL = 0; private static final int SIMPSONS = 1; private static final String ALGORITHMS[] = {"Trapezoidal", "Simpson's"}; /** * Main program. * @param args the array of runtime arguments */ public static void main(String args[]) { String arg = args[0].toLowerCase(); int algorithm = arg.startsWith("tra") ? TRAPEZOIDAL : SIMPSONS; System.out.println("ALGORITHM: " + ALGORITHMS[algorithm]); System.out.println(); // The function to integrate. Function integrand = new Function() { public float at(float x) { return (float) Math.sqrt(1 - x*x); } }; integrate(algorithm, integrand); } /** * Do the integration with either the trapezoidal algorithm * or Simpson's algorithm. * @param algorithm 0 for trapezoidal, 1 for Simpson's * @param integrand the function to integrate * @return an approximation to the area */ private static void integrate(int algorithm, Function integrand) { int intervals = 1; // number of intervals float area = 0; // total area float errorPct = Float.MAX_VALUE; // % error float prevArea; float prevErrorPct; Integrator integrator = (algorithm == TRAPEZOIDAL) ? (Integrator) new TrapezoidalIntegrator(integrand) : (Integrator) new SimpsonsIntegrator(integrand); AlignRight ar = new AlignRight(); ar.print("n", 5); ar.print("pi", 15); ar.print("% Error", 15); ar.underline(); do { prevArea = area; area = integrator.integrate(FROM_LIMIT, TO_LIMIT, intervals); float pi = 4*area; prevErrorPct = errorPct; errorPct = (float) Math.abs(100*(pi - Math.PI)/Math.PI); ar.print(intervals, 5); ar.print(pi, 15); ar.print(errorPct, 15); ar.println(); intervals *= 2; } while ((errorPct > TOLERANCE) && (errorPct < prevErrorPct) && (intervals < MAX_INTERVALS)); } } Output: ALGORITHM: Trapezoidal n pi % Error ----------------------------------- 1 2.0 36.338024 2 2.732051 13.036119 4 2.995709 4.643623 8 3.089819 1.648008 16 3.1232529 0.5837735 32 3.1351023 0.20659526 64 3.1392965 0.073087834 128 3.1407807 0.025845688 256 3.1413052 0.009149671 512 3.1414926 0.003184639 1024 3.1415575 0.0011204039 2048 3.1415784 4.525632E-4 4096 3.1415865 1.9453382E-4 8192 3.1415906 6.551914E-5 16384 3.1415942 4.8317346E-5 32768 3.1415963 1.1661924E-4 The program begins with one interval ( n = 1) and doubles the number for each iteration. Notice that it has three stopping criteria: (1) the number of intervals has exceeded the maximum, (2) the error percentage is below the tolerance, or (3) the error percentage has begun to rise. The error percentage starts to rise when the accumulated roundoff errors overwhelm the computation. Screen 7-1a shows screen shots of the first three iterations of the interactive version [2] of Program 7 §C1 using the trapezoidal algorithm, with one, two, and four intervals. This program also uses polynomial interpolation (see Chapter 6) to allow you to define an arbitrary integrand.
Screen 7-1a. The first three iterations of the interactive version of Program 7 §C1 using the trapezoidal algorithm. |
. 
3.1415926) in the following way. The area of the unit circle (its radius is 1) is p , and its equation is 
 |
| Top |
EAN: 2147483647
Pages: 141