The weights of classical Gauss–Radau quadratures with double end-point

The weights of classical Gauss–Radau quadratures with double end-point

Applied Numerical Mathematics 60 (2010) 574–586 Contents lists available at ScienceDirect Applied Numerical Mathematics www.elsevier.com/locate/apnu...

436KB Sizes 4 Downloads 11 Views

Applied Numerical Mathematics 60 (2010) 574–586

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

The weights of classical Gauss–Radau quadratures with double end-point B.D. Welfert School of Mathematical & Statistical Sciences, Arizona State University, Tempe, AZ 85287-1804, United States

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 18 January 2008 Received in revised form 10 March 2010 Accepted 12 March 2010 Available online 17 March 2010

Explicit expressions for the weights of classical Gauss–Radau quadratures with double end-point are derived. Applications to Jacobi and generalized Laguerre polynomials are considered. Asymptotic estimates of interior weights are also suggested and verified numerically. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.

Keywords: Numerical quadrature Sturm–Liouville problem Orthogonal polynomials

1. Introduction Generalized Gauss–Radau quadrature formulae

b

f (x) dλ(x) ≈ w 0 f (a) + w 0 f  (a) +

N 

w j f (x j ),

(1)

j =1

a

where dλ(x) is a positive measure on [a, b] are of practical interest in the solution of boundary value problems via spectral methods as well as in the determination of moment-preserving spline approximations [6]. It is well known that there exist a unique set of nodes {x1 , . . . , x N } and weights { w 0 , w 0 , w 1 , . . . , w N } such that (1) is exact whenever f is in P2N +1 , the set of polynomials of degree 2N + 1 or less. The nodes are the zeros of the (normalized) polynomial R N of degree N which is orthogonal to P N −1 with respect to the modified measure (x − a)2 dλ(x), while the weights can be expressed as integrals of interpolating functions. Both can be obtained numerically by solving a tridiagonal matrix eigenvalue problem. A computational procedure for general continuous or discrete measures dλ(x) is described in [2]. A Matlab implementation (gradau.m) is available [5]. A more stable version (gradau_jacobi.m) is also available for the Jacobi weight function [8]. We focus here on classical quadrature formulae of the form (1), where dλ(x) = w (x) dx is such that the Sturm–Liouville problem

  − q(x) w (x)u  = λ w (x)u ,

x ∈ [a, b],

(2)

with

q(x) = q (a)(x − a) +

q 2

(x − a)2

E-mail address: [email protected] URL: http://math.asu.edu/~bdw. 0168-9274/$30.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2010.03.005

(3)

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

575

admits polynomial solutions. Note that

q(a) = 0.

(4)

There are essentially two cases: 1. [Jacobi] q(x) = 1 − x2 and w (x) = (1 + x)α (1 − x)β with α , β > −1 in (a, b) = (−1, 1). 2. [Generalized Laguerre] q(x) = x and w (x) = xα e −x with α > −1 in (a, b) = (0, ∞). Consideration of the more general case (3) will enable us to directly apply results obtained in [13] using a similar machinery, and obtain common expressions for the quadrature weights for both cases. The more unified approach can also possibly serve as a template for dealing with quadrature rules based on weight functions arising in generalized Sturm–Liouville problems with infinite sequences of orthogonal polynomial solutions [11,10]. Numerical tests conducted in [3] (on Gauss–Lobatto quadrature rules using f (b) rather than f  (a) in (1)) indicate that explicit expressions of the quadrature weights associated to end-points tend to produce numerically more accurate values than those obtained via the eigenvalue problem. These weights are typically determined by substituting specific polynomials into the quadrature rule and solving a linear system. However the conditioning of this system worsens with N, to the point where an explicit solution may become necessary. Explicit expressions for boundary weights in the case of the Radau quadrature (1) are only known in a few cases: for the Legendre weight function w (x) = 1 [6, Examples 3.10 and 3.13] and the four Chebyshev weight functions w (x) = 1

1

(1 − x)± 2 (1 + x)± 2 [1]. Other than that, little is known.

In this work we derive explicit formulae for all weights of the quadrature. While the boundary weight w 0 and interior weights w j , j = 1, . . . , N, can be easily obtained from general formulae for Radau weights [13,4], the derivation of the boundary weight w 0 on the other hand requires a more subtle approach, described in Section 2.5, based on the judicious choice of a polynomial function f for which the right-hand side simply vanishes. Applications to the Jacobi and generalized Laguerre polynomials are considered. In addition asymptotic estimates of the weights are given. 2. Quadrature weights Let f (x) = f (a) + (x − a) g (x) in (1). Then

