Jump to ContentJump to Main Navigation
Spectral/hp Element Methods for Computational Fluid Dynamics$

George Karniadakis and Spencer Sherwin

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

(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

( 1 x ) α ( 1 + x ) β , α 1 , β 1
can be numerically evaluated using Gauss-type integration in a discrete summation of the form
(B.0.1)
1 1 ( 1 x ) α ( 1 + x ) β u ( x ) d x = i = 0 Q 1 w i α , β u ( x i ) + R ( u ) ,
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)
u ( x ) = i = 0 Q 1 u ( x i ) h i ( x ) + s ( x ) r ( x ) ,
where
h i ( x j ) = δ i j , h ( x ) P Q 1 , s ( x i ) = 0 , s ( x ) P Q , r ( x ) P Q k ,
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)
w i α , β = 1 1 ( 1 x ) α ( 1 + x ) β h i ( x ) d x ,
(B.0.4)
R ( u = ) 1 1 ( 1 x ) α ( 1 + x ) β s ( x ) r ( x ) d x .
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
R ( u ) = 1 1 ( 1 x ) α ( 1 + x ) β P Q α , β ( x ) r ( x ) d x .
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)
P P α , β ( x i , P α , β ) = 0 , i = 0 , 1, , P 1 ,
where
x 0 , P α , β x 1 , P α , β x P 1 , P α , β .

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:

x i = x i , Q α , β , i = 0 , , Q 1 , w i α , β = H i , Q α , β , i = 0 , , Q 1 , R ( u ) = 0 , u ( x ) P 2 Q 1 ( [ 1 , 1 ] ) ,
(p.596)
H i , Q α , β = 2 α + β + 1 Γ ( α + Q + 1 ) Γ ( β + Q + 1 ) Q ! Γ ( α + β + Q + 1 ) [ 1 ( x i ) 2 ] [ d d x ( P Q α , β ( x ) ) | x = x i ] 2 , R ( u ) 2 α + β + 2 Q + 1 Q ! Γ ( α + Q + 1 ) Γ ( β + Q + 1 ) Γ ( α + β + Q + 1 ) ( 2 Q ) ! ( α + β + 2 Q + 1 ) Γ ( α + β + 2 Q + 1 ) 2 M 2 Q ( M 2 Q = sup 1 x | u 2 Q ( x ) | ) .

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

s ( x ) = ( 1 + x ) s 1 ( x ) ,
where
s 1 ( x ) P Q 1 ( [ 1 , 1 ] ) .
Therefore the relation for R(u) given by eqn (B.0.4) becomes
R ( u ) = 1 1 ( 1 x ) α ( 1 + x ) β + 1 s 1 ( x ) r ( x ) d x .
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:
x i = { 1 , i = 0 , x i 1 , Q 1 α , β + 1 , i = 1 , , Q 1 , w i α , β = { ( β + 1 ) B 0 , Q 1 α , β , i = 0 , B i , Q 1 α , β , i = 0 , , Q 1 , R ( u ) = 0 u ( x ) P 2 Q 2 ( [ 1 , 1 ] ) ,
B i , Q 1 α , β = 2 α + β Γ ( α + Q ) Γ ( β + Q ) ( 1 x i ) ( Q 1 ) ! ( β + Q ) Γ ( α + β + Q + 1 ) [ P ( Q 1 ) α , β ( x i ) ] 2 , R ( u ) 2 α + β + 2 Q ( Q 1 ) ! Γ ( α + Q ) Γ ( β + Q + 1 ) Γ ( α + β + Q + 1 ) ( 2 Q 1 ) ! ( α + β + 2 Q ) [ Γ ( α + β + 2 Q ) ] 2 M 2 Q 1 ( M 2 Q 1 = sup 1 x 1 | u 2 Q 1 ( x ) | ) .

(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

s ( x ) = ( 1 x ) ( 1 + x ) s 2 ( x ) ,
where
s 2 ( x ) P Q 2 ( [ 1 , 1 ] ) ,
and the relation for R(u) given by eqn (B.0.4) becomes
R ( u ) = 1 1 ( 1 x ) α + 1 ( 1 + x ) β + 1 s 2 ( x ) r ( x ) d x .
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:
x i = { 1 , i = 0 , x i 1 , Q 2 α + 1 , β + 1 , i = 1 , , Q 2 , 1 , i = Q 1 , w i α , β = { ( β + 1 ) C 0 , Q 2 α , β , i = 0 , C i , Q 2 α , β , i = 1 , , Q 2 , ( α + 1 ) C Q 1 , Q 2 α , β , i = Q 1 , R ( u ) = 0 , u ( x ) P 2 Q 3 ( [ 1 , 1 ] ) ,
C i , Q 2 α , β = 2 α + β + 1 Γ ( α + Q ) Γ ( β + Q ) ( Q 1 ) ( Q 1 ) ! Γ ( α + β + Q + 1 ) [ P Q 1 α , β ( x i ) 2 ] , R ( u ) 2 α + β + 2 Q 1 ( Q 2 ) ! Γ ( α + Q ) Γ ( β + Q ) Γ ( α + β + Q + 1 ) ( 2 Q 2 ) ! ( α + β + 2 Q 1 ) [ Γ ( α + β + 2 Q 1 ) ] 2 M 2 Q 2 ( M 2 Q 2 = sup 1 x 1 | u 2 Q 2 ( x ) | ) .

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

x i = x i , m α , β , P m α , β ( x i , m α , β ) = 0 , i = 0 , 1 , , m 1 .

(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:

x i , m 1 / 2 , 1 / 2 = cos ( 2 i + 1 2 m π ) , i = 0 , , m 1 ,
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

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

Noting that

f m n ( x ) f m n ( x ) = P m α , β ( x ) [ P m α , β ( x ) ] P m α , β ( x ) i = 0 n 1 1 / ( x x i ) ,
a root-finding algorithm to determine the m roots of P m α , β ( x ) using the Newton–Raphson iteration with polynomial deflation is as follows:
d o k = 0 , m 1 d o r = x k , m 1 / 2 , 1 / 2 if ( k 0 ) r = ( r + x k 1 ) / 2 d o j = 1 , stop s = i = 0 k 1 1 / ( r x i ) δ = P m α , β ( r ) / ( [ P m α , β ( r ) ] P m α , β ( r ) s ) r = r + δ if ( δ ε ) exit loop continue x k = r continue
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).