- Email: [email protected]
Weighted integration of periodic functions on the real line
Applied Mathematics and Computation 128 (2002) 365–378 www.elsevier.com/locate/amc
Weighted integration of periodic functions on the real line q Giuseppe Mastroianni a, Gradimir V. Milovanovic
b,*
a
b
Dipartimento di Matematica, Universit a degli Studi della Basilicata, 85100 Potenza, Italy Faculty of Electronic Engineering, Department of Mathematics, University of Nis, P.O. Box 73, 18000 Nis, Serbia, Yugoslavia
Abstract Integration of periodic functions on the real line with an even rational weight function is considered. A transformation method of such integrals to the integrals on ð1; 1Þ with respect to the Szeg} o–Bernstein weights and a construction of the corresponding Gaussian quadrature formulas are given. The recursion coefficients in the three-term recurrence relation for the corresponding orthogonal polynomials were obtained in an analytic form. Numerical examples are also included. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Gauss-type quadratures; Error term; Convergence; Orthogonal polynomials; Nonnegative measure; Weights; Chebyshev weight; Szeg} o–Bernstein weights; Nodes; Modified moments; Chebyshev polynomials
1. Introduction We consider integrals of ð2pÞ-periodic functions over the real line R,
Iðf Þ ¼
Z
f ðtÞwðtÞ dt;
ð1:1Þ
R
q Research of the authors supported by GNIM Project 2000: Teoria costruttiva delle funzioni e problemi connessi. * Corresponding author. E-mail address: [email protected] (G.V. Milovanovic).
0096-3003/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 1 ) 0 0 0 8 0 - 7
366
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
with a given even rational weight function of the form wðtÞ ¼
P ðt2 Þ ; Qðt2 Þ
ð1:2Þ
where n Y
QðtÞ ¼
t þ b2k ;
0 < b1 6 b2 6 6 bn ;
k¼1
and P ðtÞ is a polynomial of degree at most m < n, which is nonnegative on the half line ½0; þ1Þ. Problem (1.1) can be simplified by first obtaining the partial fraction decomposition of (1.2) in the form wðtÞ ¼
rj m X X j¼1
m¼1
Cjm m; ðt2 þ b2j Þ
where the sum is over all pairs of conjugate complexPpoles ibj of Qðt2 Þ, with m corresponding multiplicities rj ðj ¼ 1; . . . ; mÞ. Here, j¼1 rj ¼ n. Thus, without loss of generality we can consider only weights of the form wm ðtÞ ¼ wm ðt; bÞ ¼
1 m ðt2 þ b2 Þ
ðm P 1Þ;
ð1:3Þ
i.e., integrals Im ðf Þ ¼ Im ðf ; bÞ ¼
Z
f ðtÞ
R
dt m ðt2 þ b2 Þ
ðb > 0; m P 1Þ:
ð1:4Þ
The paper is organized as follows. In Section 2 we develop a transformation method for reducing the previous integrals Im ðf ; bÞ to the integrals on ð1; 1Þ with respect to the Szeg} o–Bernstein weights (SBWs). Section 3 is devoted to the corresponding Gaussian formulas. For appropriate values of m we obtain the explicit expressions for the recursion coefficients in the three-term-recurrence relation for the corresponding (monic) orthogonal polynomials. Finally, some numerical examples are considered in Section 4.
2. Reduction of integrals to a finite interval In this section we will show how to reduce the integral (1.4) to an integral on the finite interval. For this purpose we need the sum of the following series Wm ðsÞ ¼ Wm ðs; bÞ ¼
þ1 X k¼1
wm ð2kp þ sÞ ¼
þ1 X k¼1
h
1 2
ð2kp þ sÞ þ b2
im :
ð2:1Þ
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
367
Since (cf. [6, p. 685]) þ1 X
1 2
k¼1
ðk þ aÞ þ b
2
¼
p sinh 2pb ; b cosh 2pb cos 2pa
in the simplest, but the most important case m ¼ 1, for 2pa ¼ s and 2pb ¼ b, we obtain W1 ðsÞ ¼ W1 ðs; bÞ ¼
sinh b 1 : 2b cosh b cos s
ð2:2Þ
In a general case we can prove: Lemma 2.1. Let wm be given by (1.3), n ¼ ðs ibÞ=ð2pÞ, and f ¼ nþ . Then 12m
ð2pÞ Wm ðsÞ ¼ 2ðm 1Þ!
(
" #
) dm1 cot pz dm1 cot pz lim m1 þ lim m1 : m z!n dz ðz þ fÞ z!nþ dz ðz þ fÞm ð2:3Þ
The proof of this result can be done by an integration of the function z 7! gðzÞ ¼ p cotðpzÞwm ð2pz þ sÞ over the rectangular contour CN with vertices at the points ðN þ ð1=2ÞÞð 1 iÞ, where N 2 N is such that the poles n of the function g are inside of CN . Then, taking N ! þ1, the corresponding integral over CN tends to zero, because wm ðzÞ ¼ Oð1=z2m Þ when z ! 1. Then, by Cauchy’s residue theorem, we get Wm ðsÞ ¼
þ1 X
wm ð2kp þ sÞ ¼ Resþ gðzÞ þ Res gðzÞ ;
k¼1
z¼n
z¼n
i.e., (2.3). For m ¼ 1, (2.3) reduces to (2.2). When m ¼ 2 we have W2 ðsÞ ¼
b cosh b sinh b cos s þ a ; 2 3 4b ðcosh b cos sÞ
where a¼
sinh 2b 2b : 2b cosh b 2 sinh b
368
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
Using MA T H E M A T I C A package: In½1 :¼ g½z ;t ; b ;n :¼ Pi Cot½Pi z =ðð2Pi z þ tÞ^2 þ b^2Þ^n In½2 :¼ suma½t ;b ; n :¼ ðResidue½g½z;t;b; n ;z;ðt þ I bÞ=ð2PiÞ þ Residue½g½z;t; b;n ;z;ðt IbÞ=ð2PiÞ Þ In½3 :¼ pol½t ; b ;n :¼ ComplexExpand½ suma½t;b;n
ðCosh½b Cos½t Þ^n==Simplify we can suspect the following form of our sum Wm ðsÞ ¼
pm ðcos sÞ ; ðcosh b cos sÞm
where pm ðxÞ is an algebraic polynomial. Indeed, we can prove the following result: Theorem 2.2. Let x ¼ cos s and c ¼ cosh b. Then Wm ðsÞ ¼ Wm ðs; bÞ ¼
pm ðxÞ m ðc xÞ
ðm ¼ 1; 2; . . .Þ;
ð2:4Þ
where pm ðxÞ ¼ pm ðx; bÞ is a nonnegative polynomial on ½1; 1 of degree m 1. These polynomials satisfy the recurrence relation pffiffiffiffiffiffiffiffiffiffiffiffiffi 1 opm ðxÞ m c2 1pm ðxÞ ðc xÞ ; ð2:5Þ pmþ1 ðxÞ ¼ 2bm ob pffiffiffiffiffiffiffiffiffiffiffiffiffi where p1 ðxÞ ¼ c2 1=ð2bÞ. Proof. We start with (2.2) written in the form ðc xÞW1 ðsÞ ¼ p1 ðxÞ, where pffiffiffiffiffiffiffiffiffiffiffiffiffi sinh b c2 1 ; x ¼ cos s; c ¼ cosh b: ¼ p1 ðxÞ ¼ 2b 2b Thus, the formula (2.4) is true for m ¼ 1. Suppose that (2.4) holds for some mð P 1Þ. Then, differentiating m
ðc xÞ Wm ðsÞ ¼ pm ðxÞ with respect to b, we get mðc xÞ
m1
dc opm ðxÞ m oWm ðsÞ Wm ðsÞ þ ðc xÞ ¼ ; db ob ob
from which it follows pffiffiffiffiffiffiffiffiffiffiffiffiffi opm ðxÞ m mþ1 ; m c2 1ðc xÞ Wm ðsÞ 2bmðc xÞ Wmþ1 ðsÞ ¼ ðc xÞ ob
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
369
i.e., ðc xÞ
mþ1
Wmþ1 ðsÞ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffi 1 opm ðxÞ m c2 1pm ðxÞ ðc xÞ ¼: pmþ1 ðxÞ: 2bm ob
Thus, the result is proved.
We are ready now to give a transformation of the integral (1.3) to one on a finite interval. Putting t ¼ 2kp þ s and using the periodicity of the function f, f ðtÞ ¼ f ð2kp þ sÞ ¼ f ðsÞ; we have Im ðf Þ ¼ Im ðf ; bÞ ¼
þ1 Z X k¼1
¼ ¼
Z
f ðtÞwm ðtÞ dt
ð2k1Þp
þ1 Z X
k¼1
ð2kþ1Þp
p
p
f ðsÞ p
f ðsÞwm ð2kp þ sÞ ds
p þ1 X
! wm ð2kp þ sÞ ds;
k¼1
because of the uniform convergence of the series (2.1). Thus, Z p Im ðf Þ ¼ f ðsÞWm ðsÞ ds; p
where Wm ðsÞ is defined by (2.1) and given by (2.4). We see that Wm ðsÞ ¼ Wm ðsÞ, i.e., Wm is an even weight function. Because of the last property of the weight function, we have Z 0 Z p f ðsÞWm ðsÞ ds þ f ðsÞWm ðsÞ ds Im ðf Þ ¼ Im ðf ; bÞ ¼ p 0 Z p ¼ ð f ðsÞ þ f ð sÞÞWm ðsÞ ds: 0
Changing the variables cos s ¼ x and putting f ðsÞ þ f ðsÞ ¼ F ðcos sÞ;
ð2:6Þ
we get the following result: Theorem 2.3. The integral (1.3) can be transformed to the form Z 1 pm ðxÞ dx Im ðf Þ ¼ Im ðf ; bÞ ¼ F ðxÞ m pffiffiffiffiffiffiffiffiffiffiffiffiffi ; ðc xÞ 1 x2 1
ð2:7Þ
370
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
where c ¼ cosh b, pm ðxÞ is a polynomial determined by the recurrence relation (2.5), and F is defined by (2.6).
3. Gaussian type formulae for Szeg} o–Bernstein weights In order to evaluate the integral (2.7) it would seem more natural and simpler to apply the Gauss–Chebyshev quadrature formula, i.e., taking m x 7! UðxÞ ¼ F ðxÞpm ðxÞ=ðc xÞ ðc > 1Þ as an integrating function with respect to the Chebyshev weight (ChW) v0 ðxÞ ¼ ð1 x2 Þ1=2 . when for some r P 1 the function F satisfies the condition ffi pffiffiffiffiffiffiffiffiffiffiffiffi R 1 In ðrÞthis case, 2 Þðr1Þ dx < þ1, the error R ðUÞ F ðxÞð 1 x n v0 of the n-point Gauss– 1 Chebyshev quadrature can be estimated as follows (see [5])
Z A 1 dr F ðxÞpm ðxÞ ðr1Þ=2 ð1 x2 Þ dx; jRn ðUÞv0 j 6 r n 1 dxr ðc xÞm where A > 0 is a constant independent on U and n. Hence, when c > 1 is very close to 1, even if the integrand is bounded, it gives a very large bound. 1=2 m On the contrary, if we take vm ðxÞ ¼ ð1 x2 Þ =ðc xÞ as a weight function (Szeg} o–Bernstein weight), then the error of the corresponding Gaussian formula is bounded as follows Z B 1 dr dx ðr1Þ=2 ½ F ðxÞpm ðxÞ ð1 x2 Þ jRn ðWÞj 6 r m; n 1 dxr ðc xÞ where B > 0 is a constant independent on W ðWðxÞ ¼ F ðxÞpm ðxÞÞ and n. It is clear that the last integral is much smaller then the previous one. Also, some numerical evidences confirm this argument (see Section 4). Thus, for evaluating the integral (2.7) it is more convenient to construct the Gaussian quadratures Z 1 n X ðnÞ ðnÞ WðxÞ dkm ðxÞ ¼ Ak Wðxk Þ þ Rn ðWÞvm ; Rn ðP2n1 Þvm 0; ð3:1Þ 1
k¼1
for the measure dkm ðxÞ ¼ vm ðxÞ dx ¼
dx pffiffiffiffiffiffiffiffiffiffiffiffiffi m ðc xÞ 1 x2
ðm P 1Þ;
ð3:2Þ
where the function W includes the algebraic polynomial pm ðxÞ, i.e., WðxÞ ¼ F ðxÞpm ðxÞ. Here, P2n1 denotes the set of all polynomials of degree at most 2n 1. It is well known that the corresponding orthogonal polynomials pn;m ðxÞ for the measure (3.2) can be calculated explicitly provided m < 2n (cf. [7, p. 31]). On
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
371
the other side, there is a nonlinear algorithm to produce the recursion coefficients in the three-term recurrence relation for the monic polynomials pn;m ðxÞ, n P 0; pnþ1;m ðxÞ ¼ ðx anðmÞ Þpn;m ðxÞ bðmÞ n pn1;m ðxÞ; Z 1 ðmÞ ðmÞ b0 , m0 ¼ dkm ðxÞ p0;m ðxÞ ¼ 1; p1;m ðxÞ
ð3:3Þ
1
in terms of ones for the polynomials pn;m1 ðxÞ orthogonal with respect to the measure dkm1 ðxÞ ¼ dkm ðxÞ=ðc xÞ. However, such an algorithm is quite numerically unstable unless c is very close to the support interval of the measure (see [3, p. 102]). Two numerical algorithms for this purpose were also discussed in [1]. Our goal in this paper is to find analytic expressions for the recursion coefficients for some appropriate values of m. ðmÞ ðmÞ Knowing these coefficients, ak , bk , k P 0, one can easily obtain the n-point ðnÞ Gaussian quadrature formula (3.1) for any n. The nodes xk , indeed, are the eigenvalues of the symmetric (tridiagonal) Jacobi matrix qffiffiffiffiffiffiffi 3 2 ðmÞ b1 O am0 7 6 7 6 qffiffiffiffiffiffiffi qffiffiffiffiffiffiffi 7 6 ðmÞ ðmÞ ðmÞ 7 6 b1 a b 1 2 7 6 7 6 qffiffiffiffiffiffiffi 7 6 . ðmÞ ðmÞ 7 6 . Jn ðdkm Þ ¼ 6 b2 . a2 7 7 6 7 6 ffiffiffiffiffiffiffiffi ffi q .. .. 6 ðmÞ 7 7 6 . b . n1 7 6 5 4 qffiffiffiffiffiffiffiffiffi ðmÞ ðmÞ O bn1 an1 ðnÞ
ðnÞ
ðmÞ
while the weights (Christoffel numbers) Ak are given by Ak ¼ b0 v2k;1 in terms of the first components vk;1 of the corresponding normalized eigenvectors (cf. [2, Section 5.1; 4]). Firstly, we introduce the modified moments for dkm ðxÞ by the orthogonal polynomials pn;m1 ðxÞ, Z 1 Z 1 pn;m1 ðxÞ dx pffiffiffiffiffiffiffiffiffiffiffiffiffi ðn P 0Þ: ¼ p ðxÞ dk ðxÞ ¼ mðmÞ ð3:4Þ n;m1 m n m ðc xÞ 1 x2 1 1 Notice that pn;0 ðxÞ are the monic Chebyshev polynomials of the first kind T^n ðxÞ (T^0 ðxÞ ¼ 1, T^n ðxÞ ¼ 21n cosðn arccos xÞ, n P 1). It is easy to prove the following auxiliary result: Lemma 3.1. For the first moment we have Z 1 dx pQm1 ðcÞ ðmÞ pffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ m0 ¼ ; m 2 1x ðc2 1Þm1=2 1 ðc xÞ
372
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
where Qm ðcÞ ¼
1 ð2m 1ÞcQm1 ðcÞ ðc2 1ÞQ0m1 ðcÞ ; m
Q0 ðcÞ ¼ 1:
Thus, we find 1 3 Q2 ðcÞ ¼ c2 þ ; Q3 ðcÞ ¼ c3 þ c; 2 2 3 4 2 Q4 ðcÞ ¼ c þ 3c þ ; etc: 8 According to [6, p. 415] we have c m=2 Qm ðcÞ ¼ ðc2 1Þ Pm pffiffiffiffiffiffiffiffiffiffiffiffiffi ; c2 1 Q1 ðcÞ ¼ c;
where Pm is the Legendre polynomial of the order m. In order to get connection with the Chebyshev measure dk0 ðxÞ it is conveð0Þ nient to put Q1 ðcÞ ¼ ðc2 1Þ1=2 . Then it gives m0 ¼ p. Now, we can prove: Theorem 3.2. The polynomials pn;m ðxÞ can be expressed in terms of polynomials fpk;m1 ðxÞg in the form pn;m ðxÞ ¼ pn;m1 ðxÞ qðmÞ n pn1;m1 ðxÞ;
ð3:5Þ
ðmÞ
ðmÞ ðmÞ ðm1Þ where qðmÞ and bðm1Þ n n ¼ mn =mn1 and the moments mn are given by (3.4). If an are the recursion coefficients in (3.3) for polynomials fpn;m1 ðxÞg, and ðm1Þ
rm ¼
m0
ðmÞ m0
¼ ðc2 1Þ
Qm2 ðcÞ ; Qm1 ðcÞ
ð3:6Þ
where the polynomials Qm ðcÞ are defined in Lemma 3.1, then ðmÞ
ðm1Þ
q1 ¼ c a0
rm ;
ðmÞ
qnþ1 ¼ c aðm1Þ n
bnðm1Þ ðmÞ
qn
ðn P 1Þ:
The coefficients in (3.3) are given by ðmÞ
ðm1Þ
a0 ¼ a 0 ðmÞ
ðmÞ
þ q1 ;
ðmÞ
anðmÞ ¼ anðm1Þ þ qnþ1 qðmÞ n
ðn P 1Þ
ðmÞ
and b0 , m0 ¼ pQm1 ðcÞ=ðc2 1Þm1=2 , h i ðm1Þ ðmÞ ðm1Þ ðmÞ ðm1Þ ðmÞ bðmÞ ¼ b þ q a þ q q a n1 nþ1 n n n n n Alternatively, ðm1Þ
bðmÞ n ¼ bn1
qðmÞ n ðmÞ
qn1
ðn P 2Þ:
ðn P 1Þ:
ð3:7Þ
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
373
Proof. Putting pn;m ðxÞ ¼ pn;m1 ðxÞ
n1 X
ðmÞ
qn;k pk;m1 ðxÞ
k¼0
and using the inner product with respect to the measure dkm1 , Z 1 ðf ; gÞm1 ¼ f ðxÞgðxÞ dkm1 ; 1
because of orthogonality, we obtain that for each 0 6 i 6 n 2, ðmÞ
ðpn;m ; pi;m1 Þm1 ¼ qn;i ðpi;m1 ; pi;m1 Þm1 and ðpn;m ; pi;m1 Þm1 ¼
Z
1
ðc xÞpn;m ðxÞpi;m1 ðxÞ dkm ðxÞ Z 1 pn;m ðxÞð xpi;m1 ðxÞÞ dkm ðxÞ ¼ 0: ¼ 1
1 ðmÞ qn;i
Thus, we conclude that ¼ 0 for such values of i and the formula (3.5) is ðmÞ true, where we put qn;n1 qðmÞ n . From (3.5), because of orthogonality 0 ¼ ðpn;m ; 1Þm ¼ ðpn;m1 ; 1Þm qðmÞ n ðpn1;m1 ; 1Þm ; ðmÞ
ðmÞ we get qðmÞ n ¼ mn =mn1 , where the modified moments are defined by (3.4). Using the recurrence relation for polynomials fpn;m1 ðxÞg we find that Z 1 ðmÞ ðm1Þ ðmÞ ðmÞ mnþ1 ¼ ðc aðm1Þ Þm b m pn;m1 ðxÞ dkm1 ; n1 n n n 1
which gives ðmÞ
ðm1Þ
m1 ¼ ðc a0
ðmÞ
ðm1Þ
Þm0 m0
and ðmÞ
ðmÞ
ðm1Þ ÞmðmÞ mn1 mnþ1 ¼ ðc aðm1Þ n n bn
ðn P 1Þ:
These equalities give (3.7). Finally, changing pk;m ðxÞ ðk ¼ n 1; n; n þ 1Þ in the recurrence relation (3.3) by (3.5) and using the corresponding relation for polynomials fpk;m1 ðxÞg we get for n P 2 ! " ! ðmÞ ðmÞ ðmÞ ðmÞ pnþ1;m1 ðxÞ ¼ x aðmÞ pn;m1 ðxÞ bðmÞ n þ qnþ1 qn n qn an " ! " ðm1Þ ðm1Þ ðmÞ ðmÞ ðmÞ þ qðmÞ a ðxÞ þ b q q p b pn2;m1 ðxÞ: n1;m1 n1 n1 n1 n n n
374
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
Comparing with the recurrence relation for fpn;m1 ðxÞg we obtain formulas for the recursion coefficients. The case n ¼ 1 should be considered separately. Notice that in Chebyshev case ðm ¼ 0Þ we have anð0Þ ¼ 0
ð0Þ
ðn P 0Þ;
b0 , p;
1 ð0Þ b1 ¼ ; 2
bnð0Þ ¼
1 4
ðn P 2Þ:
Also, from (3.6) it follows r1 ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffi c 2 1;
r4 ¼
ðc2 þ ð1=2ÞÞðc2 1Þ ; etc: c3 þ ð3=2Þc
r2 ¼
c2 1 ; c
r3 ¼
cðc2 1Þ ; c2 þ ð1=2Þ
Now, using the previous theorem we give explicit expressions for recursion coefficients for some important special cases, where c ¼ cosh b. ð1Þ
Case m ¼ 1. Here we have q1 ¼ eb ; qnð1Þ ¼ ð1=2Þeb ðn P 2Þ, and the recursion coefficients ð1Þ
1 ð1Þ a1 ¼ eb ; anð1Þ ¼ 0 ðn P 2Þ; 2 p 1 1 ð1Þ ; b1 ¼ 1 e2b ; bnð1Þ ¼ ðn P 2Þ: , sinh b 2 4
a0 ¼ eb ; ð1Þ
b0
ð2Þ
b Case m ¼ 2. Here, q1 ¼ eb tanh b, qð2Þ ðn P 2Þ, and n ¼ ð1=2Þe ð2Þ
a0 ¼ ð2Þ
b0 ,
1 ; cosh b p cosh b ; sinh3 b
ð2Þ
a1 ¼ eb tanh b; ð2Þ
b1 ¼
anð2Þ ¼ 0
1 1 e2b tanh2 b; 2
ðn P 2Þ; ð2Þ
b2 ¼
1 1 þ e2b ; 4
and bnð2Þ ¼ ð1=4Þ ðn P 3Þ. ð3Þ
ð3Þ
Case m ¼ 3. Here, q1 ¼ sinh b tanh b=ð2 þ coshð2bÞÞ, q2 ¼ e2b cosh b, ¼ 12eb ðn P 3Þ, and 3 cosh b sinh b ð3Þ ð3Þ ; a1 ¼ e2b cosh b þ eb þ tanh b; a0 ¼ 2 þ coshð2bÞ 2 þ coshð2bÞ
qnð3Þ
1 ð3Þ a2 ¼ e3b ; anð3Þ ¼ 0 ðn P 3Þ; 2 4 p cosh2 b þ ð1=2Þ ð1 e2b Þ ð3Þ ð3Þ ; b b0 , ¼ ; 1 sinh5 b 2ð1 4e2b þ e4b Þ2
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378 ð3Þ
b2 ¼
1 1 þ 3e2b 3e4b e6b ; 4
bð3Þ n ¼
1 4
375
ðn P 3Þ:
As we can see the recurrence coefficients for polynomials pn;m ðxÞ reduce to the corresponding coefficients for Chebyshev polynomials for n P n0 ðn 2 NÞ. Precisely, calculations show that
mþ1 ð0Þ aðmÞ ¼ a ¼ 0; n P þ1 n n 2 and 1 ð0Þ bðmÞ n ¼ bn ¼ ; 4
nP
hmi 2
þ 1:
4. Numerical examples In order to illustrate the presented transformation method, we consider in this section a few numerical examples. All computations were done in Darithmetic on the WO R K S T A T I O N DI G I T A L UL T I M A T E AL P H A 533au2 (with machine precision 2:22 1016 ). Example 4.1. Consider integrals of the form Z þ1 2 sin 2t 1 e cos 2t 2 Im ðf ; bÞ ¼ dt 2 m 1 3 þ 2 cos 3t ðt þ b Þ
ðm P 1Þ:
The function f ðtÞ ¼
2 sin 2t 1 cos 2t e 3 þ 2 cos 3t
is ð2pÞ-periodic and its graph on the interval ½p; p is displayed in Fig. 1.
Fig. 1. The periodic function f ðtÞ (a). The function F ðxÞ obtained by transformation (2.6) (b).
376
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
Since f ðsÞ þ f ðsÞ ¼
2e cos 2s ; 3 þ 2 cos 3s
putting x ¼ cos s and using (2.6) we find 2
F ðxÞ ¼ F ðcos sÞ ¼ f ðsÞ þ f ðsÞ ¼
2e12x 3 6x þ 8x3
and according to (2.7), Z 1 pm ðxÞ dx Im ðf ; bÞ ¼ F ðxÞ m pffiffiffiffiffiffiffiffiffiffiffiffiffi ; ðc xÞ 1 x2 1 where c ¼ cosh b. Let m ¼ 1. Applying Gaussian quadratures with the ChW for n ¼ 5ð5Þ50 and taking b ¼ 10m ðm ¼ 2; 1; 0Þ we get approximations of I1 ðf ; bÞ with relative errors given in Table 1. Numbers in parentheses indicate decimal exponents. Taking Gaussian quadratures for n ¼ 5ð5Þ50, with SBW, v1 ðxÞ ¼ 1=2 ð1 x2 Þ =ðc xÞ, the corresponding errors are also presented in the same table. The corresponding exact values of I1 ðf ; bÞ are obtained using Gaussian quadratures with SBW in Q-arithmetic (machine precision 1:93 1034 ): I1 ðf ; 0:01Þ ¼ 0:2586588216241823127882 . . . 102 ðc ¼ 1:0000500 . . .Þ; I1 ðf ; 0:10Þ ¼ 0:4968012877996286228355 . . . 101 ðc ¼ 1:0050041 . . .Þ; I1 ðf ; 1:00Þ ¼ 0:1673215409745331112726 . . . 101 ðc ¼ 1:5430806 . . .Þ: Table 1 Relative errors in Gaussian approximations of the integral I1 ðf ; bÞ with respect to the Chebyshev weight (ChW) and the Szeg} o–Bernstein weight (SBW) b
b ¼ 0:01
n
ChW
SBW
ChW
SBW
ChW
SBW
5 10 15 20 25 30 35 40 45 50
8:4 ð1Þ 8:0 ð1Þ 7:6 ð1Þ 7:2 ð1Þ 6:7 ð1Þ 6:3 ð1Þ 5:9 ð1Þ 5:5 ð1Þ 5:2 ð1Þ 4:8 ð1Þ
1:3 ð2Þ 2:4 ð4Þ 1:1 ð5Þ 9:0 ð7Þ 1:6 ð8Þ 7:4 ð10Þ 5:9 ð11Þ 1:0 ð12Þ 4:8 ð14Þ 4:7 ð15Þ
2:2 ð1Þ 1:1 ð1Þ 4:3 ð2Þ 1:6 ð2Þ 6:0 ð3Þ 2:2 ð3Þ 8:2 ð4Þ 3:0 ð4Þ 1:1 ð4Þ 4:1 ð5Þ
6:3 ð2Þ 1:5 ð3Þ 4:2 ð5Þ 4:4 ð6Þ 9:7 ð8Þ 2:8 ð9Þ 2:9 ð10Þ 6:4 ð12Þ 1:8 ð13Þ 1:9 ð14Þ
1:2 ð2Þ 2:7 ð3Þ 1:5 ð4Þ 1:0 ð6Þ 1:8 ð7Þ 9:8 ð9Þ 6:7 ð11Þ 1:2 ð11Þ 6:5 ð13Þ 5:2 ð15Þ
5:5 ð2Þ 3:5 ð3Þ 7:0 ð5Þ 3:5 ð6Þ 2:3 ð7Þ 4:6 ð9Þ 2:3 ð10Þ 1:5 ð11Þ 3:1 ð13Þ 1:6 ð14Þ
b ¼ 0:1
b ¼ 1:0
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
377
As can be seen, for smaller values of b (c is close to 1) the Gauss–Chebyshev quadratures ChW cannot be used directly. When b increases the both quadratures become comparable. However, by writing I1 ðf ; bÞ in the form Z p sinh b 1 F ðcÞ F ðxÞ dx pffiffiffiffiffiffiffiffiffiffiffiffiffi ; I1 ðf ; bÞ ¼ F ðcÞ 2b 2b c x 1 x2 1 the Gauss–Chebyshev quadratures can be applied directly. Consider now the case m ¼ 2, with the functions /k and the corresponding weights vk ðk ¼ 0; 1; 2Þ, where /k ðxÞ ¼
F ðxÞp2 ðxÞ ðc xÞ
2k
;
vk ðxÞ ¼
1 pffiffiffiffiffiffiffiffiffiffiffiffiffi : ðc xÞ 1 x2 k
Applying the Gaussian quadratures with the Chebyshev weight ChW (k ¼ 0) and the Szeg} o–Bernstein weights SBW1 (k ¼ 1) and SBW2 (k ¼ 2) we get approximations of the integral I2 ðf ; bÞ. The exact values of this integral for some selected b are: I2 ðf ; 0:01Þ ¼ 0:1156183821140487028202 . . . 106 ; I2 ðf ; 0:10Þ ¼ 0:1214706913588412300593 . . . 103 : The relative errors in Gaussian approximations for n ¼ 5ð5Þ50 are presented in Table 2. The advantage of quadrature formulas for k ¼ 2 (in this case m ¼ 2) is evident. When b increases all quadratures give similar results. Example 4.2. Consider now the integral (1.4), with a nonanalytic function 2a f ðtÞ ¼ j cosðt=2Þj (a > 0). After the transformation we obtain the integral Table 2 Relative errors in Gaussian approximations of the integral I2 ðf ; bÞ with respect to the Chebyshev weight (ChW) and to the Szeg} o–Bernstein weights (SBW1 and SBW2 ) b
b ¼ 0:01
n
ChW
SBW1
SBW2
ChW
SBW1
SBW2
5 10 15 20 25 30 35 40 45 50
1:0 ð0Þ 1:0 ð0Þ 1:0 ð0Þ 9:9 ð1Þ 9:9 ð1Þ 9:8 ð1Þ 9:7 ð1Þ 9:6 ð1Þ 9:5 ð1Þ 9:3 ð1Þ
9:1 ð1Þ 8:3 ð1Þ 7:5 ð1Þ 6:8 ð1Þ 6:1 ð1Þ 5:5 ð1Þ 5:0 ð1Þ 4:5 ð1Þ 4:1 ð1Þ 3:7 ð1Þ
5:5 ð7Þ 1:0 ð7Þ 4:7 ð9Þ 1:9 ð11Þ 4:6 ð12Þ 2:6 ð12Þ 2:3 ð12Þ 2:3 ð12Þ 2:3 ð12Þ 2:3 ð12Þ
8:9 ð1Þ 6:2 ð1Þ 3:4 ð1Þ 1:6 ð1Þ 7:4 ð2Þ 3:2 ð2Þ 1:3 ð2Þ 5:6 ð3Þ 2:3 ð3Þ 9:2 ð4Þ
3:7 ð1Þ 1:4 ð1Þ 5:0 ð2Þ 1:9 ð2Þ 6:8 ð3Þ 2:5 ð3Þ 9:2 ð4Þ 3:4 ð4Þ 1:2 ð4Þ 4:6 ð5Þ
1:1 ð3Þ 6:7 ð5Þ 4:3 ð6Þ 6:4 ð8Þ 4:4 ð9Þ 2:8 ð10Þ 4:3 ð12Þ 3:1 ð13Þ 1:4 ð15Þ 2:0 ð14Þ
b ¼ 0:1
378
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
Fig. 2. Relative errors in Gaussian approximations with n ¼ 5 (upper curve) and n ¼ 20 nodes (lower curve) for 0 6 a < 4:5.
J ðaÞ ¼ A
Z
1
1
1þx 2
a
dx pffiffiffiffiffiffiffiffiffiffiffiffiffi ; ðc xÞ 1 x2
where A ¼ 2p1 ðxÞ ¼ sinh b=b. In order to evaluate this integral, we apply Gaussian rule in n points with SBW ðm ¼ 1Þ. A typical behavior of the relative error rn of Gaussian approximations with respect to the parameter a ð0 6 a < 4:5Þ is displayed in Fig. 2 in the log-scale. Two cases for n ¼ 5 and n ¼ 20 are given, whereas b ¼ 0:01. It is clear that the rapidly increasing of accuracy achieves when the parameter a tends to an integer (i.e., when f becomes an analytic function).
References [1] B. Fischer, G. Golub, How to generate unknown orthogonal polynomials out of known orthogonal polynomials, J. Comput. Appl. Math. 43 (1992) 99–115. [2] W. Gautschi, A survey of Gauss–Christoffel quadrature formulae, in: P.L. Butzer, F. Feher (Eds.), E.B. Christoffel, Birkh€auser, Basel, 1981, pp. 72–147. [3] W. Gautschi, Orthogonal polynomials: applications and computation, Acta Numerica (1996) 45–119. [4] G. Golub, J.H. Welsch, Calculation of Gauss quadrature rules, Math. Comp. 44 (1969) 221– 230. [5] G. Mastroianni, Generalized Christoffel functions and error of positive quadrature, Numer. Algorithms 10 (1995) 113–126. [6] A.P. Prudnikov, Yu.A. Brychkov, O.I. Marichev, Integrals and Series. Elementary Functions, Nauka, Moscow, 1981 (in Russian). [7] G. Szeg} o, Orthogonal Polynomials, Amer. Math. Soc. Colloq. Publ., vol. 23, fourth ed., Amer. Math. Soc., Providence, RI, 1975.
Weighted integration of periodic functions on the real line q Giuseppe Mastroianni a, Gradimir V. Milovanovic
b,*
a
b
Dipartimento di Matematica, Universit a degli Studi della Basilicata, 85100 Potenza, Italy Faculty of Electronic Engineering, Department of Mathematics, University of Nis, P.O. Box 73, 18000 Nis, Serbia, Yugoslavia
Abstract Integration of periodic functions on the real line with an even rational weight function is considered. A transformation method of such integrals to the integrals on ð1; 1Þ with respect to the Szeg} o–Bernstein weights and a construction of the corresponding Gaussian quadrature formulas are given. The recursion coefficients in the three-term recurrence relation for the corresponding orthogonal polynomials were obtained in an analytic form. Numerical examples are also included. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Gauss-type quadratures; Error term; Convergence; Orthogonal polynomials; Nonnegative measure; Weights; Chebyshev weight; Szeg} o–Bernstein weights; Nodes; Modified moments; Chebyshev polynomials
1. Introduction We consider integrals of ð2pÞ-periodic functions over the real line R,
Iðf Þ ¼
Z
f ðtÞwðtÞ dt;
ð1:1Þ
R
q Research of the authors supported by GNIM Project 2000: Teoria costruttiva delle funzioni e problemi connessi. * Corresponding author. E-mail address: [email protected] (G.V. Milovanovic).
0096-3003/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 1 ) 0 0 0 8 0 - 7
366
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
with a given even rational weight function of the form wðtÞ ¼
P ðt2 Þ ; Qðt2 Þ
ð1:2Þ
where n Y
QðtÞ ¼
t þ b2k ;
0 < b1 6 b2 6 6 bn ;
k¼1
and P ðtÞ is a polynomial of degree at most m < n, which is nonnegative on the half line ½0; þ1Þ. Problem (1.1) can be simplified by first obtaining the partial fraction decomposition of (1.2) in the form wðtÞ ¼
rj m X X j¼1
m¼1
Cjm m; ðt2 þ b2j Þ
where the sum is over all pairs of conjugate complexPpoles ibj of Qðt2 Þ, with m corresponding multiplicities rj ðj ¼ 1; . . . ; mÞ. Here, j¼1 rj ¼ n. Thus, without loss of generality we can consider only weights of the form wm ðtÞ ¼ wm ðt; bÞ ¼
1 m ðt2 þ b2 Þ
ðm P 1Þ;
ð1:3Þ
i.e., integrals Im ðf Þ ¼ Im ðf ; bÞ ¼
Z
f ðtÞ
R
dt m ðt2 þ b2 Þ
ðb > 0; m P 1Þ:
ð1:4Þ
The paper is organized as follows. In Section 2 we develop a transformation method for reducing the previous integrals Im ðf ; bÞ to the integrals on ð1; 1Þ with respect to the Szeg} o–Bernstein weights (SBWs). Section 3 is devoted to the corresponding Gaussian formulas. For appropriate values of m we obtain the explicit expressions for the recursion coefficients in the three-term-recurrence relation for the corresponding (monic) orthogonal polynomials. Finally, some numerical examples are considered in Section 4.
2. Reduction of integrals to a finite interval In this section we will show how to reduce the integral (1.4) to an integral on the finite interval. For this purpose we need the sum of the following series Wm ðsÞ ¼ Wm ðs; bÞ ¼
þ1 X k¼1
wm ð2kp þ sÞ ¼
þ1 X k¼1
h
1 2
ð2kp þ sÞ þ b2
im :
ð2:1Þ
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
367
Since (cf. [6, p. 685]) þ1 X
1 2
k¼1
ðk þ aÞ þ b
2
¼
p sinh 2pb ; b cosh 2pb cos 2pa
in the simplest, but the most important case m ¼ 1, for 2pa ¼ s and 2pb ¼ b, we obtain W1 ðsÞ ¼ W1 ðs; bÞ ¼
sinh b 1 : 2b cosh b cos s
ð2:2Þ
In a general case we can prove: Lemma 2.1. Let wm be given by (1.3), n ¼ ðs ibÞ=ð2pÞ, and f ¼ nþ . Then 12m
ð2pÞ Wm ðsÞ ¼ 2ðm 1Þ!
(
" #
) dm1 cot pz dm1 cot pz lim m1 þ lim m1 : m z!n dz ðz þ fÞ z!nþ dz ðz þ fÞm ð2:3Þ
The proof of this result can be done by an integration of the function z 7! gðzÞ ¼ p cotðpzÞwm ð2pz þ sÞ over the rectangular contour CN with vertices at the points ðN þ ð1=2ÞÞð 1 iÞ, where N 2 N is such that the poles n of the function g are inside of CN . Then, taking N ! þ1, the corresponding integral over CN tends to zero, because wm ðzÞ ¼ Oð1=z2m Þ when z ! 1. Then, by Cauchy’s residue theorem, we get Wm ðsÞ ¼
þ1 X
wm ð2kp þ sÞ ¼ Resþ gðzÞ þ Res gðzÞ ;
k¼1
z¼n
z¼n
i.e., (2.3). For m ¼ 1, (2.3) reduces to (2.2). When m ¼ 2 we have W2 ðsÞ ¼
b cosh b sinh b cos s þ a ; 2 3 4b ðcosh b cos sÞ
where a¼
sinh 2b 2b : 2b cosh b 2 sinh b
368
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
Using MA T H E M A T I C A package: In½1 :¼ g½z ;t ; b ;n :¼ Pi Cot½Pi z =ðð2Pi z þ tÞ^2 þ b^2Þ^n In½2 :¼ suma½t ;b ; n :¼ ðResidue½g½z;t;b; n ;z;ðt þ I bÞ=ð2PiÞ þ Residue½g½z;t; b;n ;z;ðt IbÞ=ð2PiÞ Þ In½3 :¼ pol½t ; b ;n :¼ ComplexExpand½ suma½t;b;n
ðCosh½b Cos½t Þ^n==Simplify we can suspect the following form of our sum Wm ðsÞ ¼
pm ðcos sÞ ; ðcosh b cos sÞm
where pm ðxÞ is an algebraic polynomial. Indeed, we can prove the following result: Theorem 2.2. Let x ¼ cos s and c ¼ cosh b. Then Wm ðsÞ ¼ Wm ðs; bÞ ¼
pm ðxÞ m ðc xÞ
ðm ¼ 1; 2; . . .Þ;
ð2:4Þ
where pm ðxÞ ¼ pm ðx; bÞ is a nonnegative polynomial on ½1; 1 of degree m 1. These polynomials satisfy the recurrence relation pffiffiffiffiffiffiffiffiffiffiffiffiffi 1 opm ðxÞ m c2 1pm ðxÞ ðc xÞ ; ð2:5Þ pmþ1 ðxÞ ¼ 2bm ob pffiffiffiffiffiffiffiffiffiffiffiffiffi where p1 ðxÞ ¼ c2 1=ð2bÞ. Proof. We start with (2.2) written in the form ðc xÞW1 ðsÞ ¼ p1 ðxÞ, where pffiffiffiffiffiffiffiffiffiffiffiffiffi sinh b c2 1 ; x ¼ cos s; c ¼ cosh b: ¼ p1 ðxÞ ¼ 2b 2b Thus, the formula (2.4) is true for m ¼ 1. Suppose that (2.4) holds for some mð P 1Þ. Then, differentiating m
ðc xÞ Wm ðsÞ ¼ pm ðxÞ with respect to b, we get mðc xÞ
m1
dc opm ðxÞ m oWm ðsÞ Wm ðsÞ þ ðc xÞ ¼ ; db ob ob
from which it follows pffiffiffiffiffiffiffiffiffiffiffiffiffi opm ðxÞ m mþ1 ; m c2 1ðc xÞ Wm ðsÞ 2bmðc xÞ Wmþ1 ðsÞ ¼ ðc xÞ ob
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
369
i.e., ðc xÞ
mþ1
Wmþ1 ðsÞ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffi 1 opm ðxÞ m c2 1pm ðxÞ ðc xÞ ¼: pmþ1 ðxÞ: 2bm ob
Thus, the result is proved.
We are ready now to give a transformation of the integral (1.3) to one on a finite interval. Putting t ¼ 2kp þ s and using the periodicity of the function f, f ðtÞ ¼ f ð2kp þ sÞ ¼ f ðsÞ; we have Im ðf Þ ¼ Im ðf ; bÞ ¼
þ1 Z X k¼1
¼ ¼
Z
f ðtÞwm ðtÞ dt
ð2k1Þp
þ1 Z X
k¼1
ð2kþ1Þp
p
p
f ðsÞ p
f ðsÞwm ð2kp þ sÞ ds
p þ1 X
! wm ð2kp þ sÞ ds;
k¼1
because of the uniform convergence of the series (2.1). Thus, Z p Im ðf Þ ¼ f ðsÞWm ðsÞ ds; p
where Wm ðsÞ is defined by (2.1) and given by (2.4). We see that Wm ðsÞ ¼ Wm ðsÞ, i.e., Wm is an even weight function. Because of the last property of the weight function, we have Z 0 Z p f ðsÞWm ðsÞ ds þ f ðsÞWm ðsÞ ds Im ðf Þ ¼ Im ðf ; bÞ ¼ p 0 Z p ¼ ð f ðsÞ þ f ð sÞÞWm ðsÞ ds: 0
Changing the variables cos s ¼ x and putting f ðsÞ þ f ðsÞ ¼ F ðcos sÞ;
ð2:6Þ
we get the following result: Theorem 2.3. The integral (1.3) can be transformed to the form Z 1 pm ðxÞ dx Im ðf Þ ¼ Im ðf ; bÞ ¼ F ðxÞ m pffiffiffiffiffiffiffiffiffiffiffiffiffi ; ðc xÞ 1 x2 1
ð2:7Þ
370
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
where c ¼ cosh b, pm ðxÞ is a polynomial determined by the recurrence relation (2.5), and F is defined by (2.6).
3. Gaussian type formulae for Szeg} o–Bernstein weights In order to evaluate the integral (2.7) it would seem more natural and simpler to apply the Gauss–Chebyshev quadrature formula, i.e., taking m x 7! UðxÞ ¼ F ðxÞpm ðxÞ=ðc xÞ ðc > 1Þ as an integrating function with respect to the Chebyshev weight (ChW) v0 ðxÞ ¼ ð1 x2 Þ1=2 . when for some r P 1 the function F satisfies the condition ffi pffiffiffiffiffiffiffiffiffiffiffiffi R 1 In ðrÞthis case, 2 Þðr1Þ dx < þ1, the error R ðUÞ F ðxÞð 1 x n v0 of the n-point Gauss– 1 Chebyshev quadrature can be estimated as follows (see [5])
Z A 1 dr F ðxÞpm ðxÞ ðr1Þ=2 ð1 x2 Þ dx; jRn ðUÞv0 j 6 r n 1 dxr ðc xÞm where A > 0 is a constant independent on U and n. Hence, when c > 1 is very close to 1, even if the integrand is bounded, it gives a very large bound. 1=2 m On the contrary, if we take vm ðxÞ ¼ ð1 x2 Þ =ðc xÞ as a weight function (Szeg} o–Bernstein weight), then the error of the corresponding Gaussian formula is bounded as follows Z B 1 dr dx ðr1Þ=2 ½ F ðxÞpm ðxÞ ð1 x2 Þ jRn ðWÞj 6 r m; n 1 dxr ðc xÞ where B > 0 is a constant independent on W ðWðxÞ ¼ F ðxÞpm ðxÞÞ and n. It is clear that the last integral is much smaller then the previous one. Also, some numerical evidences confirm this argument (see Section 4). Thus, for evaluating the integral (2.7) it is more convenient to construct the Gaussian quadratures Z 1 n X ðnÞ ðnÞ WðxÞ dkm ðxÞ ¼ Ak Wðxk Þ þ Rn ðWÞvm ; Rn ðP2n1 Þvm 0; ð3:1Þ 1
k¼1
for the measure dkm ðxÞ ¼ vm ðxÞ dx ¼
dx pffiffiffiffiffiffiffiffiffiffiffiffiffi m ðc xÞ 1 x2
ðm P 1Þ;
ð3:2Þ
where the function W includes the algebraic polynomial pm ðxÞ, i.e., WðxÞ ¼ F ðxÞpm ðxÞ. Here, P2n1 denotes the set of all polynomials of degree at most 2n 1. It is well known that the corresponding orthogonal polynomials pn;m ðxÞ for the measure (3.2) can be calculated explicitly provided m < 2n (cf. [7, p. 31]). On
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
371
the other side, there is a nonlinear algorithm to produce the recursion coefficients in the three-term recurrence relation for the monic polynomials pn;m ðxÞ, n P 0; pnþ1;m ðxÞ ¼ ðx anðmÞ Þpn;m ðxÞ bðmÞ n pn1;m ðxÞ; Z 1 ðmÞ ðmÞ b0 , m0 ¼ dkm ðxÞ p0;m ðxÞ ¼ 1; p1;m ðxÞ
ð3:3Þ
1
in terms of ones for the polynomials pn;m1 ðxÞ orthogonal with respect to the measure dkm1 ðxÞ ¼ dkm ðxÞ=ðc xÞ. However, such an algorithm is quite numerically unstable unless c is very close to the support interval of the measure (see [3, p. 102]). Two numerical algorithms for this purpose were also discussed in [1]. Our goal in this paper is to find analytic expressions for the recursion coefficients for some appropriate values of m. ðmÞ ðmÞ Knowing these coefficients, ak , bk , k P 0, one can easily obtain the n-point ðnÞ Gaussian quadrature formula (3.1) for any n. The nodes xk , indeed, are the eigenvalues of the symmetric (tridiagonal) Jacobi matrix qffiffiffiffiffiffiffi 3 2 ðmÞ b1 O am0 7 6 7 6 qffiffiffiffiffiffiffi qffiffiffiffiffiffiffi 7 6 ðmÞ ðmÞ ðmÞ 7 6 b1 a b 1 2 7 6 7 6 qffiffiffiffiffiffiffi 7 6 . ðmÞ ðmÞ 7 6 . Jn ðdkm Þ ¼ 6 b2 . a2 7 7 6 7 6 ffiffiffiffiffiffiffiffi ffi q .. .. 6 ðmÞ 7 7 6 . b . n1 7 6 5 4 qffiffiffiffiffiffiffiffiffi ðmÞ ðmÞ O bn1 an1 ðnÞ
ðnÞ
ðmÞ
while the weights (Christoffel numbers) Ak are given by Ak ¼ b0 v2k;1 in terms of the first components vk;1 of the corresponding normalized eigenvectors (cf. [2, Section 5.1; 4]). Firstly, we introduce the modified moments for dkm ðxÞ by the orthogonal polynomials pn;m1 ðxÞ, Z 1 Z 1 pn;m1 ðxÞ dx pffiffiffiffiffiffiffiffiffiffiffiffiffi ðn P 0Þ: ¼ p ðxÞ dk ðxÞ ¼ mðmÞ ð3:4Þ n;m1 m n m ðc xÞ 1 x2 1 1 Notice that pn;0 ðxÞ are the monic Chebyshev polynomials of the first kind T^n ðxÞ (T^0 ðxÞ ¼ 1, T^n ðxÞ ¼ 21n cosðn arccos xÞ, n P 1). It is easy to prove the following auxiliary result: Lemma 3.1. For the first moment we have Z 1 dx pQm1 ðcÞ ðmÞ pffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ m0 ¼ ; m 2 1x ðc2 1Þm1=2 1 ðc xÞ
372
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
where Qm ðcÞ ¼
1 ð2m 1ÞcQm1 ðcÞ ðc2 1ÞQ0m1 ðcÞ ; m
Q0 ðcÞ ¼ 1:
Thus, we find 1 3 Q2 ðcÞ ¼ c2 þ ; Q3 ðcÞ ¼ c3 þ c; 2 2 3 4 2 Q4 ðcÞ ¼ c þ 3c þ ; etc: 8 According to [6, p. 415] we have c m=2 Qm ðcÞ ¼ ðc2 1Þ Pm pffiffiffiffiffiffiffiffiffiffiffiffiffi ; c2 1 Q1 ðcÞ ¼ c;
where Pm is the Legendre polynomial of the order m. In order to get connection with the Chebyshev measure dk0 ðxÞ it is conveð0Þ nient to put Q1 ðcÞ ¼ ðc2 1Þ1=2 . Then it gives m0 ¼ p. Now, we can prove: Theorem 3.2. The polynomials pn;m ðxÞ can be expressed in terms of polynomials fpk;m1 ðxÞg in the form pn;m ðxÞ ¼ pn;m1 ðxÞ qðmÞ n pn1;m1 ðxÞ;
ð3:5Þ
ðmÞ
ðmÞ ðmÞ ðm1Þ where qðmÞ and bðm1Þ n n ¼ mn =mn1 and the moments mn are given by (3.4). If an are the recursion coefficients in (3.3) for polynomials fpn;m1 ðxÞg, and ðm1Þ
rm ¼
m0
ðmÞ m0
¼ ðc2 1Þ
Qm2 ðcÞ ; Qm1 ðcÞ
ð3:6Þ
where the polynomials Qm ðcÞ are defined in Lemma 3.1, then ðmÞ
ðm1Þ
q1 ¼ c a0
rm ;
ðmÞ
qnþ1 ¼ c aðm1Þ n
bnðm1Þ ðmÞ
qn
ðn P 1Þ:
The coefficients in (3.3) are given by ðmÞ
ðm1Þ
a0 ¼ a 0 ðmÞ
ðmÞ
þ q1 ;
ðmÞ
anðmÞ ¼ anðm1Þ þ qnþ1 qðmÞ n
ðn P 1Þ
ðmÞ
and b0 , m0 ¼ pQm1 ðcÞ=ðc2 1Þm1=2 , h i ðm1Þ ðmÞ ðm1Þ ðmÞ ðm1Þ ðmÞ bðmÞ ¼ b þ q a þ q q a n1 nþ1 n n n n n Alternatively, ðm1Þ
bðmÞ n ¼ bn1
qðmÞ n ðmÞ
qn1
ðn P 2Þ:
ðn P 1Þ:
ð3:7Þ
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
373
Proof. Putting pn;m ðxÞ ¼ pn;m1 ðxÞ
n1 X
ðmÞ
qn;k pk;m1 ðxÞ
k¼0
and using the inner product with respect to the measure dkm1 , Z 1 ðf ; gÞm1 ¼ f ðxÞgðxÞ dkm1 ; 1
because of orthogonality, we obtain that for each 0 6 i 6 n 2, ðmÞ
ðpn;m ; pi;m1 Þm1 ¼ qn;i ðpi;m1 ; pi;m1 Þm1 and ðpn;m ; pi;m1 Þm1 ¼
Z
1
ðc xÞpn;m ðxÞpi;m1 ðxÞ dkm ðxÞ Z 1 pn;m ðxÞð xpi;m1 ðxÞÞ dkm ðxÞ ¼ 0: ¼ 1
1 ðmÞ qn;i
Thus, we conclude that ¼ 0 for such values of i and the formula (3.5) is ðmÞ true, where we put qn;n1 qðmÞ n . From (3.5), because of orthogonality 0 ¼ ðpn;m ; 1Þm ¼ ðpn;m1 ; 1Þm qðmÞ n ðpn1;m1 ; 1Þm ; ðmÞ
ðmÞ we get qðmÞ n ¼ mn =mn1 , where the modified moments are defined by (3.4). Using the recurrence relation for polynomials fpn;m1 ðxÞg we find that Z 1 ðmÞ ðm1Þ ðmÞ ðmÞ mnþ1 ¼ ðc aðm1Þ Þm b m pn;m1 ðxÞ dkm1 ; n1 n n n 1
which gives ðmÞ
ðm1Þ
m1 ¼ ðc a0
ðmÞ
ðm1Þ
Þm0 m0
and ðmÞ
ðmÞ
ðm1Þ ÞmðmÞ mn1 mnþ1 ¼ ðc aðm1Þ n n bn
ðn P 1Þ:
These equalities give (3.7). Finally, changing pk;m ðxÞ ðk ¼ n 1; n; n þ 1Þ in the recurrence relation (3.3) by (3.5) and using the corresponding relation for polynomials fpk;m1 ðxÞg we get for n P 2 ! " ! ðmÞ ðmÞ ðmÞ ðmÞ pnþ1;m1 ðxÞ ¼ x aðmÞ pn;m1 ðxÞ bðmÞ n þ qnþ1 qn n qn an " ! " ðm1Þ ðm1Þ ðmÞ ðmÞ ðmÞ þ qðmÞ a ðxÞ þ b q q p b pn2;m1 ðxÞ: n1;m1 n1 n1 n1 n n n
374
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
Comparing with the recurrence relation for fpn;m1 ðxÞg we obtain formulas for the recursion coefficients. The case n ¼ 1 should be considered separately. Notice that in Chebyshev case ðm ¼ 0Þ we have anð0Þ ¼ 0
ð0Þ
ðn P 0Þ;
b0 , p;
1 ð0Þ b1 ¼ ; 2
bnð0Þ ¼
1 4
ðn P 2Þ:
Also, from (3.6) it follows r1 ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffi c 2 1;
r4 ¼
ðc2 þ ð1=2ÞÞðc2 1Þ ; etc: c3 þ ð3=2Þc
r2 ¼
c2 1 ; c
r3 ¼
cðc2 1Þ ; c2 þ ð1=2Þ
Now, using the previous theorem we give explicit expressions for recursion coefficients for some important special cases, where c ¼ cosh b. ð1Þ
Case m ¼ 1. Here we have q1 ¼ eb ; qnð1Þ ¼ ð1=2Þeb ðn P 2Þ, and the recursion coefficients ð1Þ
1 ð1Þ a1 ¼ eb ; anð1Þ ¼ 0 ðn P 2Þ; 2 p 1 1 ð1Þ ; b1 ¼ 1 e2b ; bnð1Þ ¼ ðn P 2Þ: , sinh b 2 4
a0 ¼ eb ; ð1Þ
b0
ð2Þ
b Case m ¼ 2. Here, q1 ¼ eb tanh b, qð2Þ ðn P 2Þ, and n ¼ ð1=2Þe ð2Þ
a0 ¼ ð2Þ
b0 ,
1 ; cosh b p cosh b ; sinh3 b
ð2Þ
a1 ¼ eb tanh b; ð2Þ
b1 ¼
anð2Þ ¼ 0
1 1 e2b tanh2 b; 2
ðn P 2Þ; ð2Þ
b2 ¼
1 1 þ e2b ; 4
and bnð2Þ ¼ ð1=4Þ ðn P 3Þ. ð3Þ
ð3Þ
Case m ¼ 3. Here, q1 ¼ sinh b tanh b=ð2 þ coshð2bÞÞ, q2 ¼ e2b cosh b, ¼ 12eb ðn P 3Þ, and 3 cosh b sinh b ð3Þ ð3Þ ; a1 ¼ e2b cosh b þ eb þ tanh b; a0 ¼ 2 þ coshð2bÞ 2 þ coshð2bÞ
qnð3Þ
1 ð3Þ a2 ¼ e3b ; anð3Þ ¼ 0 ðn P 3Þ; 2 4 p cosh2 b þ ð1=2Þ ð1 e2b Þ ð3Þ ð3Þ ; b b0 , ¼ ; 1 sinh5 b 2ð1 4e2b þ e4b Þ2
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378 ð3Þ
b2 ¼
1 1 þ 3e2b 3e4b e6b ; 4
bð3Þ n ¼
1 4
375
ðn P 3Þ:
As we can see the recurrence coefficients for polynomials pn;m ðxÞ reduce to the corresponding coefficients for Chebyshev polynomials for n P n0 ðn 2 NÞ. Precisely, calculations show that
mþ1 ð0Þ aðmÞ ¼ a ¼ 0; n P þ1 n n 2 and 1 ð0Þ bðmÞ n ¼ bn ¼ ; 4
nP
hmi 2
þ 1:
4. Numerical examples In order to illustrate the presented transformation method, we consider in this section a few numerical examples. All computations were done in Darithmetic on the WO R K S T A T I O N DI G I T A L UL T I M A T E AL P H A 533au2 (with machine precision 2:22 1016 ). Example 4.1. Consider integrals of the form Z þ1 2 sin 2t 1 e cos 2t 2 Im ðf ; bÞ ¼ dt 2 m 1 3 þ 2 cos 3t ðt þ b Þ
ðm P 1Þ:
The function f ðtÞ ¼
2 sin 2t 1 cos 2t e 3 þ 2 cos 3t
is ð2pÞ-periodic and its graph on the interval ½p; p is displayed in Fig. 1.
Fig. 1. The periodic function f ðtÞ (a). The function F ðxÞ obtained by transformation (2.6) (b).
376
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
Since f ðsÞ þ f ðsÞ ¼
2e cos 2s ; 3 þ 2 cos 3s
putting x ¼ cos s and using (2.6) we find 2
F ðxÞ ¼ F ðcos sÞ ¼ f ðsÞ þ f ðsÞ ¼
2e12x 3 6x þ 8x3
and according to (2.7), Z 1 pm ðxÞ dx Im ðf ; bÞ ¼ F ðxÞ m pffiffiffiffiffiffiffiffiffiffiffiffiffi ; ðc xÞ 1 x2 1 where c ¼ cosh b. Let m ¼ 1. Applying Gaussian quadratures with the ChW for n ¼ 5ð5Þ50 and taking b ¼ 10m ðm ¼ 2; 1; 0Þ we get approximations of I1 ðf ; bÞ with relative errors given in Table 1. Numbers in parentheses indicate decimal exponents. Taking Gaussian quadratures for n ¼ 5ð5Þ50, with SBW, v1 ðxÞ ¼ 1=2 ð1 x2 Þ =ðc xÞ, the corresponding errors are also presented in the same table. The corresponding exact values of I1 ðf ; bÞ are obtained using Gaussian quadratures with SBW in Q-arithmetic (machine precision 1:93 1034 ): I1 ðf ; 0:01Þ ¼ 0:2586588216241823127882 . . . 102 ðc ¼ 1:0000500 . . .Þ; I1 ðf ; 0:10Þ ¼ 0:4968012877996286228355 . . . 101 ðc ¼ 1:0050041 . . .Þ; I1 ðf ; 1:00Þ ¼ 0:1673215409745331112726 . . . 101 ðc ¼ 1:5430806 . . .Þ: Table 1 Relative errors in Gaussian approximations of the integral I1 ðf ; bÞ with respect to the Chebyshev weight (ChW) and the Szeg} o–Bernstein weight (SBW) b
b ¼ 0:01
n
ChW
SBW
ChW
SBW
ChW
SBW
5 10 15 20 25 30 35 40 45 50
8:4 ð1Þ 8:0 ð1Þ 7:6 ð1Þ 7:2 ð1Þ 6:7 ð1Þ 6:3 ð1Þ 5:9 ð1Þ 5:5 ð1Þ 5:2 ð1Þ 4:8 ð1Þ
1:3 ð2Þ 2:4 ð4Þ 1:1 ð5Þ 9:0 ð7Þ 1:6 ð8Þ 7:4 ð10Þ 5:9 ð11Þ 1:0 ð12Þ 4:8 ð14Þ 4:7 ð15Þ
2:2 ð1Þ 1:1 ð1Þ 4:3 ð2Þ 1:6 ð2Þ 6:0 ð3Þ 2:2 ð3Þ 8:2 ð4Þ 3:0 ð4Þ 1:1 ð4Þ 4:1 ð5Þ
6:3 ð2Þ 1:5 ð3Þ 4:2 ð5Þ 4:4 ð6Þ 9:7 ð8Þ 2:8 ð9Þ 2:9 ð10Þ 6:4 ð12Þ 1:8 ð13Þ 1:9 ð14Þ
1:2 ð2Þ 2:7 ð3Þ 1:5 ð4Þ 1:0 ð6Þ 1:8 ð7Þ 9:8 ð9Þ 6:7 ð11Þ 1:2 ð11Þ 6:5 ð13Þ 5:2 ð15Þ
5:5 ð2Þ 3:5 ð3Þ 7:0 ð5Þ 3:5 ð6Þ 2:3 ð7Þ 4:6 ð9Þ 2:3 ð10Þ 1:5 ð11Þ 3:1 ð13Þ 1:6 ð14Þ
b ¼ 0:1
b ¼ 1:0
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
377
As can be seen, for smaller values of b (c is close to 1) the Gauss–Chebyshev quadratures ChW cannot be used directly. When b increases the both quadratures become comparable. However, by writing I1 ðf ; bÞ in the form Z p sinh b 1 F ðcÞ F ðxÞ dx pffiffiffiffiffiffiffiffiffiffiffiffiffi ; I1 ðf ; bÞ ¼ F ðcÞ 2b 2b c x 1 x2 1 the Gauss–Chebyshev quadratures can be applied directly. Consider now the case m ¼ 2, with the functions /k and the corresponding weights vk ðk ¼ 0; 1; 2Þ, where /k ðxÞ ¼
F ðxÞp2 ðxÞ ðc xÞ
2k
;
vk ðxÞ ¼
1 pffiffiffiffiffiffiffiffiffiffiffiffiffi : ðc xÞ 1 x2 k
Applying the Gaussian quadratures with the Chebyshev weight ChW (k ¼ 0) and the Szeg} o–Bernstein weights SBW1 (k ¼ 1) and SBW2 (k ¼ 2) we get approximations of the integral I2 ðf ; bÞ. The exact values of this integral for some selected b are: I2 ðf ; 0:01Þ ¼ 0:1156183821140487028202 . . . 106 ; I2 ðf ; 0:10Þ ¼ 0:1214706913588412300593 . . . 103 : The relative errors in Gaussian approximations for n ¼ 5ð5Þ50 are presented in Table 2. The advantage of quadrature formulas for k ¼ 2 (in this case m ¼ 2) is evident. When b increases all quadratures give similar results. Example 4.2. Consider now the integral (1.4), with a nonanalytic function 2a f ðtÞ ¼ j cosðt=2Þj (a > 0). After the transformation we obtain the integral Table 2 Relative errors in Gaussian approximations of the integral I2 ðf ; bÞ with respect to the Chebyshev weight (ChW) and to the Szeg} o–Bernstein weights (SBW1 and SBW2 ) b
b ¼ 0:01
n
ChW
SBW1
SBW2
ChW
SBW1
SBW2
5 10 15 20 25 30 35 40 45 50
1:0 ð0Þ 1:0 ð0Þ 1:0 ð0Þ 9:9 ð1Þ 9:9 ð1Þ 9:8 ð1Þ 9:7 ð1Þ 9:6 ð1Þ 9:5 ð1Þ 9:3 ð1Þ
9:1 ð1Þ 8:3 ð1Þ 7:5 ð1Þ 6:8 ð1Þ 6:1 ð1Þ 5:5 ð1Þ 5:0 ð1Þ 4:5 ð1Þ 4:1 ð1Þ 3:7 ð1Þ
5:5 ð7Þ 1:0 ð7Þ 4:7 ð9Þ 1:9 ð11Þ 4:6 ð12Þ 2:6 ð12Þ 2:3 ð12Þ 2:3 ð12Þ 2:3 ð12Þ 2:3 ð12Þ
8:9 ð1Þ 6:2 ð1Þ 3:4 ð1Þ 1:6 ð1Þ 7:4 ð2Þ 3:2 ð2Þ 1:3 ð2Þ 5:6 ð3Þ 2:3 ð3Þ 9:2 ð4Þ
3:7 ð1Þ 1:4 ð1Þ 5:0 ð2Þ 1:9 ð2Þ 6:8 ð3Þ 2:5 ð3Þ 9:2 ð4Þ 3:4 ð4Þ 1:2 ð4Þ 4:6 ð5Þ
1:1 ð3Þ 6:7 ð5Þ 4:3 ð6Þ 6:4 ð8Þ 4:4 ð9Þ 2:8 ð10Þ 4:3 ð12Þ 3:1 ð13Þ 1:4 ð15Þ 2:0 ð14Þ
b ¼ 0:1
378
G. Mastroianni, G.V. Milovanovic / Appl. Math. Comput. 128 (2002) 365–378
Fig. 2. Relative errors in Gaussian approximations with n ¼ 5 (upper curve) and n ¼ 20 nodes (lower curve) for 0 6 a < 4:5.
J ðaÞ ¼ A
Z
1
1
1þx 2
a
dx pffiffiffiffiffiffiffiffiffiffiffiffiffi ; ðc xÞ 1 x2
where A ¼ 2p1 ðxÞ ¼ sinh b=b. In order to evaluate this integral, we apply Gaussian rule in n points with SBW ðm ¼ 1Þ. A typical behavior of the relative error rn of Gaussian approximations with respect to the parameter a ð0 6 a < 4:5Þ is displayed in Fig. 2 in the log-scale. Two cases for n ¼ 5 and n ¼ 20 are given, whereas b ¼ 0:01. It is clear that the rapidly increasing of accuracy achieves when the parameter a tends to an integer (i.e., when f becomes an analytic function).
References [1] B. Fischer, G. Golub, How to generate unknown orthogonal polynomials out of known orthogonal polynomials, J. Comput. Appl. Math. 43 (1992) 99–115. [2] W. Gautschi, A survey of Gauss–Christoffel quadrature formulae, in: P.L. Butzer, F. Feher (Eds.), E.B. Christoffel, Birkh€auser, Basel, 1981, pp. 72–147. [3] W. Gautschi, Orthogonal polynomials: applications and computation, Acta Numerica (1996) 45–119. [4] G. Golub, J.H. Welsch, Calculation of Gauss quadrature rules, Math. Comp. 44 (1969) 221– 230. [5] G. Mastroianni, Generalized Christoffel functions and error of positive quadrature, Numer. Algorithms 10 (1995) 113–126. [6] A.P. Prudnikov, Yu.A. Brychkov, O.I. Marichev, Integrals and Series. Elementary Functions, Nauka, Moscow, 1981 (in Russian). [7] G. Szeg} o, Orthogonal Polynomials, Amer. Math. Soc. Colloq. Publ., vol. 23, fourth ed., Amer. Math. Soc., Providence, RI, 1975.