Home > Reference > Numerical Methods > Numerical Integration
Index of this pageGaussian Quadrature Integration
Gaussian Quadrature IntegrationIn 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 outofprint 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 GaussLegendre 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 GaussLegendre 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 onedimensional argument (x), it also applies to multidimensional 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 curvefit 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 highdegree polynomial to allow the freedom to change slope rapidly. Highdegree polynomials, however, tend to oscillate wildly around the discontinuity. Lowdegree polynomials can not change slope fast enough to model the discontinuity. The result is that no polynomial, lowdegree or highdegree, can sufficiently model a function with a step or slope discontinuity. Therefore, no Gaussian quadrature can produce an accurate, reliable result.
LinearDomain GaussLegendre QuadraturesThe following onedimensional 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.
LinearDomain Gauss Quadrature weights and abcissa
TriangularDomain GaussLegendre QuadraturesAt 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 11291148, 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 zeta3 for the 7point 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 nondegenerate straightsided 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.
TriangularDomain Gauss Quadrature weights and abscissa
TetrahedralDomain GaussLegendre QuadraturesFinally, there are available formulas for integration over tetrahedra. See, for example, P. Keast, "ModerateDegree 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.
TetrahedralDomain Gauss Quadrature weights and abscissa
Quadratures Over Circles and SpheresI 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 equallyspaced 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.
ClenshawCurtis IntegrationAlthough fewer people use ClenshawCurtis 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, ClenshawCurtis 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 ClenshawCurtis 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 stillsuspicious 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 8point and a 16point integration! So why would one not always use ClenshawCurtis 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 ClenshawCurtis integration is that every evaluation requires the computation of cos(x) where x varies in equallyspaced 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 ClenshawCurtis 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 ClenshawCurtis. My experience with the ClenshawCurtis 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 ClenshawCurtis algorithm from the following sources:
A list of abscissa and weights are provided here: Onedimensional ClenshawCurtis integration

Brought to you by Dr. Chris Bishop, PE. Contact me
Copyright 2003 All Rights Reserved