This document in PDF

Review of Gaussian Quadrature method

Nasser M. Abbasi

June 20, 2014

1 The problem

Figure 1: Integrating a function

To find a numerical value for the integral of a real valued function of a real variable over a specific range over the real line. This means to evaluate

I =    f (x) dx


Geometrically, this integral represents the area under f(x)  from a  to b

2 Solution

Figure 2: Numerical integration

We can always approximate the area by dividing it in equal width strips and then sum the areas of all the strips.

In general, there will always be an error in the estimate of the area using this method. The error will become smaller the more strips we use (which implies a smaller strip width). Hence we can write

∫ b           ( ∑N           )
   f (x ) dx =      Δx  f (x )  + E
 a              i=1

Where E  is the error between the actual area and the approximated area using the above method of numerical integration. N  above is the number of strips or can also be refereed to as the number of integration points.

Instead of keep referring to the 'width of the strip' all the time, we will call this quantity the weight wi  that we will multiply the value of the function with to obtain the area. Hence the above becomes

              (            )
∫b             ∑N
  f (x) dx =       wi f (xi)   + E

Using implied summation on indices the above becomes


  f (x ) dx = wi f (xi) + E

In the above we divided the range of the integration (the distance between the upper and lower limits of integration) into equal intervals. We can decide to evaluate f (x )
    i  at the middle of the strip or at the start of the strip or at the end of the strip. In the diagram above we evaluated the f (x )  at the left end of the strip.

Our goal is to evaluate this integral such as the error E  is minimum and using the smallest number of integration points. In a sense this can be considered an optimization problem with constraints: minimize the error of integration using the smallest possible number of points.

To be able to do this minimization, we need to consider what are the variables involved. We see that there are two degrees of freedom to this problem. One is the width of the strip or wi  . We do not have to use a fixed value of width, we can use different width for different strips if the resulting integral gives better approximation.

The second degree of freedom is the point at which to evaluate f(xi)  associated with each strip. In the example above we choose to evaluate f (xi)  at left end point of the strip. We can choose to select a different xi  point if this will result in a better approximation.

This is the main idea of Gauss Quadrature numerical integration. It is to be able to choose specific values for these two degrees of freedom, the wi  and the xi  .

It turns out that if the function f(x )  is a polynomial, then there is an optimal solution. There is an optimal {wi, xi} for each polynomial of order n  .

We can determine these degrees of freedom such that the error E  is zero, and with the least possible number of integration points. We are able to tabulate these two degrees of freedom for each polynomial of specific order. In other words, if the function f(x)  is a polynomial of order n  then we know before the computation starts what these 2 degrees of freedom should be. We know the locations of xi  and we know weight wi  that we need to multiply f(xi)  with to obtain the area with minumum error.

You might ask how can this method of integration know the locations of the integration points xi  beforehand without being given the integration range of the function to integrate? It turns out that we will map f (x)  into a new known and specific range of integration (from − 1  to + 1  ) for the method that we will now discuss.

2.1 Gauss Quadrature

From now on we will assume the function f(x)  to be integrated is a polynomial in x  of some order n  .

Gauss quadrature is optimal when the function is a polynomial

The main starting point is to represent the function f (x)  as a combination of linearly independent basis.

Instead of using strips of equal width, we assume the width can vary from one strip to the next. Let us call the width of the strip wi  . Instead of taking the height of the  th
i  strip to be the value of the function at the left edge of the strip, let us also be flexible on the location of the x  associated with strip wi  and call the height of the strip wi  as f(xi)  where xi  is to be determined. Hence the above integration becomes


So our goal is to determine wi  and xi  such as the error E  is minimized in the above equation. We would really like to find wi  and xi  such that the error is zero with the smallest value for N  .

It seems as if this is a very hard problem. We have 2N  unknown quantities to determine. N  different widths, and N  associated x  points to evaluate the height of each specific strip at. And we only have as an input f (x )  and the limits of integration, and we need to determine these 2N  quantities such that the error in integration is zero.

In other words, the problem is to find wi,xi  such that


I =   f (x) dx =  w1f (x1) + w2f (x2) + ⋅⋅⋅wN f (xN )