b

b f (x) dλ(x) = f (a)

a

b g (x) dλr (x)

dλ(x) + a

a

with dλ (x) ≡ (x − a) dλ(x). The quadrature rule (1) can be obtained by applying an ( N + 1)-point Gauss–Radau quadrature rule of order 2N to the second integral: r

b

b dλ(x) + w r0 g (a) +

f (x) dλ(x) ≈ f (a) a

 b =

N 

w rj g (x j )

j =1

a

dλ(x) −

N  w rj j =1

a

 f (a) + w r0 f  (a) +

xj − a

N  w rj j =1

xj − a

f (x j ).

Alternately, letting f (x) = f (a) + (x − a) f  (a) + (x − a)2 h(x) in (1) and applying an N-point Gauss rule to the integral involving h yields

b

b f (x) dλ(x) ≈ f (a)

a

b =



b

dλ(x) + f (a) a

dλ(x) −

j =1

a

w rrj h(x j )

j =1

a N 

N 

(x − a) dλ(x) + b



w rrj

(x j − a)2

f (a) +

(x − a) dλ(x) −

N  w rrj j =1

a

xj − a

 f  (a) +

N  j =1

w rrj

(x j − a)2

f (x j ).

As a result the nodes {x1 , . . . , x N } in (1) are the same nodes arising from the Gauss–Radau quadrature applied to g w.r.t. the modified measure dλr (x) = (x − a) dλ(x) and from the Gauss quadrature applied to h w.r.t. the modified measure dλrr = (x − a)2 dλ. These nodes are the zeros of the orthogonal polynomial of degree N w.r.t. dλrr . Also,

b w0 =

dλ(x) − a

N  w rj j =1

xj − a

b =

dλ(x) − a

N  j =1

w rrj

(x j − a)2

,

(5a)

576

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

b

w 0 = w r0 =

dλr (x) −

N  w rrj j =1

a

wj =

w rj xj − a

=

w rrj

(x j − a)2

xj − a

(5b)

,

(5c)

.

For classical measures dλ(x) (5b), (5c) can be used to obtain explicit formulae for w 0 and w j from known expressions for Gauss–Radau weights w rj [13,4]. On the other hand w 0 =

b a

dλ(x) −

N

j =1

w j (obtained from (5) or by substituting

f (x) = 1 into (1)) provides a simple, but unstable, way of computing w 0 , due to cancellations (w 0  different approach for determining w 0 shall be used instead.

b a

dλ(x)) and a

2.1. Preliminaries In the following we shall assume that the measure dλ(x) is classical, i.e.: 1. The measure dλ(x) = w (x)dx with w positive on (a, b). 2. The Sturm–Liouville eigenvalue problem (2) with q of the form (3) such that (qw )(a) = (qw )(b) = 0 admits a polynomial solution P n of degree n. In particular,

q(x) P n + (x) P n + λn P n = 0 with





(x) ≡ q (x) + q

w

(6)



w

(x).

(7)

3. The function  is a linear polynomial. Then

λn = −

n(n − 1)  q − n  2

(8)

(e.g. differentiate (6) n times). 4. For convenience P n is normalized such that

P n (a) = 1.

(9)

In particular P 0 = 1. The normalization (9) is nonstandard but facilitates the derivation of an explicit expression for the weight w 0 . The sequence { P 0 , P 1 , . . .} forms an orthogonal set w.r.t. to the inner product

b ( f , g) ≡

f (x) g (x) w (x) dx,

(10)

a

i.e.,

(Pn, Pk) =

0

if k = n,

(11)

γn if k = n.

If P −1 ≡ 0 the sequence { P 1 , P 2 , . . .} can be generated from a three-term recurrence relation

xP n = an,n+1 P n+1 + an,n P n + an,n−1 P n−1 ,

n > 0.

(12)

The normalization (9) then implies that

(x − a) P n = an,n+1 ( P n+1 − P n ) − an,n−1 ( P n − P n−1 ).

(13)

2.2. Induced polynomials Define the induced polynomial

αn Q n ≡

P n +1 − P n x−a

of degree n, where

αn = 0 is a scaling factor such that Q n (a) = 1. Then:

(14)

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

577

1. The relation (13) yields

P n = an,n+1 αn Q n − an,n−1 αn−1 Q n−1 .

(15)

Note that the scaling of P n and Q n implies

