424 10. Orthogonal Polynomials in Approximation Theory

10.3 Chebyshev Integration and Interpolation

If Gaussian quadratures are considered with respect to the Chebyshev

weight w(x) = (1 ’ x2 )’1/2 , Gauss nodes and coe¬cients are given by

(2j + 1)π π

xj = ’ cos , 0 ¤ j ¤ n,

, ±j = (10.20)

2(n + 1) n+1

while Gauss-Lobatto nodes and weights are

πj π

xj = ’ cos , 0 ¤ j ¤ n, n ≥ 1,

, ±j = (10.21)

n dj n

where d0 = dn = 2 and dj = 1 for j = 1, . . . , n ’ 1. Notice that the Gauss

nodes (10.20) are, for a ¬xed n ≥ 0, the zeros of the Chebyshev polynomial

Tn+1 ∈ Pn+1 , while, for n ≥ 1, the internal nodes {¯j , j = 1, . . . , n ’ 1}

x

are the zeros of Tn , as anticipated in Remark 10.2.

Denoting by ΠGL f the polynomial of degree n + 1 that interpolates f

n,w

at the nodes (10.21), it can be shown that the interpolation error can be

bounded as

¤ Cn’s f

f ’ ΠGL f for s ≥ 1,

s,w , (10.22)

w

n,w

where · w is the norm in L2 de¬ned in (10.9), provided that for some

w

s ≥ 1 the function f has derivatives f (k) of order k = 0, . . . , s in L2 . In

w

such a case

1

s 2

f (k) 2

f = . (10.23)

s,w w

k=0

Here and in the following, C is a constant independent of n that can assume

di¬erent values at di¬erent places. In particular, for any continuous function

f the following pointwise error estimate can be derived (see Exercise 3)

f (x) ’ ΠGL f (x) ¤ Cn1/2’s f s,w . (10.24)

∞

n,w

Thus, ΠGL f converges pointwise to f as n ’ ∞, for any f ∈ C 1 ([’1, 1]).

n,w

The same kind of results (10.22) and (10.24) hold if ΠGL f is replaced with

n,w

G

the polynomial Πn f of degree n that interpolates f at the n+1 Gauss nodes

xj in (10.20). (For the proof of these results see, for instance, [CHQZ88],

p. 298, or [QV94], p. 112). We have also the following result (see [Riv74],

p.13)

2

—

f ’ ΠG f ¤ (1 + Λn )En (f ), with Λn ¤ log(n + 1) + 1, (10.25)

∞

n

π

—

where ∀n, En (f ) = inf f ’ p is the best approximation error for f

∞

p∈Pn

in Pn and Λn is the Lebesgue constant associated with the Chebyshev

10.3 Chebyshev Integration and Interpolation 425

nodes (10.20). As far as the numerical integration error is concerned, let

us consider, for instance, the Gauss-Lobatto quadrature rule (10.18) with

nodes and weights given in (10.21). First of all, notice that

1

f (x)(1 ’ x2 )’1/2 dx = lim In,w (f )

GL

n’∞

’1

for any function f whose left integral is ¬nite (see [Sze67], p. 342). If,

moreover, f s,w is ¬nite for some s ≥ 1, we have

1

f (x)(1 ’ x2 )’1/2 dx ’ In,w (f ) ¤ Cn’s f

GL

s,w . (10.26)

’1

This result follows from the more general one

|(f, vn )w ’ (f, vn )n | ¤ Cn’s f ∀vn ∈ Pn ,

vn w, (10.27)

s,w

where the so-called discrete scalar product has been introduced

n

GL

(f, g)n = ±j f (xj )g(xj ) = In,w (f g). (10.28)

j=0

Actually, (10.26) follows from (10.27) setting vn ≡ 1 and noticing that

√

1/2

1 2 ’1/2

vn w = ’1 (1 ’ x ) dx = π. Thanks to (10.26) we can thus

conclude that the (Chebyshev) Gauss-Lobatto formula has degree of ex-

actness 2n ’ 1 and order of accuracy (with respect to n’1 ) equal to s,

provided that f s,w < ∞. Therefore, the order of accuracy is only limited

by the regularity threshold s of the integrand function. Completely similar

considerations can be drawn for (Chebyshev) Gauss formulae with n + 1

nodes.

˜

Let us ¬nally determine the coe¬cients fk , k = 0, . . . , n, of the interpolat-

ing polynomial ΠGL f at the n + 1 Gauss-Lobatto nodes in the expansion

n,w

with respect to the Chebyshev polynomials (10.10)

n

˜

ΠGL f (x) = fk Tk (x). (10.29)

n,w

k=0

Notice that ΠGL f coincides with the discrete truncation of the Chebyshev

n,w

—