One way to make some progress is to expand f (x )  as a series. We can approximate f(x)  as convergent power series for example. If f(x)  happens to be a polynomial instead, we can represent it exactly using a finite sequence of Legendre polynomials. It is in this second case where this method makes the most sense to use due to the advantages we make from the second representation.

We show both methods below.

Expanding f (x )  as convergent power series over the range a,b  gives

f (x) = a0 + a1x + a2x2 + ⋅⋅⋅amxm  + ⋅⋅⋅

Substituting (2) into (1) gives

    ∫b (                                )
I =     a0 + a1x + a2x2 + ⋅⋅⋅amxm  + ⋅⋅⋅  dx =  w1f (x1) + w2f (x2) + ⋅ ⋅⋅wN f (xN )



∫b (                                )       ∫b       ∫ b         ∫b                ∫b
    a0 + a1x + a2x2 + ⋅⋅⋅amxm  + ⋅⋅⋅  dx =    a0dx +    a1x dx +   a2x2 dx +  ⋅⋅⋅ +  amxm   dx + ⋅⋅⋅

a                                           a         a          a                 a
                           (b2-−-a2)-    (b3-−-a3)           (bm+1-−-am+1-)
          = a0 (b − a ) + a1  2     + a2    3     + ⋅⋅⋅ + am    m  + 1     + ⋅⋅⋅

Substituting the above into (3) results in

                2    2        3    3              m+1    m+1
a0 (b − a) + a1(b-−-a-) + a2(b-−--a-)+  ⋅⋅⋅ + am (b---−-a----) + ⋅⋅⋅
                  2             3                   m +  1
              =  w1f (x1) + w2f (x2) + ⋅ ⋅⋅ + wN f (xN)

But from (2) we see that


Substituting the above into (4) gives

                2    2        3    3              m+1    m+1
a (b − a) + a (b-−--a-)+  a (b-−-a--)+ ⋅⋅⋅ + a  (b----−-a----)+  ⋅⋅⋅ =
 0           1    2(       2    3             m     m)+ 1
                w1  a0 + a1x1 + a2x2+  ⋅⋅⋅amxm  + ⋅⋅⋅
                   (               12         1m       )
               +w2   a0 + a1x2 + a2x 2 + ⋅⋅⋅amx 2 + ⋅⋅⋅
                                 ⋅ ⋅⋅
                   (                                   )
              +wN   a0 + a1xN +  a2x2N + ⋅⋅⋅amxmN +  ⋅⋅⋅

Rearranging gives

              (b2 − a2)     (b3 − a3)           (bm − am )
a0 (b − a) + a1---------+  a2---------+ ⋅⋅ ⋅ + am----------+  ⋅⋅⋅ =
                  2             3                   m
                     a0 (w1 + w2 + ⋅⋅⋅ + wN )
                 +a1 (w1x1 + w2x2 +  ⋅⋅⋅ + wN xN )
                     (    2      2             2)
                 +a2  w1x 1 + w2x 2 + ⋅⋅⋅ + wN x N
                               ⋅ ⋅⋅
                         m       m             m
                +am  (w1x 1 +  w2x2 +  ⋅⋅⋅ + wN xN )

Equating coefficients on both sides results in


Since we have 2N  unknown quantities to solve for (N  weights wi  and N  points xi  ) we need 2N  equations. In other words, we need to have m  = 2N  . The above set of simultaneous 2N  equations can now be solved for the unknown wi,xi  and this will give us the numerical integration we wanted.

The above assumed that the function f(x)  can be represented exactly by the power series expansion with m  terms.

We now consider the representation of f(x)  as a series of Legendre polynomials. We do this since when f (x )  itself is a polynomial. We can represent f (x)  exactly by a finite number of Legendre polynomials. But since Legendre polynomials Pn (x)  are defined over [− 1,1]  we start by transforming f (x)  to this new range and then we can expand the mapped f (x )  (which we will call g(ζ)  ) in terms of the Legendre polynomials.


An easy way to find this mapping is to align the ranges over each others and take the ratio between as the scale factor. This diagram shows this for a general case where we map f (x )  defined over [a, b]  to a new range defined over [c1,c2]


We see from the diagram that

ζ = c1 + dζ