1 = an,n+1 αn − an,n−1 αn−1 .

(16)

2. Q n is orthogonal to Pn−1 w.r.t. the inner product

b f (x) g (x) w r (x) dx,

( f , g )r ≡

(17)

a

with

w r (x) = (x − a) w (x). Indeed,





αn ( Q n , p )r = αn (x − a) Q n , p = ( P n+1 − P n , p ) = 0 for all polynomials p of degree less than n. Therefore



( Q n , Q k )r = The quantity

if k = n,

0

(18)

γnr if k = n.

γnr can be expressed in terms of γn as follows: 



γnr = ( Q n , Q n )r = Q n , (x − a) Q n =

( Q n , P n +1 − P n )

αn ( P n + an,n−1 αn−1 Q n−1 , P n ) γn =− =− . an,n+1 αn2 an,n+1 αn2

=−

( Q n, Pn)

αn (19)

Remark 1. (19) implies an,n+1 < 0. 3. Q n satisfies a Sturm–Liouville problem

  − q(x) w r (x) Q n = λnr w r (x) Q n ,

(20)

i.e.,

qr (x) Q n + r (x) Q n + λnr Q n = 0, with

qr (x) = q(x),



r (x) = q (x) + q

( w r )

(x) = (x) +

wr

= (x) + q (a) +

(21)



q 2

q(x) x−a

(x − a)

(22)

(r is again a linear polynomial). In particular

r (a) = (a) + q (a)

(23)

and

λnr = −

  q λn + λn+1 +  n(n − 1)  r  q − n r = λn − n = . 2 2 2

(24)

4. The sequence { Q 0 , Q 1 , . . .} satisfies a three-term recurrence relation

xQ n = anr ,n+1 Q n+1 + anr ,n Q n + anr ,n−1 Q n−1 . The coefficient anr ,n+1 can be related to an+1,n+2 as follows (LC ( p ) denotes the leading coefficient of a polynomial p):

anr ,n+1 =

LC (xQ n ) LC ( Q n+1 )

=

LC ((x − a) Q n ) LC ((x − a) Q n+1 )

=

αn+1 LC ( P n+1 ) αn+1 LC (xP n+1 ) αn+1 = = an+1,n+2 . αn LC ( P n+2 ) αn LC ( P n+2 ) αn

(25)

578

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

5. The polynomial

Rn ≡

Q n +1 − Q n

(26)

x−a

has degree n and defines a sequence of polynomials induced by { Q 0 , Q 1 , . . .} which are orthogonal w.r.t. the inner product

b f (x) g (x) w rr (x) dx

( f , g )rr ≡

(27)

a

with

w rr (x) = (x − a) w r (x) = (x − a)2 w (x). Indeed,





αn+1 αn ( R n , p )rr = αn ( P n+2 − P n+1 ) − αn+1 ( P n+1 − P n ), p = 0 for all polynomials p of degree less than n. 6. The quadrature nodes {x1 , . . . , x N } are the zeros of R N , i.e., the zeros of Q N +1 − Q N . From (15) and (16) we obtain

P N +1 (x j ) = Q N (x j ) = Q N +1 (x j ).

(28)

2.3. The boundary weight w 0 According to (5) the boundary weight w 0 can be obtained from the weight w r0 of an ( N + 1)-point Gauss–Radau quadrature of order 2N. An explicit expression for w r0 was derived in [13]:

w r0 = r (a) where arN , N +1 , (29) yields

w 0 =

γNr arN , N +1 (λrN − λrN +1 )

(29)

,

γNr , r and λrN are quantities associated to the sequence { Q n }n0 . Substituting (25), (19), (24) and (23) into 2[(a) + q (a)]γ N

a N , N +1 a N +1, N +2 α N α N +1 (λ N +2 − λ N )

.

The expression (A.2) in Appendix A yields the alternate form

w 0 =

2[(a)]2 [(a) + q (a)]γ N a N , N +1 a N +1, N +2 (λ N +2 − λ N +1 )(λ N +1 − λ N )(λ N +2 − λ N )

.

(30)

2.4. The interior weights w j , j = 1, . . . , N The internal weights { w 1 , . . . , w N } can be determined from the weights w rj of the Gauss–Radau quadrature as well. The relation [13, expression (23)]

w rj =

w r0

qr ( x j )

(31)

r (a)(x j − a) [ Q N (x j )]2

yields, by virtue of (21), (23), and (5b), (5c),

wj =

w 0

q (x j )

