Electromagnetics.Biz

Numerical Integration

Home > Reference > Numerical Methods > Numerical Integration

 

To Home Page  

To Home Page

     
History of Electromagnetics  

History of Electromagnetics

   

 

 
Education and Theory  

Education and Theory

   
Reference Information  

Reference Information

   
Tools and Calculators  

Tools and Calculators

   
Applications  

Applications

   
         
Companies, Products, and Services  

Companies

   
Organizations  

Organizations

   
Job Sites and Resources  

Jobs

   
Product Development  

Product Development

   
         
Publishing Information  

Publishing Information

   
Books and Magazines  

Books and Magazines

   
         
Other Resources  

Other Resources

   
         
   

About Me

   
   

About This Site

   
   

Contact Me

   
         
   

Site Map

   
   

Site Index

   
         
   

Spawn New Browser

   
         

Index of this page

        Gaussian Quadrature Integration

        Clenshaw-Curtis Integration

 

Gaussian Quadrature Integration

        In my own work, I have extensively used Gaussian quadratures for numerical integration over various domains.  The derivation of the formulas are available from a variety of sources including the papers listed below.  A sketch of one approach is given in the series called Numerical Recipes by Press, Flannery, Teukolsky and Vetterling (Cambridge University Press).  My own copy is the out-of-print Numerical Recipes in Pascal.  Formulas are also available in Handbook of Mathematical Functions 9th ed. by Abramowitz and Stegun (Dover Publications).

        Often people think that Gaussian quadrature is the only quadrature scheme and the best method of numerical integration.  The first assertion is false, and the second assertion may be false depending on your application.  There are a variety of quadrature schemes (see the Numerical Recipes reference for a list of a few of them).  In general, a quadrature scheme allows one to approximate the integral of a function, f(x), times another function, w(x), over some domain.  In Gauss-Legendre quadrature, w(x) = 1.  Other schemes allow w(x) to take on other values.

        If you are trying to integrate the product of an unknown function, f(x) and a known function, w(x), then you have two choices.  First, you could use Gauss-Legendre quadrature to integrate over the entire function, f(x)*w(x).  Second, you could also spend the time to develop or research a Gauss quadrature scheme based on the known function, w(x), and apply the quadrature formula to f(x).  In this manner, you could save computation time by not having to evaluate w(x) for each abscissa.  Although the preceding argument was applied to a one-dimensional argument (x), it also applies to multi-dimensional arguments.

        As for being the best method of numerical integration, Gaussian quadratures are perhaps the best method of integrating polynomials because they guarantee that they are exact for polynomials less than a specified degree.  For any other function besides a polynomial -- WATCH OUT!  Quadratures guarantee nothing!  When integrating with Gaussian quadratures, what is essentially happening is this.  First, you are curve fitting the "best" possible polynomial to the actual function.  Then, you are computing the area under the polynomial exactly.  This means that if you take your original function and manually curve-fit a polynomial to it, and if the curve fit does not match the original function very well, then your integration result will not be very accurate.  On the other hand, if the two agree, then the integration result will be accurate.

        This perspective points out why quadratures do not integrate over discontinuities very well.  To accurately model a discontinuity with a polynomial, one must begin with a high-degree polynomial to allow the freedom to change slope rapidly.  High-degree polynomials, however, tend to oscillate wildly around the discontinuity.  Low-degree polynomials can not change slope fast enough to model the discontinuity.  The result is that no polynomial, low-degree or high-degree, can sufficiently model a function with a step or slope discontinuity.  Therefore, no Gaussian quadrature can produce an accurate, reliable result.

 

Linear-Domain Gauss-Legendre Quadratures

        The following one-dimensional weights and abscissa come from Abramowitz and Stegun.  This table assumes that the domain extends over the interval [0,1].  The usual scheme is to have the domain vary over [-1,1].  The usual scheme can easily be recovered from the following table by deleting the divide by 2's and the additions by 0.5.  d is the highest degree of polynomial that the formula can integrate exactly.  N is the number of points.  w represents the weights.  x represents the abscissa.

 

        Linear-Domain Gauss Quadrature weights and abcissa

 

Triangular-Domain Gauss-Legendre Quadratures

        At least three authors provide formula for Gaussian integration over triangular domains:

                G. R. Cowper, "Gaussian Quadrature Formulas for Triangles", International Journal of Numerical Methods in Engineering, Vol 7, pp 405 - 408, 1973.

                J. N. Lyness and D. Jespersen, "Moderate Degree Symmetric Gaussian Quadrature Rules for the Triangle", International Journal for Numerical Methods in Engineering, Vol 21, pp 1129-1148, 1985.

                D. A. Dunavant, "High Degree Efficient Symmetrical Gaussian Quadrature Rules for the Triangle", International Journal of Numerical Methods in Engineering, Vol 21, pp 1129 - 1148, 1985.

Also of interest is

                P. C. Hammer and A. H. Stroud, "Numerical Integration over Simplexes and Cones", Mathematical Tables        and other Aides to Computation, Vol X, No 54, pp 130 - 139, April 1956.

 

When using Cowper, beware that zeta-3 for the 7-point formula, degree 5 differs from the result given in reference [5] of that paper.  Compare 0.05971... to 0.05961....  Of the most practical interest to users of quadratures is Lyness and Jespersen.  The weights and abscissa below come from this source.  With regard to Hammer and Stroud, a "simplex" is, for example, a triangle or a

