Print publication date: 2005

Print ISBN-13: 9780198528692

Published to Oxford Scholarship Online: September 2007

DOI: 10.1093/acprof:oso/9780198528692.001.0001

Show Summary Details
Page of

PRINTED FROM OXFORD SCHOLARSHIP ONLINE (www.oxfordscholarship.com). (c) Copyright Oxford University Press, 2017. All Rights Reserved. Under the terms of the licence agreement, an individual user may print out a PDF of a single chapter of a monograph in OSO for personal use (for details see http://www.oxfordscholarship.com/page/privacy-policy). Subscriber: null; date: 25 February 2017

(p.594) APPENDIX B GAUSS-TYPE INTEGRATION

Source:
Spectral/hp Element Methods for Computational Fluid Dynamics
Publisher:
Oxford University Press

Integration of u(x) in the interval [−1, 1] with respect to the function

$Display mathematics$
can be numerically evaluated using Gauss-type integration in a discrete summation of the form
(B.0.1)
$Display mathematics$
where R(u) = 0 if u(x) ∈ P 2Qk([−1, 1]). The value of k is determined by the type of quadrature used, which can be either classical Gauss (k = 1), Gauss–Radau (k = 2), or Gauss–Lobatto (k = 3).

To derive eqn (B.0.1) we consider the case where u(x) is a polynomial of order 2Qk (that is, u(x) ∈ P 2Qk([−1, 1])) and decompose u(x) into

(B.0.2)
$Display mathematics$
where
$Display mathematics$
and δ‎ij denotes the Kronecker delta. The function u(x) has been decomposed into h i(x j) which is Lagrangian polynomial through the Q nodal points x i, s(x) which is a polynomial of order Q with roots at the nodal points x i, and the remainder polynomial r(x) of order Qk. Integrating eqn (B.0.2) with respect to (1 − x)α‎ (1 + x)β‎ in the interval [−1, 1] and comparing with eqn (B.0.1), we deduce
(B.0.3)
$Display mathematics$
(B.0.4)
$Display mathematics$
If s(x) is specified then the nodal points x i will be determined as the roots of s(x) such that s(x i) = 0. Specifying s (x) is preferable to specifying the nodal (p.595) points x i because an appropriate choice of s(x) will make R(u) = 0. For example, if $s ( x ) = P Q α , β ( x )$, where $P Q α , β ( x )$ is a Jacobi polynomial (see Appendix A), then
$Display mathematics$
This integral is the orthogonality relationship for $P Q α , β ( x )$ (see eqn (A.1.7), so, if k = 1, then r(x) is a polynomial of order Q − 1 and therefore R(u) = 0. This choice of s(x) determines the classical Gauss quadrature and therefore the nodal points are determined by $x i = x i , P α , β$, where $x i , P α , β$ denote the m zeros of the Jacobi polynomial $P P α , β$ such that
(B.0.5)
$Display mathematics$
where
$Display mathematics$

Gauss–Radau and Gauss–Lobatto integration both require that the nodal points include one or both of the end-points x = ±1. This restriction modifies the form of s(x) and has the effect of altering the weight function (1 − x)α‎ (1+x)β‎. Therefore the nodal points are determined by a Jacobi polynomial with different values of α‎ and β‎.

B.1 Jacobi formulae

In this section we state the nodal values x i and weights $w i α , β$ for Gaussian integration with respect to the function (1 − x)α‎ (1 + x)β‎, where α‎, β‎ 〉 −1. The Legendre case, when α‎ = β‎ = 0, was introduced in Section 2.4.1 and can be considered as a special case of the following formulae. The nodal points x i are determined in terms of the roots of Jacobi polynomials $x i , m α , β$ defined by eqn (B.0.5). A complete derivation of the weight formulae can be found in Ghizzetti and Ossicini [182].

Gauss–Jacobi

For this type of quadrature no restrictions are placed on the nodal points and the nodal points are the roots of $P Q α , β ( x )$. The Gauss–Jacobi formulae are as follows:

$Display mathematics$
(p.596)
$Display mathematics$

Gauss–Radau–Jacobi—Bouzitat formulae of the first kind

Gauss–Radau-type integration requires that the nodal points include one of the end-points of the integration interval x = ±1. If we choose to include x = −1 then the value of s(x) in eqn (B.0.2) has the form

$Display mathematics$
where
$Display mathematics$
Therefore the relation for R(u) given by eqn (B.0.4) becomes
$Display mathematics$
If we let $s 1 ( x ) = P Q − 1 α , β + 1 ( x )$ then R(u) = 0 if k = 2. The Gauss–Radau–Jacobi formulae are as follows:
$Display mathematics$
$Display mathematics$

(p.597) Gauss–Lobatto–Jacobi—Bouzitat formulae of the second kind

Gauss–Lobatto-type integration requires that the nodal points include both of the end-points of the integration interval x = ±1. For this case the value of s(x) in eqn (B.0.2) has the form

$Display mathematics$
where
$Display mathematics$
and the relation for R(u) given by eqn (B.0.4) becomes
$Display mathematics$
Therefore, if we let $s 2 ( x ) = P Q − 2 α + 1 , β + 1 ( x )$ then R(u) = 0 if k = 3. The Gauss–Lobatto–Jacobi formulae are as follows:
$Display mathematics$
$Display mathematics$

B.2 Evaluation of the zeros of Jacobi polynomials

The formulae for the weights in Section B.1 have a closed form in terms of the nodal points x i. In general, however, there are no explicit formulae for the nodes. These are defined in terms of the roots of the Jacobi polynomial such that

$Display mathematics$

(p.598) The zeros $x i , m α , β$ can be numerically evaluated using an iterative technique such as Newton–Raphson. However, we note that the zeros of the Chebyshev polynomial $( α = β = − 1 2 )$ do have an explicit form:

$Display mathematics$
and so we can use $x i , m − 1 / 2 , − 1 / 2$ as an initial guess to the iteration.

To ensure that we find a new root at each search we can apply polynomial deflation, where the known roots are factored out of the initial polynomial once they have been determined. This means that the root-finding algorithm is applied to the polynomial

$Display mathematics$
where x i (i = 0,…, n − 1) are the known roots of $P m α , β ( x )$.

Noting that

$Display mathematics$
a root-finding algorithm to determine the m roots of $P m α , β ( x )$ using the Newton–Raphson iteration with polynomial deflation is as follows:
$Display mathematics$
Here ε‎ is a specified tolerance. Numerically, we find that a better approximation for the initial guess is given by the average of $r = x k , m − 1 / 2 , − 1 / 2$ and $x k − 1$. The values of $P m α , β ( x )$ and $[ P m α , β ( x ) ] ′$ can be generated using the recursion relationships (A.1.3) and (A.1.4).