[(a) + q (a)](x

j

− a)2

[ Q N (x j )]2

(32)

,

which can easily be generated from w 0 and the nodes x j . From (28) an alternate expression is

wj =

w 0

q (x j )

[(a) + q (a)](x

j

− a)2

[ P N +1 (x j )]2

.

(33)

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

579

2.5. The boundary weight w 0 The weight w 0 is the most difficult one to evaluate. Rather than using interpolation we apply the quadrature (1) to the polynomial f = P N +1 Q N of degree 2N + 1. The left-hand side vanishes by orthogonality (11), while the right-hand side can be expressed as a linear combination of w 0 and w 0 thanks to (32). We obtain, using P N +1 (a) = Q N (a) = 1, (32), and (28),

 





0 = w 0 + w 0 P N +1 (a) + Q N (a) +

N 



q (x j )

[(a) + q (a)](x j − a)2

j =1

,

i.e., using (3),

w0 w 0

 



= − P N +1 (a) + Q N (a) +

N 

 



= − P N +1 (a) + Q N (a) +

j =1



q (x j )

[(a) + q (a)](x j − a)2 

1

(a) + q (a)

Nq 2



− q (a)

N  j =1

1



a − xj

.

(34)

The quantity P N +1 (a) can be evaluated directly using the Sturm–Liouville problem (2), see (A.1) in Appendix A. The term Q N (a) can be evaluated using the Sturm–Liouville problem (20) or directly using the definition (14), see (A.3) in Appendix A. Finally, the sum appearing in (34) can be recognized as the value at x = a of the logarithmic derivative of RN: N  j =1



1

=

a − xj

=

N 

 ln |x − x j |

j =1

x=a

  = ln |(x − x1 ) · · · (x − xN )| x=a

[(x − x1 ) · · · (x − xN )] (x − x1 ) · · · (x − xN )

x=a

=

R N (a) R N (a)

,

which can be evaluated from the definition (26), see (A.4) in Appendix A. Substitution into (34) yields an expression which simplifies to (after a tedious but straightforward calculation)



 (2 + q )((a) + q (a))  q ((a) + q (a)) 2 . = − N + N + w 0 (a)((a) + 2q (a)) (a)((a) + 2q (a)) (a) w0

The positivity of all weights w i , i = 1, . . . , N, w 0 , and w 0 was proved in [7]. 3. Applications In this section

z! = Γ ( z + 1), where Γ is the standard Gamma function, and

z

y

=

z! y !( z − y )!

denotes the binomial coefficient. 3.1. Gauss–Radau–Jacobi quadrature with double endpoint (α ,β)

The Jacobi polynomials P n = P n

w (x) = w

(α ,β)

are orthogonal w.r.t. the weight function

(x) = (1 − x)α (1 + x)β

in the interval (a, b) = (−1, 1) and satisfy (6) with

q(x) = 1 − x2 ,

(x) = β − α − x(α + β + 2),

λn = n(n + α + β + 1).

(35)

580

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

Table 1 Boundary weights obtained from (36) and (37) for selected Gauss–Radau–Jacobi rules. The case (α , β) = (1, 0) is used in Section 3.1.2.

α

1/ w 0

β

w 0 / w 0

1

0

( N +1)( N +2) ( N +3)

0

0

( N +1)2 ( N +2)2

− 12

− 12

1 2

2

Reference

2 3 2 3

N2 +

( N +1)(2N +1)(2N +3)

6 5

N +

1 2

( N +1)( N +2)( N +3)(2N +3)(2N +5)

10 21

N2 +

40 21

N +1

− 12

1 2

( N +1)( N +2)(2N +1)(2N +3)(2N +5)

10 21

N +

10 7

2 3

1 2

− 12

( N +1)( N +2)(2N +3)

6 5

16

8



45π

45π



(α ,β)

When normalized by P n

an,n+1 = −

8 3

N+

3 2

N 2 + 2N + 1 2

12 5

2

N2 +

18 5

[6, p. 164]

N +1

N+

[1, (2.18-19)]

N +2

[1, (2.38-39)] [1, (2.41-42)] [1, (2.44-45)]

(−1) = 1 these polynomials satisfy the recurrence (12) with

2(n + 1 + β)(n + 1 + β + α ) , (2n + 1 + α + β)(2n + 2 + α + β)

α2 − β 2

an,n = −

, (2n + α + β)(2n + α + β + 2) 2n(n + α ) an,n−1 = − (2n + α + β)(2n + 1 + α + β) for n > 0. Furthermore