tetrahedron; that is, a simplex is the simplest non-degenerate straight-sided object that can be made in the allowed vector space.  This paper is written in mathematese, which makes it harder for us mortal types.

        Integration over triangular domains is usually carried out in normalized coordinates.  To perform the integration, first map one vertex (vertex 1) to the origin, the second vertex (vertex 2) to point (1,0), and the third vertex (vertex 3) to (0,1).  This transformation is most easily accomplished by defining vector ra to point from vertex 1 to vertex 2 and vector rb to point from vertex 1 to vertex 3.  The x and y components of ra are given by rax and ray, respectively.  Then the substitution

        x = rax (z) + rbx (h)

        y = ray (z) + rby (h)

is made where z and h are coordinates in the normalized coordinate system.  Taking the Jacobian, J, one finds that

        J = rax rby - ray rbx

which just happens to equal twice the area of the triangle.  Hence, integration in the original and normalized coordinate systems can be related by

        Integral ( ) dx dy = 2A Integral ( ) dz dh

where A is the area of the triangle.

        Hence, the following weights and abscissa are presented for the normalized coordinates just described, and abscissa, x and y, vary from 0 to 1 subject to the restriction that x + y <= 1.

 

                Triangular-Domain Gauss Quadrature weights and abscissa

 

Tetrahedral-Domain Gauss-Legendre Quadratures

        Finally, there are available formulas for integration over tetrahedra.  See, for example,

                P. Keast, "Moderate-Degree Tetrahedral Quadrature Formulas", Computer Methods in Applied Mechanics and        Engineering, Vol 55, pp 339 - 348, 1986.

 

In this reference, the first weight for N = 8 appears to be in error by a sign.  The correct weight is -0.3932....

 

        The normalized coordinate system is similar to the one for integration over triangles, except that one begins with

        x = rax (z) + rbx (h) + rcx (x)

        y = ray (z) + rby (h) + rcy (x)

        z = raz (z) + rbz (h) + rcz (x)

With this definition, the Jacobian happens to be 6 times the volume, V, of the original tetrahedron with the result that

        Integral ( ) dx dy dz = 6V Integral ( ) dz dh dx

        Hence, the following weights and abscissa are presented for the normalized coordinates just described, and abscissa, x, y, and z vary from 0 to 1 subject to the restriction that x + y + z <= 1.

 

                Tetrahedral-Domain Gauss Quadrature weights and abscissa

 

Quadratures Over Circles and Spheres

        I am not aware of any published quadrature formulas over the area of a circle.  Integration over the circumference of a circle is easy -- just use points equally-spaced in angle and weighted uniformly.  If any one knows where to find information about integration over the area of a circle, the surface of a sphere, or the volume of a sphere, please email me.

 

Clenshaw-Curtis Integration

        Although fewer people use Clenshaw-Curtis integration than quadrature integration, there are some good reasons to consider them.  First, error bounds are available to estimate how far off the result is from the actual result.  This is in stark contrast to quadrature schemes which promise nothing unless the equation being integrated is a polynomial of less than a specified degree.  Second, Clenshaw-Curtis integration provides the ability to reuse previous results.

        For example, a common method of applying Gaussian quadratures is to integrate using, say, 8 points of evaluation, then repeating the integration with 16 points (twice as many as 8).  If the two results match, then the answer is taken as accurate.  The problem is that 24 total evaluation points are required.  With Clenshaw-Curtis integration, 8 points of integration can be used initially with the result compared to an error bound.  If the error is small enough, one might stop at that point.  A still-suspicious individual might compute the result for 16 points, but only 8 more evaluation points are required because the first 8 points are reused in the second calculation!  This means that the suspicious individual has to compute only 16 total evaluation points to perform both an 8-point and a 16-point integration!

        So why would one not always use Clenshaw-Curtis integration?  The first argument usually given is that with a given number of evaluation points, N, it can exactly integrate a polynomial of degree N - 1.  This is as compared to Gaussian integration which can exactly integrate a polynomial of degree 2N - 1.  I do not find this a very good argument, however, because if one ever had to integrate a polynomial, he would do it by hand and not by numerical integration techniques.  Another reason to avoid Clenshaw-Curtis integration is that every evaluation requires the computation of cos(x) where x varies in equally-spaced intervals from 0 to pi.  The computation time of cos(x) is very long.  One can, however, store the results in a table ahead of time.  Finally, perhaps the best reason to avoid the Clenshaw-Curtis scheme is that it is more complex than the Gaussian quadrature scheme.  It is very easy to write reliable Gaussian integration code, and much more difficult to program Clenshaw-Curtis.  My experience with the Clenshaw-Curtis procedure is that it did not really save any time because of the added complexity.

        With all that said, one can find information about the Clenshaw-Curtis algorithm from the following sources:

 

A list of abscissa and weights are provided here:

        One-dimensional Clenshaw-Curtis integration

 

To Index at Top of Page To Index at Top of Page

Back to Previous Page Back to Numerical Methods

To Home Page To Home Page

 

 

Brought to you by Dr. Chris Bishop, PE.             Contact me

Copyright 2003 All Rights Reserved