dx  is the same ratio as -b−a-


dx    b − a
---=  -------
dζ    c2 − c1

The above is called the Jacobin of the transformation. Now, From the diagram we see that

dx  = x − a


dζ =  ζ − c1

Hence (6) becomes

x-−-a- = -b −-a-
ζ − c1   c2 − c1

Hence we obtain that

ζ =  x −-a(c  − c ) + c
     b − a  2    1    1


     -b −-a-
x =  c2 − c1 (ζ − c1) + a

For the specific case when c =  − 1
 1  and c  = +1
 2  the above expressions become




Before we proceed further, It will be interesting to see the effect of this transformation on the shape of some functions. Below I plotted some functions under this transformation. The left plots are the original functions plotted over some range, in this case [4,10 ]  and the right side plots show the new shape (the function g (ζ )  ) over the new range [− 1,1]


We must remember that in the following analysis, we are integrating now the function g (ζ)  over [− 1, 1]  and not f (x )  over [a,b]  . Hence to obtain the required integral we need to transform back the final result. We will show how to do this below.

We can approximate any function f (x)  as a series expansion in terms of weighted sums of a set of basis functions. We do this for example when we use Fourier series expansion.

Hence we write

f (x) =     a Φ (x)
             i i

We can give an intuitive justification to the above formulation as follows. If we think in terms of vectors. A vector V  in an n-dimensional space is written in terms of its components as follows


Where ei  is the basis vectors in that space.

If we now consider a general infinite dimensional vector space where each point in that space is a function, then we see that we can also represent that function as a weighted series of a basis functions just as we did for a normal vector.

There are many sets of basis functions we can choose to represent the function f (x )  with. The main requirements for the basis functions is that they are complete (This means they span the whole space) and there is defined an inner product on them.

For our purposes here, we are interested in the class of function f (x)  that are polynomials in x  . The basis we will select to use are the Legendre basis as explained above. To do this, we transform f (x)  to g (ζ)  as shown above and then now our integral becomes

∫ b          ∫ 1        (          )
   f (x )dx =    f (x (ζ))  (b-−-a)dζ
 a           −1

This is because we found that dx = (b−a)dζ
       2  from above.

If we call f (x(ζ ))  as g (ζ)  the integral becomes

∫b           ∫1
  f (x)dx =    (b-−-a)g (ζ) dζ
a           − 1

Since (b−a)-
 2   is the Jacobian of the transformation, we write the integral as

∫b           ∫ 1
  f (x)dx →     J g(ζ) d ζ

a            −1

Since the Jacobian is constant in this case, we will only worry about  1
   g (ζ) dζ

−1  and we just need to remember to scale the result at the end by J  . This is the integral we will now numerically integrate.

Equation (7) is now written as

g(ζ) =     aiPi (ζ)

Where Pi(ζ)  is the Legendre polynomial of order i  and g(ζ )  is a polynomial in ζ  or some order m  .

Now we carry the same analysis we did when we expressed f (x)  as a power series. The difference now is that the limit of integration is symmetrical and the basis are now the Legendre polynomials instead of the polynomials xn  as the case was in the power series expansion. So now equation (1) becomes

    ∫ 1

I =    g (ζ) dζ = w1g (ζ1) + w2g (ζ2) + ⋅⋅⋅ + wN g(ζN )

And equation (2) becomes

g(ζ) = a0P0 (ζ) + a1P1(ζ) + a2P2 (ζ) + ⋅⋅⋅amPm  (ζ ) + ⋅⋅⋅

Substituting (9) into (8) we get the equivalent of equation (3)


The above occurs since the integral of any Legendre polynomial of order greater than zero will be zero over [− 1,1]

Substituting the above into (10) we obtain

I = 2a0 = w1g (ζ1) + w2g (ζ2) + ⋅⋅⋅ + wN g(ζN )

But from (9) we see that


Substituting the above in (11) gives


Rearranging results in


Equating coefficients gives


If we select the points ζi  to be the roots of Pi− 1   we can write the above as


3 References

Mathematica Structural Mechanics help page
MIT course 16.20 lecture notes. MIT open course website.
Theory of elasticity by S. P. Timoshenko and J. N. Goodier. chapter 10