1

γn =

(α ,β) Pn

2 an,n−1 (x) w (α ,β) (x) dx = γn−1 = · · · = an−1,n

−1

=

2α +β+1 (β!)2 n!(n + α )!

(n + β)!(n + α + β)!(2n + 1 + α + β)

β!(α + β + 1)!n!(n + α )! γ0 α !(n + β)!(n + α + β)!(2n + 1 + α + β)

.

3.1.1. Explicit weights Then (30) yields

1 w 0

=

2+β



22+α +β

N +2+β



N +2+β +α



N +α

N

(36)

. (α ,β+1)

The resulting expression of w 0 matches the Gauss–Radau weight λ0 obtained from [4, formula (3.10)], in accordance with (5b). For selected values of (α , β) (in particular integers/half-integers) the formula (36) can be simplified. Table 1 √ contains the Legendre and Chebyshev cases (the latter uses the value (−1/2)! = Γ (1/2) = π ).  The ratio w 0 / w 0 from (35) becomes

w0 w 0

=

β +2 α+β +2 . N ( N + α + β + 3) + (β + 1)(β + 3) 2(β + 1)

(37)

This expression matches available results from the literature, listed in Table 1. Because of orthogonality properties and normalization at a = −1 the relations (α ,β+1)

Q n = Pn

,

(α ,β+2)

R n = R n (−1) P n

(38)

hold. The interior weights (32) therefore simplify to

wj =

(1 − x j ) w 0 (α ,β+1)

2(β + 2)(1 + x j )[ P N

(x j )]2

(39)

. (α ,β)

From (33) an equivalent expression in terms of P N +1 is

wj =

(1 − x j ) w 0 (α ,β)

2(β + 2)(1 + x j )[ P N +1 (x j )]2

.

(40)

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

581

3.1.2. Numerical illustration As an illustration we carry out a numerical comparison of three approaches on the integral

1 −1

with

1−x

(1 +

x)2



2

dx = (2/ω) arctan(2/ω) −

1 2



ln 1 + (2/ω)2



(41)

ω = 0.001, α = 1, β = 0, f (x) = 1/((1 + x)2 + ω2 ) ( f  (−1) = 0):

1. Using gauss.m [5] to define interior nodes x j and weights w rrj associated to the weight function w (1,2) satisfying

dλrr = (1 + x)2 dλ = (1 − x)(1 + x)2 dx = w (1,2) (x) dx and using (5c), together with the explicit expressions (36) and (37) for w 0 and w 0 :

ab = r_jacobi(N,1,2); xw = gauss(N,ab); wp0 = 16/((N+1)*(N+2)^2*(N+3)); w0 = wp0*(2*N*(N+4)/3+3/2); w = [w0;xw(:,2)./(1+xw(:,1)).^2]; x = [-1;xw(:,1)]; f = 1./((1+x).^2+om^2); int = sum(w.*f); 2. Using gradau.m [5], or the modified version gradau_jacobi.m [8] to handle the double quadrature point a = −1:

%-------------------------- either ------------------------ab = r_jacobi(N+2,1,0); xw = gradau(N,ab,2,-1); %---------------------------- or --------------------------xw = gradau_jacobi(N,1,0,2); %----------------------------------------------------------i = [1 3:N+2]; x = xw(i,1); w = xw(i,2); f = 1./((1+x).^2+om^2); int = sum(w.*f); 3. Using the explicit expressions (36), (37) for w 0 , w 0 , together with either (39) or (40) for the interior weights w j (using (1,1)

the recurrence relation (12) to evaluate either P N

(1,0)

or P N +1 ; the nodes x j are obtained by gauss.m):

