Chapter 2
Integration on a one-dimensional standard region

In our finite element formulation we typically require a technique to evaluate, within each elemental domain, integrals of the form

∫ 1
    u(ξ)dξ,
 − 1
(2.1)

where u(ξ) may well be made up of products of polynomial bases. Since the form of u(ξ) is problem specific, we need an automated way to evaluate such integrals. This suggests the use of numerical integration or quadrature. The fundamental building block is the approximation of the integral by a finite summation of the form

∫ 1          q−∑ 1
 − 1u(ξ)dξ ≈    wiu(ξi),
             i=0

where wi are specified constants or weights and ξi represents an abscissa of q distinct points in the interval −1 ≤ ξi ≤ 1. Although there are many different types of numerical integration we shall restrict our attention to Gaussian quadrature.

2.1 Gaussian Quadrature

Gaussian quadrature is a particularly accurate method for treating integrals where the integrand, u(ξ), is smooth. In this technique the integrand is represented as a Lagrange polynomial using the q points ξi, which are to be specified, that is,

      q−1
u(ξ) = ∑ u (ξi)hi(ξ)+ ε(u),
      i=0
(2.2)

where ε(u) is the approximation error. If we substitute equation (2.2) into (2.1) we obtain a representation of the integral as a summation:

∫
  1         q∑−1
 −1u(ξ)dξ =    wiu (ξi)+ R (u),
            i=0
(2.3)

where

          ∫
            1
  wi   =   −1 hi(ξ)dξ,                         (2.4)
          ∫ 1
R (u)  =      ε(u )dξ.                          (2.5)
           −1

Equation (2.4) defines the weights wi in terms of the integral of the Lagrange polynomial but to perform this integration we need to know the location of the abscissae or zeros ξi. Since u(ξ) is represented by a polynomial of order (q − 1) we would expect the relation above to be exact if u(ξ) is a polynomial of order (q − 1) or less [that is, when u(ξ) ∈ Pq−1([−1,1]) then R(u) = 0]. This would be true if, for example, we choose the points so that they are equispaced in the interval. There is, however, a better choice of zeros which permits exact integration of polynomials of higher order than (q − 1). This remarkable fact was first recognised by Gauss and is at the heart of Gaussian quadrature.

We here consider only the result of the Gauss quadrature for integrals of the type shown in equation (2.3) known as Legendre integration. There are three different types of Gauss quadrature known as Gauss, Gauss-Radau, and Gauss-Lobatto, respectively. The difference between the three types of quadrature lies in the choice of the zeros. Gauss quadrature uses zeros which have points that are interior to the interval, −1 < ξi < 1 for i = 0,…,q − 1. In Gauss-Radau the zeros include one of the end-points of the interval, usually ξ = −1, and in Gauss-Lobatto the zeros include both end points of the interval, that is, ξ = 1.

Introducing ξi,P α,β to denote the P zeros of the Pth order Jacobi polynomial PP α,β such that

  α,β  α,β
PP  (ξi,P ) = 0,   i = 0,1,...,P − 1,

where

ξα0,,βP < ξα1,,βP < ⋅⋅⋅ < ξαP,−β1,P,

we can define zeros and weights which approximate the integral

∫ 1         q∑−1
   u (ξ)dξ =    wiu (ξi)+ R (u),
 −1         i=0

as:

(1) Gauss-Legendre

           0,0
   ξi  =  ξi,q      i = 0,...,q − 1

  0,0         2     [ d           ]− 2
 wi    =  --------2  ---(Lq(ξ))|ξ=ξi        i = 0,...,q − 1
          [1− (ξi) ] dξ

R (u)  =  0    if u(ξ) ∈ P2q− 1([− 1,1])

(2) Gauss-Radau-Legendre

          (
          |{  − 1     i = 0
   ξi =   |
          (  ξ0,i−11,q− 1 i = 1,...,q − 1
             (1 − ξ )
 w0i,0 =   -2------i--2     i = 0,...,q − 1
          q [Lq−1(ξi)]

R (u) =   0    if u(ξ) ∈ P2q− 2([− 1,1 ])

(3) Gauss-Lobatto-Legendre

          (|  − 1      i = 0
          ||{
   ξi  =      1,1
          |||(  ξi− 1,q−2  i = 1,...,q − 2
             1        i = q − 1

 w0,i0  =   -------2---------     i = 0,...,q − 1
           q(q − 1)[Lq −1(ξi)]2
R (u)  =  0    if u (ξ) ∈ P2q−3([− 1,1])

In all of the above quadrature formulae Lq(ξ) is the Legendre polynomial (Lq(ξ) = Pq0,0(ξ)). The zeros of the Jacobi polynomial ξi,mα,β do not have an analytic form and commonly the zeros and weights are tabulated. Tabulation of data can lead to copying errors and therefore a better way to evaluate the zeros is by the use of a numerical algorithm (see the appendix in “Spectral/hp element methods for CFD”).