series fn de¬ned in (10.4). Enforcing the equality ΠGL f (xj ) = f (xj ), j =

n,w

0, . . . , n, we ¬nd

n

kjπ ˜

f (xj ) = cos fk , j = 0, . . . , n. (10.30)

n

k=0

426 10. Orthogonal Polynomials in Approximation Theory

Recalling the exactness of the Gauss-Lobatto quadrature, it can be checked

that (see Exercise 2)

n

2 1 kjπ

˜

fk = cos f (xj ), k = 0, . . . , n. (10.31)

ndk j=0 dj n

˜

Relation (10.31) yields the discrete coe¬cients {fk , k = 0, . . . , n} in terms

of the nodal values {f (xj ), j = 0, . . . , n}. For this reason it is called

the Chebyshev discrete transform (CDT) and, thanks to its trigonomet-

ric structure, it can be e¬ciently computed using the FFT algorithm (Fast

Fourier transform) with a number of ¬‚oating-point operations of the order

of n log2 n (see Section 10.9.2). Of course, (10.30) is the inverse of the CDT,

and can be computed using the FFT.

10.4 Legendre Integration and Interpolation

As previously noticed, the Legendre weight is w(x) ≡ 1. For n ≥ 0, the

Gauss nodes and the related coe¬cients are given by

2

xj zeros of Ln+1 (x), ±j = , j = 0, . . . , n, (10.32)

(1 ’ x2 )[Ln+1 (xj )]2

j

while the Gauss-Lobatto ones are, for n ≥ 1

x0 = ’1, xn = 1, xj zeros of Ln (x), j = 1, . . . , n ’ 1 (10.33)

2 1

±j = , j = 0, . . . , n (10.34)

n(n + 1) [Ln (xj )]2

where Ln is the n-th Legendre polynomial de¬ned in (10.12). It can be

checked that, for a suitable constant C independent of n,

2 C

¤ ±j ¤ , ∀j = 0, . . . , n

n(n + 1) n

(see [BM92], p. 76). Then, letting ΠGL f be the polynomial of degree n that

n

interpolates f at the n + 1 nodes xj given by (10.33), it can be proved that

it ful¬lls the same error estimates as those reported in (10.22) and (10.24)

in the case of the corresponding Chebyshev polynomial.

Of course, the norm · w must here be replaced by the norm · L2 (’1,1) ,

while f s,w becomes

1

s 2

f (k) 2

f = . (10.35)

L2 (’1,1)

s

k=0

10.4 Legendre Integration and Interpolation 427

The same kinds of results are ensured if ΠGL f is replaced by the polynomial

n

of degree n that interpolates f at the n + 1 nodes xj given by (10.32).

Referring to the discrete scalar product de¬ned in (10.28), but taking

now the nodes and coe¬cients given by (10.33) and (10.34), we see that

(·, ·)n is an approximation of the usual scalar product (·, ·) of L2 (’1, 1).

Actually, the equivalent relation to (10.27) now reads

|(f, vn ) ’ (f, vn )n | ¤ Cn’s f ∀vn ∈ Pn

vn L2 (’1,1) , (10.36)

s

and holds for any s ≥ 1 such that f s < ∞. In particular, setting vn ≡ 1,

√

we get vn = 2, and from (10.36) it follows that

1

f (x)dx ’ In (f ) ¤ Cn’s f

GL

(10.37)

s

’1

which demonstrates a convergence of the Gauss-Legendre-Lobatto quadra-

ture formula to the exact integral of f with order of accuracy s with respect

to n’1 provided that f s < ∞. A similar result holds for the Gauss-

Legendre quadrature formulae.

Example 10.1 Consider the approximate evaluation of the integral of f (x) =

3

|x|±+ 5 over [’1, 1] for ± = 0, 1, 2. Notice that f has “piecewise” derivatives up

to order s = s(±) = ± + 1 in L2 (’1, 1). Figure 10.1 shows the behavior of the

error as a function of n for the Gauss-Legendre quadrature formula. According

to (10.37), the convergence rate of the formula increases by one when ± increases

•

by one.

2

10

0

10

’2

10

’4

10

’6

10

’8

10

’10

10 0 1 2 3

10 10 10 10

FIGURE 10.1. The quadrature error in logarithmic scale as a function of n in the

case of a function with the ¬rst s derivatives in L2 (’1, 1) for s = 1 (solid line),

s = 2 (dashed line), s = 3 (dotted line)

428 10. Orthogonal Polynomials in Approximation Theory

The interpolating polynomial at the nodes (10.33) is given by

n

˜

ΠGL f (x) = fk Lk (x). (10.38)

n

k=0

Notice that also in this case ΠGL f coincides with the discrete truncation

n

—

of the Legendre series fn de¬ned in (10.4). Proceeding as in the previous