ab = r_jacobi(N,1,2); xw = gauss(N,ab); x = xw(:,1); %-------------------------- either ------------------------P0 = 1; P1 = -x; for n = 1:N-1 P2 = -((2*n+3)*x.*P1+n*P0)/(n+3); P0 = P1; P1 = P2; end %P2 = P_N^(1,1)(x) %---------------------------- or --------------------------P0 = 1; P1 = -(3*x+1)/2; for n = 1:N nn = 2*n+1; P2 = -((1/nn+(nn+2)*x).*P1+(n*(nn+2)/nn)*P0)/(n+2); P0 = P1; P1 = P2; end %P2 = P_(N+1)^(1,0)(x) %----------------------------------------------------------wp0 = 16/((N+1)*(N+2)^2*(N+3)); w0 = wp0*(2*N*(N+4)/3+3/2); w = [w0; (1-x).*wp0./(4*(1+x).*P2.^2)]; x = [-1;x]; f = 1./((1+x).^2+om^2); int = sum(w.*f); Fig. 1 displays the relative accuracy obtained using the various strategies. All methods yield spectrally convergent approximations. The breakdown of gradau.m for N  400, first reported by the author, was later investigated by Gautschi who determined that improper scaling of the orthogonal polynomials { P n }n0 was at fault [9]. The remedy, already mentioned in [3] in the case of Gauss–Lobatto formulae, was implemented in the HOGGRL routine gradau_jacobi.m [8]. Table 2 compares values of w 0 and w 0 obtained using both gradau.m and gradau_jacobi.m with the explicit expressions (37), (36) (whose values for (α , β) = (1, 0) are also listed in Table 1).

582

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

Table 2 Comparison between the weights w 0 and w 0 obtained using gradau.m (top line), gradau_jacobi.m (middle line) and explicit expressions from (37) and (36) (bottom line) for α = 1, β = 0. N

w0

w 0

8

1.058585858585866(−1) 1.058585858585865(−1) 1.058585858585859(−1) 9.221228598390077(−3) 9.221228598390016(−3) 9.221228598391229(−3) 6.311356962063987(−4) 6.311356962064164(−4) 6.311356962055456(−4) n/a 4.037395110470794(−5) 4.037395110449350(−5)

1.616161616161623(−3) 1.616161616161622(−3) 1.616161616161616(−3) 1.198340298686251(−5) 1.198340298686246(−5) 1.198340298686319(−5) 5.602376247892514(−8) 5.602376247892348(−8) 5.602376247885541(−8) n/a 2.292287839606539(−10) 2.292287839600606(−10)

32

128

512

Fig. 1. Accuracy in evaluating (41) using various options described in the text: (1,0)

(1,1)

option (3) with P N +1 (x j ), + option (3) with P N

gradau_jacobi.m,

option (1),

option (2) with gradau.m,

option (2) with

(x j ).

For larger values of N round-off errors limit the accuracy of all stable approximations. Option (1) and option (2) using

gradau_jacobi.m deliver strikingly similar, but not identical, results. Option (3) with interior weights computed from (40) yields comparable accuracy, although somewhat worse on average (over the range 400  N  700) than options (1) or (2), as previously noted in [3] for Gauss–Lobatto rules. On the other hand, option (3) with interior weights obtained from (39) yields consistently more accurate results, with an error smaller by about two orders of magnitude. One might suspect that the simpler recurrence in the case α = β = 0 (namely an,n = 0) would account for the better round-off properties. However, a similar exercise on the integral

1

dx

(1 + x)2 + ω2

−1

= (1/ω) arctan(2/ω) (0,1)

(0,0)

instead, with α = β = 0, shows that using P N rather than P N +1 (Legendre) again yields better accuracy by two orders of magnitude. Whether this difference in accuracy is more general and why it occurs is currently under investigation. 3.1.3. Asymptotics As N → ∞ Stirling’s formula n! ∼

w 0

N →∞



C α ,β N 4 +2 β

,

w

0

N →∞



1

2π nn+ 2 e −n yields

C α ,β = 22+α +β (1 + β)!(2 + β)!.

On the other hand (37) yields

w0



(β + 2) N 2 (β + 1)(β + 3)

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

Fig. 2. Uniform convergence of

583

ϕ (x j ) = |1 − w j / w ∞ j | → 0 in [−c , c ], 0 < c < 1, for (44) (N = 2, 4, . . . , 1024). The horizontal axis represents (1 + x)/(1 − x).

so that

w0

D α ,β

N →∞



N 2 +2 β

D α ,β = 22+α +β β!(2 + β)!

,

β +2 . β +3

(42)

The distribution√of the interior nodes is known to converge, as N → ∞, to the “arcsin/Chebyshev” distribution with density function dx/(π 1 − x2 ), i.e.,

j N

N →∞

x j



dx



−1

π 1 − x2

=1−

arccos x j

π



1

=

2

+

arcsin x j



π

uniformly on any compact interval in (−1, 1) [12]. In particular

x ρ N

N →∞

∼ cos πρ

(43)

for any fixed 0 < ρ  1 (note x1 ∼ cos π / N, although 1 + x1 ∼ O (1/ N 2 )). The approximation

1

π f (x) w (x) dx

=

−1



f (cos θ) w (cos θ) 1 − cos2 θ dθ 0

