Gaussian Quadrature Methods


Gaussian Quadrature Methods

The three techniques to solve integral equations that we have implemented so far in this chapter divide the integration domain into constant- sized increments of D x . This is acceptable for smoothly varying functions but may be inefficient or even inaccurate for functions with widely varying slopes. In this case, you would want to use a variable step size with small D x increments in regions of high gradients and larger D x increments in smoothly varying regions.

The trapezoidal, Simpson's rule, and midpoint methods are really specialized cases of a more general class of integration techniques called Gaussian quadrature formulas. Gaussian quadrature formulas assume that an integral equation can be approximated by a summation of weights multiplied by function evaluations at various x locations along the integration domain, as shown in Eq. (21.21).

Equation 21.21

graphics/21equ21.gif


The x locations at which the function is evaluated are called the abscissas . The coefficients in front of the function evaluations are called the weights . The determination of the weights depends on the choice for the W ( x ) function. There are a number of classical forms for W ( x ), each designed to be used for a certain type of integral equation. In this section, we will work with a commonly used form called the Gauss-Legendre formula, which assumes that W ( x ) = 1.

The implementation of the Gauss-Legendre formula requires two methods, one to compute the weights and abscissas and a second to integrate a specified function. You could alternately use a hard-coded array of weight and abscissa data, but evaluating this data directly gives you more flexibility in choosing the number of points to be considered in the integration domain.

Without going into all of the details, the computeGaussWeights() method uses Newton's method to determine the roots of the orthogonal polynomial equation associated with the Gauss-Legendre formula. The weights and abscissas are computed from the polynomial roots and are normalized for an integration ranging from x = 0 to x = 1. The weights and abscissas can be applied to any integration domain by multiplying the resulting integral value by b a , where a and b are the limits of integration.

The computeGaussWeights() method source code follows .

 public static void computeGaussWeights(              double x[], double w[], int numPoints) {   int i,j,m;   double z, zOld, root, p3, p2, p1;   double EPS = 1.0e-10;   //  The weights and abscissa values are normalized   //  for an integration range between 0 and 1. The   //  values are symmetric about x=0.5 so you only   //  have to compute half of them.   m = (numPoints + 1)/2;   //  Compute the Legendre polynomial, p1, at point   //  z using Newton's method.   for(i=0; i<m; ++i) {     z = Math.cos(Math.PI*(i+0.75)/(numPoints+0.5));     do {       p1 = 1.0;       p2 = 0.0;       for(j=0; j<numPoints; ++j) {         p3 = p2;         p2 = p1;         p1 = ((2.0*(j+1) - 1.0)*z*p2 - j*p3)/              (j+1.0);       }       root = numPoints*(z*p1 - p2)/(z*z - 1.0);       zOld = z;       z = zOld - p1/root;     } while ( Math.abs(z-zOld) > EPS);     //  Compute the weights and abscissas. The values     //  are symmetric about the point z=0.5     x[i] = 0.5*(1.0 - z);     x[numPoints-1-i] = 0.5*(1.0 + z);     w[i] = 1.0/((1.0- z*z)*root*root);     w[numPoints-1-i] = w[i];   } } 

Once we have the weight and abscissa data, the rest of the implementation is simple. The gaussLegendre() method integrates a function using the Gauss-Legendre weighting coefficients and abscissa data. The method takes the same four input arguments as the trapezoidal() , simpsonsRule() , and midpoint() methods.

After declaring some variables , the method enters an iteration loop. Starting with 10 points in the integration domain, the computeGaussWeights() method is called to obtain the 10-point Gauss-Legendre weights and abscissas. The integral value is then obtained by summing up the weights and function evaluations at the abscissas. The number of abscissa points is then doubled, the process is repeated, and the current solution is compared against the previous solution. If the two results are within the convergence tolerance, the value is returned. Otherwise, the number of points is again doubled and the process repeats itself.

The gaussLegendre() method source code is ”

 public static double gaussLegendre(Function function,      double start, double end, double tolerance) {   double dx = end - start;   double value = 1.0e+10;   double oldValue;   int numPoints = 10;   int maxPoints = 100000;   double x[] = new double[maxPoints];   double w[] = new double[maxPoints];   //  Compute the integral value adding points   //  until the solution converges or maxPoints   //  is reached.   while(true) {     //  Calculate the Gauss-Legendre weights and     //  abscissas     computeGaussWeights(x, w, numPoints);     //  Compute the integral value by summing up the     //  weights times the function value at each     //  abscissa     oldValue = value;     value = 0.0;     for(int i=0; i<numPoints; ++i) {      value += w[i]*function.getValue(start + x[i]*dx);     }     value *= dx;     //  If the solution is converged, return the     //  value. Otherwise, double the points in the     //  integration domain and try again.     if ( Math.abs(value-oldValue) <=          tolerance*Math.abs(oldValue) ) {       return value;     }     numPoints *= 2;     if (numPoints >= maxPoints) break;   }   //  Solution did not converge   return 12345.6; } 

The Gaussian quadrature technique can be applied to both proper and improper integrals so we use the gaussLegendre() method to solve Eq. (21.10) and Eq. (21.20). The driver class is called TestGauss and its source code is:

 import TechJava.MathLib.*; public class TestGauss {   public static void main(String args[]) {     double value;     //  Integrate the function y = x*ln(x) from     //  x=2 to x=4 using 10 points in the integration     LogFunction function = new LogFunction();     value =         Integrator.gaussLegendre(function,2.0,4.0,1.0e-4);     System.out.println("value = " + value);     //  Integrate the function y = 1/sqrt(x-2) from     //  x=2 to x=3 using 10 points in the integration     SqrtFunction2 function2 = new SqrtFunction2();     value =        Integrator.gaussLegendre(function2,2.0,3.0,1.0e-4);     System.out.println("value = " + value);   } } 

Output ”

 value = 6.70406052759485 value = 1.9998289705958587 

The gaussLegendre() method successfully reproduced the exact solutions to within four decimal places for both the proper and improper integrals. The logarithmic integral converged on the first iteration. The improper integral took longer ”10 iterations ”but this was somewhat more efficient than the midpoint method that required 14 iterations to converge.



Technical Java. Applications for Science and Engineering
Technical Java: Applications for Science and Engineering
ISBN: 0131018159
EAN: 2147483647
Year: 2003
Pages: 281
Authors: Grant Palmer

flylib.com © 2008-2017.
If you may any questions please contact us: flylib@qtcs.net