Trapezoidal Algorithms

Trapezoidal algorithms are a general class of solution methods used to estimate the area under the curve f ( x ) from x = a to x = b . The algorithms we will discuss in this section will be for proper integrals whose function f ( x ) can be evaluated at all points between a and b . A typical curve is shown in Figure 21.1.

Figure 21.1. Two-point trapezoidal method

The simplest way to approximate the area under f ( x ) is to draw a line between f ( a ) and f ( b ) and compute the area under the resulting trapezoidal figure.

Equation 21.2

The D x term in Eq. (21.2) is the width of the trapezoid and for the two-point trapezoid method is given by the expression D x = b “ a . Of course for most curves, Eq. (21.2) is not a very good approximation . Looking at the curves in Figure 21.1, we are missing the area between the actual curve and the straight line. To improve the accuracy of the solution, we can divide the range of integration into two sections at the midpoint x location between x = b and x = a . In this case, we have to compute the area of two trapezoids, and the approximate value of the integral becomes the following.

Equation 21.3

In Eq. (21.3), the width of each of the two trapezoids is half as large as when only one trapezoid was considered so D x = ( b “ a )/2. You can easily see that Eq. (21.3) can be extended to a general form with an arbitrary number of divisions of D x .

Equation 21.4

In Eq. (21.4) D x = ( b “ a )/( N “ 1) where N is the number of trapezoids considered between x = a and x = b . Eq. (21.4) is called the trapezoidal method . It is a second-order method because its error is proportional to 1/ N ² . One thing to note about the trapezoidal method is that the D x increments are uniform over the entire integration range.

But how do you know how many times to divide up the integration range for a given integral equation? You keep adding points until the computed integral value changes less than a specified convergence tolerance. If you are clever about how you add points, you can reuse the values from the previous approximation. For example, consider an integration range between x = 0 and x = 1. If three points are used the expression to evaluate is the following.

Equation 21.5

Keep in mind when looking at Eq. (21.5) that D x = ( b “ a )/2 = (1 “ 0)/2 = 1/2. If you add two additional points for the next approximation, the expression becomes Eq. (21.6).

Equation 21.6

Adding four additional points for the next approximation yields Eq. (21.7).

Equation 21.7

You can see that a pattern is developing. If you start by adding one point to the two-point trapezoidal approximation and then double the number of additional points with each iteration, the next approximation to the integral is one half the previous approximation plus the factors due to the additional points. An equation for a general update to an integral approximation is shown in Eq. (21.8).

Equation 21.8

where

Equation 21.9

Now let's write a method to implement the trapezoidal method. As we did with our system of equation and differential equation solvers, the trapezoidal method will be defined as a public , static method and its name will be trapezoidal() . The method will be defined in a class named Integrator . This is a mathematical routine so the Integrator class will be placed in the TechJava.MathLib package.

The trapezoidal() method takes four arguments. The first is a Function object that represents the mathematical function being integrated. We'll discuss the Function class in more detail later in this section. The second and third arguments define the range of integration. The fourth argument defines the convergence criteria for the iteration that will be performed.

After defining some variables , the trapezoidal() method computes an approximation to the integral using the two-point trapezoidal algorithm shown in Eq. (21.2). Subsequent (and better) approximations are produced by successively adding points to the integration range according to Eq. (21.8). When the change in the computed value of the integral is less than the specified convergence tolerance, the method returns. If the iteration does not converge, the value 12345.6 is returned.

The trapezoidal() method source code is ”

 package TechJava.MathLib; public class Integrator {     public static double trapezoidal(Function function,            double start, double end, double tolerance) {     double dx = end - start;     double value, oldValue, sum, scale;     int maxIter = 20;     //  Compute an initial value for the integral using     //  the two-point trapezoidal rule.     value = 0.5*dx*( function.getValue(start) +                        function.getValue(end) );     // Subdivide the integration range by adding more     // points, first 1 point, then two more, then four     // more and so on. Recompute the integral value     // each time. If the solution converges, return the     // value.     for(int i=0; i<maxIter; ++i) {       sum = 0.0;       scale = Math.pow(0.5,i+1);       for(int j=0; j<Math.pow(2,i); ++j) {         sum += function.getValue(                        start + scale*(2*j + 1.0)*dx);       }       oldValue = value;       value = 0.5*value + scale*dx*sum;       if ( Math.abs(value-oldValue) <=            tolerance*Math.abs(oldValue) ) {          return value;       }     }     //  Solution did not converge     return 12345.6;   } } 

The trapezoidal() method makes use of a reference to a Function object. The Function class is defined to be the abstract parent class of classes that represent mathematical functions. The class declares one method named getValue() that simply returns the value of the function evaluated at a given independent variable value. The Function class listing is ”

 package TechJava.MathLib; public abstract class Function {    abstract double getValue(double u); } 

Let's test the trapezoidal() method by having it integrate two mathematical functions. Both functions have an exact, closed-form solution, so an exact value of the integral can be computed. The two functions, their integral solutions, and the numerical values of their solutions are shown in Eq. (21.10) and Eq. (2.11). In Eq. (21.11), if a = 2, y = 7.797853348.

Equation 21.10

Equation 21.11

To use the trapezoidal() method on these functions, we will have to create two Function subclasses to represent the functions. The LogFunction class encapsulates Eq. (21.10). Its code listing is shown here.

 package TechJava.MathLib; //  This class represents the function //  y = u*ln(u) public class LogFunction extends Function {    public double getValue(double u) {       return u*Math.log(u);    } } 

The SqrtFunction class represents the expression shown in Eq. (21.11) and here is its source code.

 package TechJava.MathLib; //  This class represents the function //  y = sqrt(u^2 + a^2) public class SqrtFunction extends Function {    private double a;    public SqrtFunction(double a) {       this.a = a;    }    public double getValue(double u) {       return Math.sqrt(u*u + a*a);    } } 

All that remains is to write a driver program that will send the appropriate input to the trapezoidal() method and display the results. The TestTrap class performs this function. It simply creates a Function subclass object and sends a reference to the object to the trapezoidal() method. The computed integral value is then printed to standard output. The TestTrap class source code is shown next.

 import TechJava.MathLib.*; public class TestTrap {   public static void main(String args[]) {     double value;     //  Integrate the function y = x*ln(x) from     //  x=2 to x=4 with an error tolerance of 1.0e-4     LogFunction function = new LogFunction();     value = Integrator.trapezoidal(                          function,2.0,4.0,1.0e-4);     System.out.println("value = "+value);     //  Integrate the function y = sqrt(x^2 + 2^2) from     //  x=0 to x=3 with an error tolerance of 1.0e-4     SqrtFunction function2 = new SqrtFunction(2.0);     value = Integrator.trapezoidal(                          function2,0.0,3.0,1.0e-4);     System.out.println("value = "+value);   } } 

Output ”

 value = 6.704116936052852 value = 7.798005701125523 

The results from the trapezoidal() method are the same as the exact results to within four decimal places, which is the value to which the convergence tolerance was set.