N →∞



N π

N



f (cos θ j ) w (cos θ j ) 1 − cos θ 2j

j =1

with θ j = j π / N then suggests



wj

N →∞



π w (x j ) 1 − x2j N

1

=

1

π (1 − x j )α + 2 (1 + x j )β+ 2 N

≡ w∞ j ,

j = ρ N .

(44)

It is interesting to note that (44) is consistent with (42) in the sense that it yields x j ∼ O (1/ N 2+2β ) upon substituting 1 + x j ∼ O (1/ N 2 ) (occurring for small j). The convergence is illustrated in Fig. 2. 3.2. Gauss–Radau–Laguerre quadrature with double endpoint (α )

The (generalized) Laguerre polynomials P n = L n

w (x) = w

(α )

(x) = xα e −x

in the interval (a, b) = (0, ∞) and satisfy (6) with

q(x) = x,

(x) = α + 1 − x,

λn = n.

are orthogonal w.r.t. the weight function

584

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

(α )

When normalized by L n (0) = 1 these polynomials satisfy the recurrence (12) with

an,n+1 = −(n + 1 + α ),

an,n = 2n + 1 + α ,

an,n−1 = −n

for n > 0. Furthermore

γn =

∞

(α )

2

L n (x)

w (α ) (x) dx =

0

an,n−1 an−1,n

γn−1 = · · · =

(α !)2n! α !n! . γ0 = (n + α )! (n + α )!





Then (30) yields

1

( N + 2 + α )! 1 = (1 + α )!(2 + α )! N ! (1 + α )!

=

w 0

N +2+α N

(α +1)

The resulting expression of w 0 matches the Gauss–Radau weight λ0 with (5b). The ratio w 0 / w 0 from (35) becomes

w0

2(α + 2) 1 N+ . (α + 1)(α + 3) α+1

=

w 0

(45) (α +1)

The interior weights are, using (32) and noting that Q N = L N

wj =

obtained from [4, formula (5.5)], in accordance

w

0

(α + 2)x j [ L (Nα +1) (x j )]2

because of orthogonality properties and scaling,

,

or, using (33),

wj =

w 0 (α )

(α + 2)x j [ L N +1 (x j )]2

.

No reference is apparently available for Gauss–Radau–Laguerre quadratures with double endpoint, despite the relative simplicity of the formulae compared to the Jacobi case. 3.2.1. Asymptotics As N → ∞ Stirling’s formula yields

w 0



N →∞



N 2 +α

C α = (1 + α )!.

,

On the other hand (45) yields

w0

2(α + 2) N

N →∞



w

(α + 1)(α + 3)

0

so that

w0



N →∞



N 1 +α

D α = 2α !

,

α+2 . α+3

(46)

An adaptation of the results in [12] yields the estimate

j



N →∞

N

x j



0



dx





1+x

=

1 + xj − 1

π

uniformly on any compact interval in (0, ∞). In particular

x ρ √ N

N →∞

∼ πρ (2 + πρ )

(47)

for any fixed 0 < ρ and N  N 0 = ρ 2 . Numerical evidence suggests

wj

N →∞





π w (x j ) x j √

N

=

α + 12 −x j e

πxj



N

≡ w∞ j ,

j = ρ N ,

(48)

see Fig. 3. Again (48) is consistent with (44) in the sense that it yields x j ∼ O (1/ N 1+α ) upon substituting x j ∼ O (1/ N ) (which occurs for small j and can be observed in Fig. 3).

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

Fig. 3. Uniform convergence of

585

ϕ (x j ) = |1 − w j / w ∞ | → 0 in [1/c , c ], c > 1, for (48) (N = 2, 4, . . . , 1024) (note: for large N values of ϕ (x j ) associated to j

nodes x j  800 are missing because w j ≈ 0 within machine underflow).

Acknowledgements The author would like to thank the reviewers for useful suggestions and very careful check of all formulae, including the Matlab codes in Section 3.1.2. Their comments also led to the discovery of strong differences in accuracy in the evaluation of mathematically equivalent expressions for interior weights. These discrepancies are currently under investigation. Appendix A We derive expressions for the successive derivatives of P n , Q n and R n at x = a. Because q(a) = 0 the successive derivatives of P n evaluated at x = a can easily obtained by differentiating the Sturm–Liouville problem (6) multiple times and setting x = a. We obtain (with P n (a) = 1)

q P n +  P n + λn P n = 0

λn , (a)     q P n +  + q P n + λn +  P n = 0 ⇒