section, we get

n

˜

f (xj ) = fk Lk (xj ), j = 0, . . . , n, (10.39)

k=0

and also

± n

2k + 1 1

f (xj ), k = 0, . . . , n ’ 1,

n(n + 1) Lk (xj ) 2

Ln (xj )

j=0

˜

fk = (10.40)

n

1 1

f (xj ), k=n

n+1 Ln (xj )

j=0

(see Exercise 6). Formulae (10.40) and (10.39) provide, respectively, the

discrete Legendre transform (DLT) and its inverse.

10.5 Gaussian Integration over Unbounded

Intervals

We consider integration on both half and on the whole of real axis. In both

cases we use interpolatory Gaussian formulae whose nodes are the zeros of

Laguerre and Hermite orthogonal polynomials, respectively.

The Laguerre polynomials. These are algebraic polynomials, orthogonal

on the interval [0, +∞) with respect to the weight function w(x) = e’x .

They are de¬ned by

dn ’x n

Ln (x) = e n ≥ 0,

x

(e x ),

dxn

and satisfy the following three-term recursive relation

Ln+1 (x) = (2n + 1 ’ x)Ln (x) ’ n2 Ln’1 (x) n ≥ 0,

L’1 = 0, L0 = 1.

∞

For any function f , de¬ne •(x) = f (x)ex . Then, I(f ) = 0 f (x)dx =

∞ ’x

e •(x)dx, so that it su¬ces to apply to this last integral the Gauss-

0

Laguerre quadratures, to get, for n ≥ 1 and f ∈ C 2n ([0, +∞))

n

(n!)2 (2n)

I(f ) = ±k •(xk ) + • (ξ), 0 < ξ < +∞, (10.41)

(2n)!

k=1

10.6 Programs for the Implementation of Gaussian Quadratures 429

where the nodes xk , for k = 1, . . . , n, are the zeros of Ln and the weights

are ±k = (n!)2 xk /[Ln+1 (xk )]2 . From (10.41), one concludes that Gauss-

Laguerre formulae are exact for functions f of the type •e’x , where • ∈

P2n’1 . In a generalized sense, we can then state that they have optimal

degrees of exactness equal to 2n ’ 1.

Example 10.2 Using a Gauss-Laguerre quadrature formula with n = 12 to com-

pute the integral in Example 9.12 we obtain the value 0.5997 with an absolute

error with respect to exact integration equal to 2.96 · 10’4 . For the sake of com-

parison, the composite trapezoidal formula would require 277 nodes to obtain the

•

same accuracy.

The Hermite polynomials. These are orthogonal polynomials on the

2

real line with respect to the weight function w(x) = e’x . They are de¬ned

by

n

n x2 d 2

(e’x ),

Hn (x) = (’1) e n ≥ 0.

dxn

Hermite polynomials can be recursively generated as

Hn+1 (x) = 2xHn (x) ’ 2nHn’1 (x) n ≥ 0,

H’1 = 0, H0 = 1.

∞

2

As in the previous case, letting •(x) = f (x)ex , we have I(f ) = f (x)dx =

’∞

∞ 2

e’x •(x)dx. Applying to this last integral the Gauss-Hermite quadra-

’∞

tures we obtain, for n ≥ 1 and f ∈ C 2n (R)

√

∞ n

(n!) π (2n)

2

e’x •(x)dx = ξ ∈ R,

I(f ) = ±k •(xk ) + n • (ξ),

2 (2n)!

k=1

’∞

(10.42)

where the nodes √ k , for k = 1, . . . , n, are the zeros of Hn and the weights

x

n+1

n! π/[Hn+1 (xk )]2 . As for Gauss-Laguerre quadratures, the

are ±k = 2

2

Gauss-Hermite rules also are exact for functions f of the form •e’x , where

• ∈ P2n’1 ; therefore, they have optimal degrees of exactness equal to 2n’1.

More details on the subject can be found in [DR75], pp. 173-174.

10.6 Programs for the Implementation of Gaussian

Quadratures

Programs 82, 83 and 84 compute the coe¬cients {±k } and {βk }, introduced

in (10.8), in the cases of the Legendre, Laguerre and Hermite polynomials.

These programs are then called by Program 85 for the computation of nodes

430 10. Orthogonal Polynomials in Approximation Theory

and weights (10.32), in the case of the Gauss-Legendre formulae, and by

Programs 86, 87 for computing nodes and weights in the Gauss-Laguerre

and Gauss-Hermite quadrature rules (10.41) and (10.42). All the codings

reported in this section are excerpts from the library ORTHPOL [Gau94].

Program 82 - coe¬‚ege : Coe¬cients of Legendre polynomials

function [a, b] = coe¬‚ege(n)