P n (a) = −

(A.1)

λn +  λn (λn +  ) , P n (a) =  (a) + q (a) (a)((a) + q (a))     q P n +  + 2q P n + λn + 2 + q P n = 0 ⇒

P n (a) = −



P n (a) = −

λn + 2 + q  λn (λn +  )(λn + 2 + q ) . P n (a) = −  (a) + 2q (a) (a)((a) + q (a))((a) + 2q (a))

If P n+1 − P n = αn (x − a) Q n with Q n (a) = 1 then successive differentiations yield







αn = P n +1 (a) − P n (a) =

αn Q n + (x − a) Q n = P n +1 − P n

2Q n + (x − a) Q n =



Q n (a) =

λn − λn+1 , (a)

(A.2)

P n+1 − P n

αn

P n+1 (a) − P n (a) 2αn

=−

λn + λn+1 +  , 2((a) + q (a))

(A.3)

586

B.D. Welfert / Applied Numerical Mathematics 60 (2010) 574–586

3Q n + (x − a) Q n =



Q n (a) =

P n+1 − P n

αn

P n+1 (a) − P n (a) 3αn

=

λn2 + λn λn+1 + λn2+1 + (3 + q )(λn + λn+1 ) + (2 + q ) 3((a) + q (a))((a) + 2q (a))

.

If Q n+1 − Q n = (x − a) R n then successive differentiations yield

R n + (x − a) R n = Q n +1 − Q n



λn − λn+2

R n (a) = Q n +1 (a) − Q n (a) = , 2((a) + q (a))

2R n + (x − a) R n = Q n+1 − Q n Q  (a) − Q n (a) (λn − λn+2 )(λn + λn+1 + λn+2 + 3 + q ) ⇒ R n (a) = n+1 =− . 2 6((a) + q (a))((a) + 2q (a)) Therefore

R n (a) R n (a)

=−

λn + λn+1 + λn+2 + 3 + q n((n + 1)q /2 +  ) = . 3((a) + 2q (a)) (a) + 2q (a)

(A.4)

References [1] W. Gautschi, S. Li, Gauss–Radau and Gauss–Lobatto quadratures with double end points, J. Comp. Appl. Math. 34 (1991) 343–360. [2] W. Gautschi, Algorithm 726: ORTHPOL – A package routines for generating orthogonal polynomials and Gauss-type quadrature rules, ACM Trans. Math. Soft. 20 (1994) 21–62. [3] W. Gautschi, High-order Gauss–Lobatto formulae, Num. Alg. 25 (2000) 213–225. [4] W. Gautschi, Gauss–Radau formulae for Jacobi and Laguerre weight functions, Math. Comp. Sim. 54 (2000) 403–412. [5] W. Gautschi, OPQ: Matlab suite for generating orthogonal polynomials and related quadrature rules, http://www.cs.purdue.edu/archives/2002/wxg/ codes/OPQ.html, 2002. [6] W. Gautschi, Orthogonal Polynomials, Computation and Approximation, Numerical Mathematics and Scientific Computations, Oxford University Press, Oxford, 2004. [7] W. Gautschi, Generalized Gauss–Radau and Gauss–Lobatto formulae, BIT Num. Math. 44 (2004) 711–720. [8] W. Gautschi, HOGGRL: high-order generalized Gauss–Radau and Gauss–Lobatto formulae for Jacobi and Laguerre weight functions, http://www.cs. purdue.edu/archives/2002/wxg/codes/HOGGRL.html, 2009. [9] W. Gautschi, High-order generalized Gauss–Radau and Gauss–Lobatto formulae for Jacobi and Laguerre weight functions, Num. Alg. 51 (2009) 143–149. [10] D. Gómez-Ullate, N. Kamran, R. Milson, An extended class of orthogonal polynomials defined by a Sturm–Liouville problem, J. Math. Anal. Appl. 359 (2009) 352–367. [11] M. Masjed-Jamei, A basic class of symmetric orthogonal polynomials using the extended Sturm–Liouville theorem for symmetric functions, J. Math. Anal. Appl. 325 (2007) 753–773. [12] P. Nevai, Orthogonal polynomials, Mem. Amer. Math. Soc. 213 (18) (1979), Amer. Math. Soc., Providence, RI. [13] B. Welfert, A note on classical Gauss–Radau and Gauss–Lobatto quadratures, Applied Numerical Mathematics (2010), in press, doi:10.1016/j.apnum.2010.03.006.