The existence and construction of rational Gauss-type quadrature rules

Applied Mathematics and Computation 218 (2012) 10299–10320

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

The existence and construction of rational Gauss-type quadrature rules q Karl Deckers ⇑,1, Adhemar Bultheel Department of Computer Science, K.U. Leuven, Heverlee, Belgium

a r t i c l e

i n f o

Keywords: Quasi-orthogonal rational functions Generalized eigenvalue problem Positive rational interpolatory quadrature rules

a b s t r a c t Consider a hermitian positive-definite linear functional F, and assume we have m distinct nodes fixed in advance anywhere on the real line. In this paper we then study the existence and construction of nth rational Gauss–Radau ðm ¼ 1Þ and Gauss–Lobatto ðm ¼ 2Þ quadrature formulas that approximate Fff g. These are quadrature formulas with n positive weights and with the n  m remaining nodes real and distinct, so that the quadrature is exact in a ð2n  mÞ-dimensional space of rational functions. Further, we also consider the case in which the functional is defined by a positive bounded Borel measure on an interval, for which it is required in addition that the nodes are all in the support of the measure. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction The central object of study in this paper is the numerical approximation of a hermitian positive-definite linear functional F by means of an n-point positive interpolatory quadrature rule with real distinct nodes fxk gnk¼1 and positive weights fkk gnk¼1 :

Fff g 

n X kk f ðxk Þ: k¼1

When taking as nodes the zeros of an nth orthogonal polynomial (OP) (where orthogonality is with respect to the inner product defined by hf ; gi ¼ Fff ðxÞ  gðxÞg), we obtain the nth Gaussian quadrature rule, which has a maximal polynomial domain of validity; i.e., the approximation is exact for every polynomial of degree less than or equal to 2n  1. Clearly, the existence of the nth Gaussian quadrature rule depends on whether an nth OP exists. If one or two nodes in the quadrature formula are fixed in advance, and the remaining nodes are chosen to obtain exactness for at least every polynomial of degree less than or equal to 2n  2 or 2n  3 respectively, we get more general nth Gauss-type quadrature rules. In this case, the nodes are zeros of an nth quasi-orthogonal polynomial, and the existence depends on the fixed node(s) as well (i.e., for certain choices, one of the remaining nodes may lie at infinity and/or some of the corresponding weights may be non-positive). See e.g. [1,14,16,17] and the references therein. When f has singularities, it is often more appropriate to not consider a maximal polynomial domain of validity, but rather consider more general spaces of rational functions. In such a case the OPs are replaced by orthogonal rational functions (ORFs) with preassigned poles (to simulate the singularities of f).

q This work is partially supported by the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian State, Science Policy Office. The scientific responsibility rests with its authors. ⇑ Corresponding author. E-mail addresses: [email protected] (K. Deckers), [email protected] (A. Bultheel). 1 The first author is a Postdoctoral Fellow of the Research Foundation–Flanders (FWO).

0096-3003/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.04.008

10300

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

ORFs were first introduced by Dzˇrbašian in the 1960s. Most of his papers appeared in Russian literature, but an accessible survey in English can be found in [11,18]. These ORFs are a generalization of OPs in such a way that they are of increasing degree with a given sequence of poles, and the OPs result if all the poles are at infinity. So, when taking as nodes the zeros of an nth ORF instead, we obtain the nth rational Gaussian quadrature rule; see e.g. [2, Chapter 11.6] and [3,9,10,12,13,15,19,20,22–25]. Here we can also consider more general rational Gauss-type quadrature rules obtained by fixing one or two nodes in the quadrature rule. However, the existence of these quadrature rules now also depends on the sequence of poles. The aim of this paper is to investigate the existence and construction of rational Gauss-type quadrature rules. As opposed to the polynomial case, in the rational case we can distinguish two different kinds of positive rational interpolatory quadrature rules for a given set of m > 0 nodes fixed in advance. The difference lies in the domain of validity of the quadrature rule. Only when the domain of validity is fixed, the quadrature rule is (as in the polynomial case) unique whenever it exists. So far, only one kind of positive rational interpolatory quadrature rule with m ¼ 1 has been investigated in detail; see [7,8]. In this paper we will investigate the other kind for m ¼ 1, as well as one quadrature rule for the case in which m ¼ 2. Furthermore, we will provide a framework to deal with the general case m > 0 for both kinds of quadrature rules. The outline of the paper is as follows. After giving the necessary theoretical background in Section 2, in Section 3 we discuss the connection of positive rational interpolatory quadrature formulas with certain kinds of quasi-orthogonal rational functions. Next, in Section 4 we use this connection to study the existence and construction of rational Gauss-type quadrature formulas for a given (set of) node(s) fixed in advance. In special, we also consider the case in which the functional F is defined by a positive bounded Borel measure on an interval. Finally, some illustrative numerical examples are given in Section 5. 2. Preliminaries The field of complex numbers will be denoted by C and the Riemann sphere by C ¼ C [ f1g. For the real line we use the symbol R and for the extended real line R ¼ R [ f1g. Further, the positive half line will be represented by Rþ ¼ fx 2 R : x P 0g. Whenever the value zero is omitted in a set X # C, this will be represented by X 0 . Similarly, the complement of a set Y  C with respect to a set X # C will be given by X Y ; i.e., X Y ¼ ft 2 X : t R Yg. Let a 2 C, then Refag refers to the real part of a, while I mfag refers to the imaginary part, and the imaginary unit will be denoted by i. For a – 0 we will also use the short notation T fag, defined by

( Imfag jaj2

Tfag ¼ I mf1=ag ¼

0;

; a 2 C0 ; a ¼ 1:

For any complex function f, we define the involution operation or super-c conjugate by f c ðxÞ ¼ f ðxÞ. With P k we denote the space of polynomials of degree less than or equal to k, while P represents the space of all polynomials. Further, we will use the short notation ½f ðxÞx¼a to denote limx!a f ðxÞ. Let a sequence of poles A ¼ fa1 ; a2 ; . . .g  C0 be fixed, where the poles are arbitrary complex or infinite; hence, they do not have to appear in pairs of complex conjugates. The rational functions we then deal with, are of the form

fk ðxÞ ¼

ck xk þ ck1 xk1 þ    þ c0 x0 ; ð1  x=a1 Þð1  x=a2 Þ    ð1  x=ak Þ

k ¼ 1; 2; . . .

Define the factors

Z k ðxÞ :¼ Z ak ðxÞ ¼

x ; 1  x=ak

k ¼ 1; 2; . . .

ð1Þ

and the basis functions

b0 ðxÞ  1;

bk ðxÞ ¼ bk1 ðxÞZ k ðxÞ;

k ¼ 1; 2; . . .

ð2Þ

These basis functions generate the nested spaces of rational functions with poles in A defined by

L1 ¼ f0g;

Lk :¼ Lfa1 ; . . . ; ak g ¼ spanfb0 ; . . . ; bk g;

L0 ¼ C;

With L we then denote the closed linear span of all

p0 ðxÞ  1; pk ðxÞ ¼

k Y

ð1  x=aj Þ;

k ¼ 1; 2; . . . ;

j¼1

then for k P 0 we may write equivalently

bk ðxÞ ¼

xk pk ðxÞ

fbk g1 k¼0 .

and Lk ¼ fpk =pk : pk 2 P k g:

Let

k ¼ 1; 2; . . .

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

10301

½a

In the remainder we will also use the short notation Lk :¼ Lfa1 ; . . . ; ak1 ; ag (i.e., the space of rational functions with the same dimension as Lk , but with the last pole ak replaced by the pole a) and Lk ðaÞ :¼ ff ¼ ppk 2 Lk : pk ðaÞ ¼ 0g ¼ ZZak Lk1 . k

Note that Lk and L are rational generalizations of P k and P. Indeed, if aj ¼ 1 for every j P 1, the expression in (1) becomes Z k ðxÞ ¼ x and the expression in (2) becomes bk ðxÞ ¼ xk . Thus the polynomial case is automatically accounted for. With the definition of the super-c conjugate we introduce Lck ¼ ff c : f 2 Lk g. Further, with Rk;j and S k ¼ Sfa1 ; . . . ; ak g we denote the space of rational functions defined by Rk;j ¼ Lk  Lcj :¼ ff  g c : f 2 Lk and g 2 Lj g and S k ¼ Lk þ Lck :¼ ff þ g c : f 2 Lk and g 2 Lk g respectively. Consider an inner product defined by a linear functional F: c

hf ; giF ¼ Fffg g;

f ; g 2 L;

where the linear functional F is assumed to be hermitian positive-definite (HPD); i.e., c

c

Fffg g ¼ Fff c gg and Ffff g > 0 for f – 0: Orthogonalizing the basis functions fb0 ; b1 ; . . .g with respect to this inner product, we obtain a sequence of orthogonal rational functions (ORFs) f/0 ; /1 ; . . .g, with /k 2 Lk n Lk1 , so that /k ?F Lk1 ; i.e.,

1

h/k ; /j iF ¼

dk 2 C0 ;

dk;j ;

jdk j2

k; j ¼ 0; 1; . . . ;

where dk;j is the Kronecker Delta. In the special case in which uk ðxÞ ¼ dk /k ðxÞ, so that

kuk kF :¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi huk ; uk iF ¼ 1;

we say that uk is an orthonormal rational function (nORF). Put by convention a1 ¼ a0 ¼ 1. We then call a rational function fk ¼ ppk 2 Lk n Lk1 , with k P 0, exceptional (respectively k degenerate) iff pk ðak1 Þ ¼ 0 (respectively pk ðak1 Þ ¼ 0). A zero of pk at 1 means that the degree of pk is less than k. Further, for k > 1 we say that fk is singular iff pk ðaj Þ ¼ 0 for at least one j 2 f0; . . . ; k  2g. Finally, fk is called regular iff fk is not singular, not degenerate and not exceptional. With these definitions we now introduce the following spaces of rational functions: e

Lk :¼ ff 2 Lk n Lk1 : f is not exceptionalg;

d

Lk :¼ ff 2 Lk n Lk1 : f is not degenerateg;

de r

Lk :¼ e Lk \ d Lk ;

Lk :¼ ff 2 Lk n Lk1 : f is regularg:

In [23, Thm. 2.1.1] (and [2, Chapter 11.1] for the special case of all real poles) the following three-term recurrence relation has been proved for nORFs uk 2 Lk n Lk1 . Theorem 1. Consider the nORFs uj 2 Lj n Lj1 , with j ¼ k; k  1; k  2 and k > 0. Then these nORFs satisfy a three-term recurrence relation of the form



uk ðxÞ ¼ Ek Z k ðxÞ

1þ

  Dk C uk1 ðxÞ  c k uk2 ðxÞ ¼ Ek wk ðxÞ; Z k1 ðxÞ Z k2 ðxÞ

ð3Þ

þ e d with jEk j ¼ kwk k1 F 2 R0 ; C k 2 C0 and Dk 2 C iff uk 2 Lk and uk1 2 Lk1 . The initial conditions are a1 ¼ a0 ¼ 1; u1 ðxÞ  0, 1 de and ju0 ðxÞj  k1kF ¼ jE0 j. In the special case in which uk1 2 Lk1 and uk2 2 d Lk2 for k > 1, it holds that (see [5, Thm. 3.3])

Ck ¼

1  Dk =Z ck1 ðak1 Þ Ek1

2 C0 :

In this case the coefficient Dk can be expressed in terms of inner products as follows (see [5, Thms. 3.7 and 3.9]):

Dk ¼

K k;j  Lk;j Lk;j Z k1 ðak Þ

þ Zc

K k;j ðak1 Þ

k1

þ dk1;j Ek1

;

0 6 j < k;

with

K k;j ¼

1 Z ck2 ð k Þ

a

D E Z k uk2 ; uj þ dk2;j F

D E and Lk;j ¼ Ek1 Z k uk1 ; uj ; F

while for ak R R, the modulus of the coefficient Ek can also be computed as follows:

jEk j2 ¼

Tfak gjEk1 j2   : ðI mfDk g  jDk j2 Tfak1 gÞ  jEk1 j2  4Tfak1 g  Tfak2 g þ Tfak2 g

In the opposite direction as in Theorem 1, the following Favard-type theorem has been proved in [4].

10302

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

Theorem 2. Let fvk g1 k¼0 be a sequence of rational functions in L, and assume that the following conditions are satisfied: (A1) a1 ¼ a0 ¼ 1 and ak 2 C0 ; k ¼ 1; 2; . . .; (A2) vk is generated by the three-term recurrence relation (3); (A3) vk 2 Lk n Lk1 ; k ¼ 0; 1; 2; . . ., and v1  0; (A4) Dk 2 C and Ek 2 C0 ; k ¼ 1; 2; . . .; (A5)

ImfDk g ¼

Tfak g jEk j

2



Tfak2 g jEk1 j2

if ak1 2 R0 , respectively 2

RefDk g2 þ fI mfDk g  iZ k1 ðak1 Þg ¼ fZ k1 ðak1 Þg2

jEk1 j2 jEk j2



Dk Dk1

if ak1 R R, wherec Dk ¼ jEk j2  4Tfak g  Tfak1 g > 0; k ¼ 1; 2; . . ., with E0 2 C0 ; 1Dk =Z k1 ðak1 Þ (A6) C k ¼ 2 C0 ; k ¼ 1; 2; . . .Then there exists a HPD linear functional G on L  Lc so that E k1

c

hf ; giG ¼ Gffg g defines an inner product on L for which the

vk form an orthonormal system.

The above three-term recurrence relation can also be written as follows:

(

)





½u ak1 x x x u ½u u x uk1 ðxÞ  c uk2 ðxÞ ¼ ak1 1 uk2 ðxÞ þ b½k1 1 uk1 ðxÞ þ c½k1 1 u ðxÞ; ak2 ak1 ak k Z k2 ðak2 Þ

0 < k 6 n;

where ½u bk1

( ½u ck1

¼ Dk ;

¼

E1 k ;

½u ak1

and

¼ Ck ¼

) ½u  bk1 ½u 1þ c c : Z k1 ðak1 Þ k2

ð4Þ

Let In denote the n by n identity matrix, and define the matrices

0

½u 

b0

B B ½u Ba B 1 B Jn ¼ B B 0 B B .. B . @ 0 0 0 B a½u B c1 B Z 0 ða 0 Þ B B Sn ¼ B B 0 B . B . B . @ 0

½u

c0

0

½u

b1 ..

.

..

.

... .. . ..

.

..

.

...

1

C C C C C  C; Dn ¼ diag a1 ; a1 ; . . . ; a1 ; 0 1 n1 . 0 C C C ½u ½u  ½u  an2 bn2 cn2 C A ½u ½u  0 an1 bn1 1 ... ... 0 .. C .C C C .. .. C . .C C; and Bn ¼ Jn Dn þ In  Sn ; .. .. .. C C . . .C A ½u  a 0 0 Zc n1 ða Þ ½u

c1 ..

...

0 .. .

... .. . .. .

n2

n2

and column vectors

~ n ðxÞ ¼ ðu0 ðxÞu1 ðxÞ . . . un1 ðxÞÞT and ~ u en ¼ ð0    01ÞT 2 Cn : Assuming the nORFs uk 2 de Lk for k ¼ 0; . . . ; n  1, and un 2 e Ln , we obtain that



x u ~ n ðxÞ ¼ xBn u ~ n ðxÞ  c½n1 Jn u 1 un ðxÞ~en :

an

ð5Þ

The following theorem has then be proved in [21, Section. 4]. Theorem 3. Suppose the nORFs uk 2 de Lk for k ¼ 0; . . . ; n  1. Then the zeros xn;j ; j ¼ 1; . . . ; n, of a nORF un ðxÞ 2 e Ln are eigenvalues of the generalized eigenvalue problem (GEP)

10303

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

Jn~ v n;j ¼ xn;j Bn~ v n;j ;

ð6Þ

with

~ v n;j

( )1=2 n1 X 2 ~ n ðxn;j Þ; ¼ gn juk ðxn;j Þj u

jgn j ¼ 1;

k¼0

the corresponding normalized eigenvector. In the remainder of this paper we will assume that the sequence of nORFs fuk gnk¼0 exists and that uk 2 de Lk for k ¼ 0; . . . ; n.2 3. Positive rational interpolatory quadrature rules Consider the set of n nodes fxn;k gnk¼1  RAn , where An ¼ fak gnk¼1 and xn;k – xn;j if k – j, and weights fkn;k gnk¼1  Rþ 0 . We then call a quadrature formula of the form F

Fn ff g :¼

F

F

F

n X F F kn;k f ðxn;k Þ;

ð7Þ

k¼1

an nth positive rational interpolatory quadrature formula (abbreviated PRIQ) iff Rn ff g :¼ Fff g  Fn ff g ¼ 0 for every e c , where L0 # L e j  L½an  . If R ~ is the largest space for which Rn ff g ¼ 0, then it is called the domain f 2 R ~ ¼ Ln1  L nþ1

j

n1;j

n1;j

of validity. Clearly, a necessary (but not sufficient) condition for the existence of a PRIQ is that Rn1;~j ¼ Rcn1;~j ; i.e., for every function f in the domain of validity it should hold that f c is in the domain of validity too.3 With ‘ðn þ i  1 : n  i  m; An ; FÞ-PRIQ’, where i 2 f0; 1g and m 2 f1  i; . . . ; n  ig, we denote an nth PRIQ for which the space Rnþi1;nim is contained in the domain of validity (i.e., the quadrature is exact for at least every f 2 Rnþi1;nim ). To study ðn þ i  1 : n  i  m; An ; FÞ-PRIQs we will need the so-called quasi-orthogonal rational functions (qORFs), which are defined as follows. Definition 1. Let 0 6 m 6 n. We then call a rational function Q n;~amþ1 2 Ln an m-qORF of the first kind iff Q n;~amþ1 is of the form

Q n;~amþ1 ðxÞ ¼

m X aj unj ðxÞ;

~ amþ1 ¼ ða0 . . . am ÞT 2 Cmþ1 :

j¼0

Definition 2. Let 0 6 m 6 n. We then call a rational function P n;~amþ1 2 Ln an m-qORF of the second kind iff P n;~amþ1 is of the form

Pn;~amþ1 ðxÞ ¼ a0 un ðxÞ þ

m X Z n ðxÞ aj c u ðxÞ; Z nj ðxÞ nj j¼1

~ amþ1 ¼ ða0 . . . am ÞT 2 Cmþ1 :

From the previous definitions it clearly follows that Q n;~amþ1 ?F Lnðmþ1Þ and that

8 m ¼ 0; > < Ln1 ; Z nm L ð a Þ  L ; 0 < m < n; Pn;~amþ1 ?F nm n nðmþ1Þ Z cn > : m ¼ n: L0 ðan Þ  L1 ; Conversely, for every r n 2 Ln and every m 2 f0; . . . ; ng it holds that rn ?F Lnðmþ1Þ iff rn is an m-qORF of the first kind. Further, ½a

½a

½a

assume that there exists a rational function un nm 2 e L½nanm  so that un nm ?F Ln1 and kun nm kF ¼ 1. Then it holds that 





r n ?F Lnm ðan Þ iff rn is an m-qORF of the second kind. The first statement is trivial. To prove the second statement, we first need the following lemma. mþ1 Lemma 4. Let 1 6 m 6 n. Then for every sequence of constants faj gm j¼1  C there exists a sequence of constants fbj gj¼0  C so that

Z n ðxÞ

m m þ1 X X Z n ðxÞ aj unj ðxÞ  b0 un ðxÞ þ bj c unj ðxÞ: Z nj ðxÞ j¼1 j¼1

n o ½u ½u n2 ½u n1 It is easily proved by induction over k ¼ 1; . . . ; n, with the aid of (4), (5) and Theorem 1, that uk 2 de Lk for k ¼ 0; . . . ; n iff akþ1 ; ck  C0 ; fbk gk¼0 C n o k¼0 1 and det ak Jk  Bk – 0 for k ¼ 1; . . . ; n. 2

3

If this condition is not satisfied, the weights cannot be all real; see [6, Sect. 2].

10304

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

Proof. For m ¼ 1 it follows from Theorem 1 that

  1 Dn Z ðxÞ un ðxÞ ¼ Z n ðxÞ 1 þ u ðxÞ  C n c n u ðxÞ; En Z n1 ðxÞ n1 Z n2 ðxÞ n2 where it holds that

1þ

Dn Dn : ¼ En1 C n þ c Z n1 ðxÞ Z n1 ðxÞ

Consequently, the statement follows for m ¼ 1 by taking b0 ¼ E

a1

n En1 C n

2 C; b1 ¼  Ea1 DCn 2 C and b2 ¼ Ea1 2 C. n1 n

n1

So, suppose now that the statement holds true for 1 6 m  1 < n. For m we then have that

Z n ðxÞ

m m1 m X X X Z n ðxÞ aj unj ðxÞ ¼ Z n ðxÞ aj unj ðxÞ þ am Z n ðxÞunm ðxÞ  b0 un ðxÞ þ bj c unj ðxÞ þ am Z n ðxÞunm ðxÞ: Z nj ðxÞ j¼1 j¼1 j¼1

From Theorem 1 it now follows that

1 Enðm1Þ

  D Z n ðxÞ Z ðxÞ unðm1Þ ðxÞ ¼ Z n ðxÞ 1 þ nðm1Þ unm ðxÞ  C nðm1Þ c n u ðxÞ; Z nðm1Þ ðxÞ Z nm ðxÞ Z nðmþ1Þ ðxÞ nðmþ1Þ

where it holds that

1þ

Dnðm1Þ Dnðm1Þ ¼ Enm C nðm1Þ þ c Z nm ðxÞ Z nm ðxÞ

1 1 1  : ¼ Z nðm1Þ ðxÞ Z cnðm1Þ ðxÞ Z cnðm1Þ ðanðm1Þ Þ

and

As a result we find that

am Z n ðxÞunm ðxÞ ¼

m þ1 X j¼m1

with dm1 ¼ E

am nðm1Þ Enm C nðm1Þ

dj

Z n ðxÞ u ðxÞ þ em Z n ðxÞunðm1Þ ðxÞ; Z cnj ðxÞ nj am Dnðm1Þ

2 C; dm ¼  E

nm C nðm1Þ

2 C; dmþ1 ¼ Eam 2 C and em ¼  E nm

from the induction hypothesis it follows that there exist constants

Z n ðxÞunðm1Þ ðxÞ  ^c0 un ðxÞ þ

m X

^cj

j¼1

f^cj gm j¼0

am c nðm1Þ Enm C nðm1Þ Z nðm1Þ ðanðm1Þ Þ

2 C. Finally,

 C so that

Z n ðxÞ u ðxÞ: Z cnj ðxÞ nj

This ends the proof. h We are now able to prove the following. Theorem 5. Let 1 6 m 6 n, and suppose rn 2 Ln . Further, assume that

un½anm  2 e Ln½anm  : un½anm  ?F Ln1 and kun½anm  kF ¼ 1:

ð8Þ

Then it holds that rn ?F Lnm ðan Þ iff rn is an m-qORF of the second kind. Proof. From the definition of an m-qORF of the second kind it is clear that rn ?F Lnm ðan Þ if rn is an m-qORF of the second kind. So, we only need to prove that r n is an m-qORF of the second kind if r n ?F Lnm ðan Þ. The case in which m ¼ 1 has already been proved in [8, Thm. 5] under the less restrictive condition that

u½nan1  2 L½nan1  n Ln1 . Thus, suppose that 1 < m 6 n. Since rn ?F Lnm ðan Þ, it follows that Z cnm Z n r n ?F Lnðmþ1Þ .

½anm  n

Consider now the nORF u

2

e ½anm  Ln .

Then there exist constants

m X Z cnm ðxÞ bj unj ðxÞ r n ðxÞ ¼ b0 un½anm  ðxÞ þ Z n ðxÞ j¼1

and hence,

rn ðxÞ ¼ b0 ¼ b0

m X Z n ðxÞ ½anm  Z ðxÞ un ðxÞ þ bj c n u ðxÞ c Z nm ðxÞ Z nm ðxÞ nj j¼1 m1 m X X Z n ðxÞ ½anm  bj Z ðxÞ un ðxÞ þ Z n ðxÞ unj ðxÞ þ bj c n unj ðxÞ: c c Z ð a Þ Z nm ðxÞ Z nj nj ðxÞ j¼1 nm j¼1

fbj gm j¼0

Z cnm Zn rn

 C so that

½a

2 Ln nm



and

10305

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

From the previous lemma it now follows that there exist constants fcj gm j¼0  C so that

rn ðxÞ ¼ b0

m X Z n ðxÞ ½anm  Z ðxÞ un ðxÞ þ c0 un ðxÞ þ cj c n u ðxÞ: c Z nm ðxÞ Z nj ðxÞ nj j¼1 ½a

½a

Thus, it remains to prove that ZcZn un nm is an m-qORF of the second kind. Since un nm 2 e L½nanm  , we deduce from Theorem 1 nm that

1 E½nanm 





" # Z n ðxÞ ½anm  D½nanm  Z ðxÞ u ðxÞ ¼ Z n ðxÞ 1 þ u ðxÞ  C n½anm  c n u ðxÞ Z n1 ðxÞ n1 Z cnm ðxÞ n Z n2 ðxÞ n2 ¼ En1 C n½anm  Z n ðxÞun1 ðxÞ þ Dn½anm 

Z n ðxÞ Z ðxÞ u ðxÞ  C ½nanm  c n u ðxÞ: Z cn1 ðxÞ n1 Z n2 ðxÞ n2 ½a



Finally, because Z n un1 is a 2-qORF of the second kind due to the previous lemma, it follows that ZcZn un nm is an m-qORF of nm the second kind. This concludes the proof. h ½a



In everything that follows, we will assume that the nORF un nm in (8) exists. The following two theorems now give the connection between m-qORFs and ðn : n  1  m; An ; FÞ-PRIQs, where i 2 f0; 1g and m 2 f1  i; . . . ; n  ig. Theorem 6. Let 0 6 m 6 n  1, and consider a ðn : n  1  m; An ; FÞ-PRIQ of the form (7). Then there exists an m-qORF of the F first kind Q n;~amþ1 2 r Ln so that Q n;~amþ1 ðxn;k Þ ¼ 0 for k ¼ 1; . . . ; n, and ½Q

F

kn;k ¼ FfLn;k g;

k ¼ 1; . . . ; n;

ð9Þ

½Q 

with Ln;k the fundamental interpolating rational functions defined by ½Q 

Ln;k ðxÞ ¼

Q n;~amþ1 ðxÞ ð1  x=an Þ  2 Ln1 ; F F F ð1  xn;k =an Þ ðx  xn;k ÞQ 0n;~amþ1 ðxn;k Þ

k ¼ 1; . . . ; n;

where the prime denotes the derivative with respect to x. ½Q 

Proof. (The proof is similar to the proof of [2, Thm. 11.6.1].) First, note that Ln;k ðxn;j Þ ¼ dk;j ; 1 6 k; j 6 n. Consider now the auxiliary function

hðxÞ ¼ f ðxÞ 

n X

½Q 

f ðxn;k Þ  Ln;k ðxÞ:

ð10Þ

k¼1

For f 2 Rn;n1m , this can be written as

hðxÞ ¼

p2n1m ðxÞ

pn ðxÞ  pcn1m ðxÞ

;

p2n1m 2 P 2n1m :

Since hðxn;k Þ ¼ 0 for k ¼ 1; . . . ; n (interpolation property), and since the numerator of Q n;~amþ1 is a nonzero constant times the nodal polynomial ðx  xn;1 Þ . . . ðx  xn;n Þ, we can write

hðxÞ ¼ Q n;~amþ1 ðxÞ  g c ðxÞ;

g 2 Ln1m :

From (10), together with (7) and (9), it now follows that

Ffhg ¼ FfQ n;~amþ1  g c g ¼ Rn ff g: Finally, since the quadrature is exact for every f 2 Rn;n1m (i.e., Rn ff g ¼ 0 for every f 2 Rn;n1m ), it follows that Q n;~amþ1 is an m-qORF of the first kind. This ends the proof. h Theorem 7. Let 1 6 m 6 n, and consider a ðn  1 : n  m; An ; FÞ-PRIQ of the form (7). Then there exists an m-qORF of the secF ond kind Pn;~amþ1 2 r Ln so that Pn;~amþ1 ðxn;k Þ ¼ 0 for k ¼ 1; . . . ; n, and ½P

F

kn;k ¼ FfLn;k g; with

½P Ln;k

k ¼ 1; . . . ; n;

ð11Þ

the fundamental interpolating rational functions defined by

½P

Ln;k ðxÞ ¼

Pn;~amþ1 ðxÞ ð1  x=an Þ  2 Ln1 ; F F F ð1  xn;k =an Þ ðx  xn;k ÞP0n;~amþ1 ðxn;k Þ

where the prime denotes the derivative with respect to x.

k ¼ 1; . . . ; n;

10306

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320 ½P

Proof. (The proof is similar as above.) First, note that Ln;k ðxn;j Þ ¼ dk;j ; 1 6 k; j 6 n. Consider the auxiliary function

hðxÞ ¼ f ðxÞ 

n X ½P f ðxn;k Þ  Ln;k ðxÞ:

ð12Þ

k¼1

For f 2 Rn1;nm , this can be written as

hðxÞ ¼

p2n1m ðxÞ

pn1 ðxÞ  pcnm ðxÞ

;

p2n1m 2 P 2n1m :

Since hðxn;k Þ ¼ 0 for k ¼ 1; . . . ; n (interpolation property), and since the numerator of Pn;~amþ1 is a nonzero constant times the nodal polynomial ðx  xn;1 Þ    ðx  xn;n Þ, we can write

hðxÞ ¼ Pn;~amþ1 ðxÞ  g c ðxÞ; where

gðxÞ ¼

ð1  x=an Þpn1m ðxÞ ; pnm ðxÞ

pn1m 2 P n1m :

Hence, g 2 Lnm ðan Þ. From (12), together with (7) and (11), it follows that

Ffhg ¼ FfPn;~amþ1  g c g ¼ Rn ff g: Finally, since the quadrature is exact for every f 2 Rn1;nm (i.e., Rn ff g ¼ 0 for every f 2 Rn1;nm ), it follows that Pn;~amþ1 is an m-qORF of the second kind. This ends the proof. h 4. Existence of rational Gauss-type quadrature rules The rational Gaussian quadrature rule is a special kind of PRIQ that has maximal domain of validity Rn;n1 ; hence, it is a ðn þ i  1 : n  i  m; An ; FÞ-PRIQ with i ¼ 1 and m ¼ 0. It is well known that the rational Gaussian quadrature rule exists iff an 2 R0 and un 2 r Ln . Whenever it exists, it is unique, and the nodes can be computed by solving the GEP (6), while the corresponding weights can be found by means of the corresponding eigenvectors as follows:

kn;k ¼ jv n;k j2 k1k2F ; F

ð13Þ

where v n;k denotes the first entry in the normalized eigenvector ~ v n;k . See e.g. [8, Sect. 5]. From the previous two theorems we can deduce that m-qORFs, with 1 6 m 6 n  i and i 2 f0; 1g, can be used to construct ðn þ i  1 : n  i  m; An ; FÞ-PRIQs with m nodes fixed in advance. Indeed, consider an nth PRIQ of the form

Fn ff g ¼

m n     X X F F F F kn;k f xn;k þ kn;k f xn;k ; k¼1

ð14Þ

k¼mþ1

where the sequence of m distinct nodes X m ¼ fxn;k gm k¼1  RAn are fixed in advance, and the remaining nodes and corresponding weights are chosen such that the PRIQ is a ðn þ i  1 : n  i  m; An ; FÞ-PRIQ. Then these remaining nodes and corresponding weights can be found by means of the m-qORF Q n;~amþ1 (for i ¼ 1) or Pn;~amþ1 (for i ¼ 0), with the coefficients ~ amþ1 chosen in such a way that the m-qORF has zeros in X m . Clearly, the existence of the ðn þ i  1 : n  i  m; An ; FÞ-PRIQ (14) depends on the sequence of poles Anþi1 and the sequence of nodes X m (i.e., on whether the zeros of the m-qORF are all distinct real, and the corresponding weights are positive, for the given sequences Anþi1 and X m ). A special kind of PRIQs with nodes fixed in advance are the rational Gauss–Radau ðm ¼ 1Þ and Gauss–Lobatto ðm ¼ 2Þ quadrature rules. The existence of rational Gauss–Radau quadrature rules has been studied in [8, Sect. 5] for the case of F

i ¼ 0. In this reference it has been proved that, if P n;~a2 2 de Ln , there exist HPD linear functionals G, for which fuj gn1 j¼0 forms c

an orthonormal system in Ln1 with respect to the inner product hf ; giG ¼ Gffg g, and poles a 2 C0 , so that ZZan Pn;~a2 ?G Ln1 . For this reason, whenever a ðn  1 : n  1; An ; FÞ-PRIQ exists, it is a ðn : n  1; fa1 ; . . . ; an1 ; ag; GÞ-PRIQ (i.e., a rational Gaussian quadrature rule) too. Consequently, its nodes can be computed by solving a modified GEP of the form (6), with the coef½u

½u

½P

½P

ficients an1 and bn1 that appear in the last row of the matrices Jn and Sn replaced with an1;xn;1 and bn1;xn;1 (depending on the fixed node; see [8, Sect. 3]), while the corresponding weights again can be found by means of the corresponding eigenvectors. The aim of this section now is to study the existence and construction of the rational Gauss–Radau quadrature rules for i ¼ 1, and the rational Gauss–Lobatto quadrature rules for i ¼ 0. Additionally, we also consider the case in which the linear functional F is defined by a positive bounded Borel measure supported on an interval of the form ½a; b, with 1 6 a < b 6 1. 4.1. Rational Gauss–Radau quadrature rules for i ¼ 1 Recall from Section 3 that a necessary (but not sufficient) condition for the existence of a PRIQ is that for every function f in the domain of validity it should hold that f c is in the domain of validity too. Since a ðn : n  2; An ; FÞ-PRIQ should at least

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

10307

be exact for every f 2 Rn;n2 , it follows that it should also be exact for every f 2 Rn;n2 , and hence, for every f 2 Rn2;n2  Sfan1 ; an g. This leaves us with the following possible situations: (1) If fan1 ; an g  CR , then (a) for an ¼ an1 it follows that Rn2;n2  Sfan1 ; an g ¼ Rn1;n1 ð¼ Rn;n2 Þ, so that the ðn : n  2; An ; FÞ-PRIQ is a ðn  1 : n  1; An ; FÞ-PRIQ (i.e., a rational Gauss–Radau quadrature rule with i ¼ 0) too; (b) for an – an1 it follows that Rn2;n2  Sfan1 ; an g ¼ Rn;n , so that the ðn : n  2; An ; FÞ-PRIQ does not exist. (2) If an1 2 CR and an 2 R0 , it follows that Rn2;n2  Sfan1 ; an g ¼ Rn;n1 ð Rn;n2 Þ, so that the ðn : n  2; An ; FÞ-PRIQ is a ðn : n  1; An ; FÞ-PRIQ (i.e., the rational Gaussian quadrature rule) too. (3) If an1 2 R0 and an 2 CR , it follows that Rn2;n2  Sfan1 ; an g ¼ Rn;n2  Lfan g ð Rn;n2 Þ, so that the ðn : n  2; An ; FÞPRIQ is a ðn : n  1; fa1 ; . . . ; an2 ; an ; an1 g; FÞ-PRIQ (i.e., the rational Gaussian quadrature rule for the case in which the last two poles are interchanged) too. (4) If fan1 ; an g  R0 , then (a) for an ¼ an1 it follows that Rn2;n2  Sfan1 ; an g ¼ Rn1;n1 ð¼ Rn;n2 Þ, so that the ðn : n  2; An ; FÞ-PRIQ is a ðn  1 : n  1; An ; FÞ-PRIQ (i.e., a rational Gauss–Radau quadrature rule with i ¼ 0) too; (b) for an – an1 it follows that Rn2;n2  Sfan1 ; an g ¼ Rn;n2 . For the first two situations (1) and (2), the existence of the rational Gauss–Radau quadrature rule (with i ¼ 1) can easily be verified by checking the existence of the corresponding rational Gauss–Radau quadrature rule (with i ¼ 0) or rational F Gaussian quadrature rule for the given sequence of poles An and fixed node xn;1 . Thus, in what follows we can restrict ourselves to the case in which an1 2 R0 and an 2 Cf0;an1 g .4 We then have the following lemma. Lemma 8. A 1-qORF of the first kind

Q n;~a2 ðxÞ ¼ a0 un ðxÞ þ a1 un1 ðxÞ 2 de Ln satisfies a relation of the form

(" Q n;~a2 ðxÞ ¼

 E½Q n Z n ðxÞ

# ) Dn½Q  C n½Q  K n un1 ðxÞ  c 1þ u ðxÞ ; Z n1 ðxÞ Z n2 ðxÞ n2

where K n 2 C0 ,

E½Q n ¼  C ½Q n ¼

a0 En þ a1 =Z n ðan1 Þ 2 C0 ; Kn jK n j2 a0 En ½Q En1 ½a0 En

þ a1 =Z n ðan1 Þ

D½Q n ¼ 2 C0 ;

a0 En Dn þ a1 ; a0 En þ a1 =Z n ðan1 Þ ½Q 

and En1 ¼ K n En1 2 C0 :

Proof. From Theorem 1 it follows that

( Q n;~a2 ðxÞ ¼ a0 En Z n ðxÞ 1 þ

)  Dn 1 u ðxÞ  un2 ðxÞ þ a1 un1 ðxÞ; Z n1 ðxÞ n1 En1 Z cn2 ðxÞ

where a0 – 0 due to Q n;~a2 2 Ln n Ln1 . Setting

un1 ðxÞ ¼

  Z n ðxÞ 1 1 un1 ðxÞ ¼ Z n u ðxÞ; þ Z n ðxÞ Z n ðan1 Þ Z n1 ðxÞ n1

we obtain that

(" Q n;~a2 ðxÞ ¼ Z n ðxÞ

e En þ

# ) e en D C un1 ðxÞ  c n un2 ðxÞ ; Z n1 ðxÞ Z n2 ðxÞ

where

e E n ¼ a0 En þ

a1 ; Z n ðan1 Þ

e n ¼ a0 En Dn þ a1 D

and

e n ¼ a0 En : C En1

Further, for Q n;~a2 ¼ pqnn 2 de Ln we find with uk ¼ ppk that k

qn ðan1 Þ ¼ e E n an1 pn1 ðan1 Þ – 0;

4 Because of the re-ordering of the last two poles in the rational Gaussian quadrature rule for the case in which an1 2 R0 and an 2 CR , we will allow an to be arbitrary complex.

10308

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

 ½Q  ½Q  ½Q  ½Q  e e so that e E n – 0. The relation now follows by setting e E n ¼ E½Q h n K n (for an arbitrary K n 2 C0 ), D n ¼ Dn En K n , and C n ¼ En C n .

As a consequence of the previous lemma and Theorem 2, we now can prove the following. Lemma 9. Consider a 1-qORF of the first kind

Q n;~a2 ðxÞ ¼ a0 un ðxÞ þ a1 un1 ðxÞ ¼

qn ðxÞ

pn ðxÞ

2 de Ln :

Let a 2 C0 be chosen in such a way that qn ðaÞ – 0 and

 Tfag ¼ k  I mfa1 En =a0 g þ Tfan g n2 for some k 2 Rþ 0 . Then there exist HPD linear functionals G, for which fuj gj¼0 forms an orthonormal system in Ln2 and

un1 ?G Ln2 with respect to the inner product defined by hf ; giG ¼ Gffg c g, so that

v

½Q  n

:¼

a1 a0 En



1

1

an1



Za Zn

Q n;~a2 2 de L½na and

Za Zn

Q n;~a2 ?G Ln1 iff

 an 2 ð1; 1Þ.

Proof. Since an1 2 R0 , it follows from assumptions (A4) and (A5) in Theorem 2 that for every K n 2 C0 , there exist HPD linear c ½a e so that K u e de functionals G Ln1 with respect to the inner product defined by hf ; gi e ¼ Gffg g. So, let un be n n1 is a nORF in G given by

un½a ¼

E½na Z a ðxÞ Q n;~a2 ðxÞ ¼ E½na Z a ðxÞ Z ðxÞ n E½Q n

(" 1þ

) #  D½Q C ½Q n K n un1 ðxÞ  c n un2 ðxÞ ; Z n1 ðxÞ Z n2 ðxÞ

½a

where E½na 2 C0 . Then it holds that un satisfies assumptions (A1)–(A4) in Theorem 2. From assumption (A6) it now follows ½Q  1 2  that C ½Q n should be equal to ðK n En1 Þ . From the previous lemma we deduce that this can only be for jK n j ¼ 1 þ vn , which ½Q  implies that vn should be in ð1; 1Þ. Finally, we have that

ImfDn½Q  g

¼

n o I mfDn g þ I m aa0 1En ½Q

1 þ vn

¼

I mfa1 En =a0 g þ Tfan g ð1 þ vn ÞjEn j2 ½Q 



Tfan2 g jEn1 j2 ½Q 

¼

Tfag kð1 þ vn ÞjEn j2 ½Q 



Tfan2 g jEn1 j2 ½Q

;

where the second equality follows from assumption (A5) applied on Dn . Consequently, assumption (A5) in Theorem 2 is sat½Q  isfied too for every E½na 2 C0 if Tfag ¼ 0, respectively for jE½na j2 ¼ kð1 þ vn ÞjEn j2 2 Rþ h 0 if Tfag – 0. The previous lemma now allows us to prove our first main result. Theorem 10. Consider a quadrature formula of the form F

F

Fn ff g ¼ kn;1 f ðxn;1 Þ þ

n X F F kn;k f ðxn;k Þ;

ð15Þ

k¼2

with fixed node in xn;1 2 RAn . Further, suppose that the nodes fxn;k gnk¼1 are zeros of the 1-qORF of the first kind: F

F

Q n;~a2 ðxÞ ¼ a0 un ðxÞ þ a1 un1 ðxÞ;

½Q

An;1

    a1 un ðxÞ F ½Q :¼ ¼ ¼ An;1 ðxn;1 Þ a0 un1 ðxÞ x¼xF n;1

and let the corresponding weights be given by (9). Then (15) is a rational Gauss–Radau quadrature rule with Rn ff g ¼ Fff g  Fn ff g ¼ 0 for (at least) every f 2 Rn;n2 iff the following conditions are satisfied: ½Q 

½Q 

1. An;1 2 C and An;1 –  En Z n ðan1 Þ; 2.

½Q  n

 v 2 ð1; 1Þ, where v½Q n is defined as before in Lemma 9; ½Q 

3. I mfAn;1 En g ¼ Tfan g; 4. ½Q 



An;1 – 

un ðxÞ un1 ðxÞ



 ~ 2 ðAn2 [ f1gÞ \ Rfan1 ;an g : for every a

~ x¼a

If an 2 R0 , then the first and third condition are automatically satisfied when the second condition is satisfied. On the other hand, if an R R, then the first three conditions are satisfied iff ½Q 

An;1 ¼



1 1 1 2 C0 :  En an1 an

10309

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

Proof. First, note that Q n;~a2 2 de Ln iff the first condition is satisfied. Next, from Lemma 9 it follows that there exist poles a 2 C0 , constants fkn1 ; kn g  C0 , and HPD linear functionals G so that the sequence of rational functions fu0 ; . . . ; un2 ; kn1 un1 ; kn ZZan Q n;~a2 g, with ZZan Q n;~a2 2 de L½na , forms an orthonormal system in L½na with respect to the inner prodc

uct defined by hf ; giG ¼ Gffg g iff the second condition is satisfied. The third condition then follows from the expression for Tfag in Lemma 9 and the fact that the nodes of an (n)ORF are all real and distinct iff the last pole is in R0 (i.e, iff Tfag ¼ 0). Finally, we have that Q n;~a2 2 r Ln iff the last condition is satisfied too. Thus, the conditions given in the statement ensure the existence of the rational Gaussian quadrature formula

Gn ff g :¼

n X

G kG n;k f ðxn;k Þ;

k¼1

where xG ¼ xn;k for k ¼ 1; . . . ; n, and ðGff g  Gn ff gÞ ¼ 0 for every f 2 Ln½a  Lcn1 . Furthermore, because n;k F

Gff g ¼ Fff g for every f 2 Rn1;n2 ; it follows that

Fff g ¼ Gn ff g for every f 2 Rn1;n2 Ln1 :

ð16Þ

Due to the uniqueness of a ðn  1; 0; An ; FÞ-PRIQ with n fixed nodes, it follows that F

þ kn;k ¼ kG n;k 2 R0 ;

k ¼ 1; . . . ; n:

Moreover, from the proof of Theorem 6 we deduce that the equality in (16) holds for every f 2 Rn;n2 due to the fact that the nodes are zeros of a 1-qORF of the first kind. Suppose now that an 2 R0 . From the second condition it then follows that

I mfv½Q n g ¼ 0 $ I mfa1 En =a0 g

1



1

jEn j2 an1



1

an



¼ 0:

Since we assumed that an – an1 , this implies that

I mfa1 En =a0 g ¼ 0 ¼ Tfan g: On the other hand, for an R R, it follows from the second and third condition that

(

I mfa1 En =a0 gð1=an1  Ref1=an gÞ  Refa1 En =a0 gI mf1=an g ¼ 0 I mfa1 En =a0 g ¼ I mf1=an g – 0 (  I mfa1 En =a0 g 1=an1  Ref1=an g  Refa1 En =a0 g ¼ 0 () I mfa1 En =a0 g ¼ I mf1=an g – 0 ( a 1 En 1 1 Refa1 En =a0 g ¼ 1=an1  Ref1=an g () ¼  : () a0 an1 an I mfa1 En =a0 g ¼ I mf1=an g – 0

As a result,

A

½Q 

 2 n;1 v½Q 2 Rþ n ¼ En Z n ðan1 Þ ¼ jEn Z n ðan1 Þj 0.

Finally, suppose (15) is a rational Gauss–Radau quadrature with Rn ff g ¼ 0 for (at least) every f 2 Rn;n2 , while the second condition is not satisfied. We then have that

kn;1 jun1 ðxn;1 Þj2 þ F

F

n X F F kn;k jun1 ðxn;k Þj2 ¼ jkn1 j2 > 0: k¼2 ½a



Consider now a rational function un n1 2 de Ln , defined by

("

½an1  ðxÞ n

u

¼

E½nan1  Z n1 ðxÞ

# ) Dn½an1  1 kn1 un1 ðxÞ  1þ un2 ðxÞ ; Z n1 ðxÞ kn1 En1 Z ck2 ðxÞ

e so where E½nan1  2 C0 and I mfD½nan1  g ¼  jkTfaEn2 gj2 . From Theorem 2 it then follows that there exist HPD linear functionals G n1 n1

½a

that the sequence of rational functions fu0 ; . . . ; un2 ; kn1 un1 ; un n1 g forms an orthonormal system in L½nan1  with respect c e to the inner product defined by hf ; gi e ¼ Gffg g. Moreover, the quadrature rule (15) is then a ðn  1; n  1; fa1 ; . . . ; G e a ; a g; GÞ-PRIQ. Since assumption (8) is satisfied, it follows from Theorem 7 that the nodes fxF gn are zeros of a 1n1

n1

n;k k¼1

qORF of the second kind of the form

Pn;~b2 ðxÞ ¼ b0 un½an1  ðxÞ þ b1



Z n1 ðxÞ u ðxÞ 2 r Ln½an1  ; Z cn1 ðxÞ n1

10310

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

n ðxÞ ~ 2 R0 so that for every kn 2 C0 , the with ZZn1 P b2 ðxÞ  Q n;~a2 ðxÞ. In [8, Thm. 7] it has been proved then that there exists poles a ðxÞ n;~

a^ P n;~b2 g forms an orthonormal system in Ln½a with respect to inner sequence of rational functions fu0 ; . . . ; un2 ; kn1 un1 ; kn ZZn1 c ^ products defined by certain HPD linear functionals G on L  L . But this is in contradiction with Lemma 9 and our assumption

^

that the second condition is not satisfied. This concludes the proof.

h

4.2. Rational Gauss–Lobatto quadrature rules for i ¼ 0 In a similar way as above, we can now study the existence of rational Gauss–Lobatto quadrature rules for i ¼ 0. Since a ðn  1 : n  2; An ; FÞ-PRIQ should at least be exact for every f 2 Rn1;n2 , it follows that it should also be exact for every f 2 Rn2;n1 , and hence, for every f 2 Rn2;n2  Sfan1 g. This leaves us with the following possible situations: 1. If an1 R R, it follows that Rn2;n2  Sfan1 g ¼ Rn1;n1 ð Rn;n2 Þ, so that the ðn  1 : n  2; An ; FÞ-PRIQ is a ðn  1 : n  1; An ; FÞ-PRIQ (i.e., a rational Gauss–Radau quadrature rule with i ¼ 0) too. 2. If an1 2 R0 , it follows that Rn2;n2  Sfan1 g ¼ Rn1;n2 . For the first situation, the existence of the rational Gauss–Lobatto quadrature rule (with i ¼ 0) can easily be verified by checking the existence of the corresponding rational Gauss–Radau quadrature rule (with i ¼ 0) for the given sequence of F F poles An and fixed nodes X 2 ¼ fxn;1 ; xn;2 g. Thus, in what follows we can restrict ourselves to the case in which an1 2 R0 . We then have the following lemma. Lemma 11. A 2-qORF of the second kind

Pn;~a3 ðxÞ ¼ a0 un ðxÞ þ a1

Z n ðxÞ Z ðxÞ u ðxÞ þ a2 c n u ðxÞ 2 de Ln ; Z cn1 ðxÞ n1 Z n2 ðxÞ n2

with an1 2 R0 (hence, Z cn1 ðxÞ  Z n1 ðxÞ), satisfies a relation of the form

("

Pn;~a3 ðxÞ ¼ E½P n Z n ðxÞ

1þ

# ) D½P C ½P n K n un1 ðxÞ  c n un2 ðxÞ ; Z n1 ðxÞ Z n2 ðxÞ

ð17Þ

where K n 2 C0 ,

E½P n ¼ C ½P n ¼

a0 En 2 C0 ; Kn jK n j2 ½P

En1

a1 D½P ; n ¼ Dn þ a0 En ! a2 En1 ½P 2 C; and En1 ¼ K n En1 2 C0 : 1 a0 En

Proof. (The proof is similar to the proof of Lemma 8.) From Theorem 1 it follows that

(

)  Dn 1 Z n ðxÞ Z ðxÞ un1 ðxÞ  u ðxÞ þ a1 u ðxÞ þ a2 c n u ðxÞ n2 Z n1 ðxÞ n1 Z n1 ðxÞ Z n2 ðxÞ n2 En1 Z cn2 ðxÞ # (" ) en en C D e u ðxÞ  c u ðxÞ ; ¼ Z n ðxÞ E n þ Z n1 ðxÞ n1 Z n2 ðxÞ n2

Pn;~a3 ðxÞ ¼ a0 En Z n ðxÞ

1þ

where

e E n ¼ a0 En ;

e n ¼ a0 En Dn þ a1 D

and

e n ¼ a0 En  a2 : C En1

Further, for Pn;~a3 ¼ pqnn 2 de Ln we find with uk ¼ ppk that k

qn ðan1 Þ ¼ e E n an1 pn1 ðan1 Þ – 0; ½P ½P ½P ½P e e so that e E n – 0. The relation now follows by setting e E n ¼ E½P h n K n (for an arbitrary K n 2 C0 ), D n ¼ Dn En K n , and C n ¼ En C n .

As a consequence of the previous lemma and Theorem 2, we now can prove the following. Lemma 12. Consider a 2-qORF of the second kind

Pn;~a3 ðxÞ ¼ a0 un ðxÞ þ a1

Z n ðxÞ Z ðxÞ q ðxÞ u ðxÞ þ a2 c n u ðxÞ ¼ n 2 de Ln : pn ðxÞ Z cn1 ðxÞ n1 Z n2 ðxÞ n2

10311

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

Let a 2 C0 be chosen in such a way that qn ðaÞ – 0 and

( ) ( ) ! a1 En a2 En þ Tfan g  Re Tfan2 g a0 a0 En1

Tfag ¼ k  I m

n2 for some k 2 Rþ 0 . Then there exist HPD linear functionals G, for which fuj gj¼0 forms an orthonormal system in Ln2 and

un1 ?G Ln2 with respect to the inner product defined by hf ; giG ¼ Gffg c g, so that a2 En1 v½P n :¼ a0 En 2 ð1; 1Þ.

Za Zn

Pn;~a3 2 de L½na and

Za Zn

P n;~a3 ?G Ln1 iff

Proof. (The proof is similar to the proof of Lemma 9.) Since an1 2 R0 , it follows from assumptions (A4) and (A5) in Theorem e so that K u 2 that for every K 2 C , there exist HPD linear functionals G is a nORF in de L with respect to the inner n

0

n

n1

n1

c ½a e g. So, let un be given by product defined by hf ; gi e ¼ Gffg G

u½na ¼

En½a E½P n

Z a ðxÞ P ~ ðxÞ ¼ En½a Z a ðxÞ Z n ðxÞ n;a3

("

1þ

) # D½P C ½P n K n un1 ðxÞ  c n un2 ðxÞ ; Z n1 ðxÞ Z n2 ðxÞ

½a

where E½na 2 C0 . Then it holds that un satisfies assumptions (A1)–(A4) in Theorem 2. From assumption (A6) it now follows that

C ½P n

½P

should be equal to ðK n En1 Þ . From the previous lemma we deduce that this can only be for jK n j2 ¼ 1  vn , which

implies that

1

½P n

v should be in ð1; 1Þ. Finally, we have that

n o a1 En     I m þ Tfan g a0 a a Tf a g Tf a g Tfan2 g Tfan2 g 1 1 n n2 ¼ I m þ I mfD½P g ¼ ImfD g þ I m  ¼  v½P  ½P 2 n n n 2 2 2 a0 En a0 En jEn j jEn1 j jEn1 j2 jEn j jEn1 j n o a1 En En 2 ½P Im a0 þ Tfan g  j En1 j vn Tfan2 g Tfa g Tfag Tfa g n2 n2 ¼  ½P 2 ¼  ½P 2 ; kjEn j2 jEn j2 jEn1 j jEn1 j where the second equality follows from assumption (A5) applied on Dn , and the last equality is due to

   En 2 ½P a2 En   2 R: E  vn ¼ a E n1 0 n1 Consequently, assumption (A5) in Theorem 2 is satisfied too for every E½na 2 C0 if Tfag ¼ 0, respectively for jE½na j2 ¼ kjEn j2 2 Rþ h 0 if Tfag – 0. Before we can prove our second main result, we need the following lemma. Lemma 13. Suppose the nORFs uk 2 de Lk ; k ¼ 1; . . . ; n. Then the nORFs un and un1 have no common zeros. ðkÞ

ðkÞ

Proof. Set uk ¼ ppk . From Theorem 1 it then follows that there exist polynomials q1 2 P 1 and constants q0 2 C0 so that k

ðkÞ

ðkÞ

pk ðxÞ ¼ q1 ðxÞpk1 ðxÞ þ q0 ð1  x=ak1 Þð1  x=ak2 Þpk2 ðxÞ;

k ¼ 2; . . . ; n:

So, let k ¼ n and suppose that pk1 ðaÞ ¼ pk ðaÞ ¼ 0 for a 2 C (recall that ‘pk ð1Þ ¼ 0’ means that ‘pk 2 P k1 ’). Since uk 2 e Lk and uk1 2 d Lk1 , it follows that a – ak1 and a – ak2 , and hence, that pk2 ðaÞ ¼ pk1 ðaÞ ¼ 0. Thus, continuing for k ¼ n  1; . . . ; 2, we find that a R fak1 ; ak2 gnk¼2 , and that p0 ðaÞ ¼ p1 ðaÞ ¼ 0. This, however, contradicts the fact that u0  p0 2 C0 . h Finally, we are able to prove the following. Theorem 14. Consider a quadrature formula of the form F

F

F

F

Fn ff g ¼ kn;1 f ðxn;1 Þ þ kn;2 f ðxn;2 Þ þ

n X F F kn;k f ðxn;k Þ;

ð18Þ

k¼3 F F with fixed distinct nodes in fxn;1 ; xn;2 g 2 RAn . Further, let the functions að~ xÞ and bð~ xÞ be defined by

að~xÞ ¼



Z n1 ðxÞ Z n ðxÞ





 x¼~x

un ðxÞ un1 ðxÞ



x¼~x

and bð~xÞ ¼



Z n1 ðxÞ Z cn2 ðxÞ





 x¼~x

un2 ðxÞ un1 ðxÞ



: x¼~x

Suppose the nodes fxn;k gnk¼1 are zeros of the 2-qORF of the second kind: F

Pn;~a3 ðxÞ ¼ a0 un ðxÞ þ a1

 Z n ðxÞ Z ðxÞ u ðxÞ þ a2 c n un2 ðxÞ Z n1 ðxÞ  Z cn1 ðxÞ ; Z n1 ðxÞ n1 Z n2 ðxÞ

10312

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

with ½P

An;1 :¼

F F F F F F   a1 aðxn;1 Þbðxn;2 Þ  aðxn;2 Þbðxn;1 Þ  a2 aðxn;2 Þ  aðxn;1 Þ  F F F F ½P ½P ½P ¼ and An;2 :¼ ¼ ¼ An;1 ðxn;1 ; xn;2 Þ ¼ An;2 ðxn;1 ; xn;2 Þ F F F F a0 a0 bðxn;1 Þ  bðxn;2 Þ bðxn;1 Þ  bðxn;2 Þ

and let the corresponding weights be given by (11). Then (18) is a rational Gauss–Lobatto quadrature rule with Rn ff g ¼ Fff g  Fn ff g ¼ 0 for (at least) every f 2 Rn1;n2 iff the following conditions are satisfied: 1. 2.

n o F F F F bðxn;1 Þ; bðxn;2 Þ  C, and bðxn;1 Þ  bðxn;2 Þ – 0; ½P

En1 ½Z n1 ðan Þ þ Dn þ An;1 =En  ½P n

1v

– bðan Þ if an – an1 and an 2 R0 ;

v½P n is defined as before in Lemma 12; ½P 3. vn 2 ð1; 1Þ; ½P ½P n Tfan2 g ¼ Tfan g; 4. I mfAn;1 En g  An;2 EEn1

where

 ½P ~ Þ þ A½P ~ ~ 5. aða n;1 þ An;2 bðaÞ – 0 for every a 2 ðAn2 [ f1gÞ \ Rfan1 ;an g . ½P

½P

Proof. First, note that Pn;~a3 R de Ln if a0 ¼ 0. Hence, we may assume that a0 – 0 so that the coefficients An;1 and An;2 are well F defined. We then have that xn;j ; j ¼ 1; 2, are zeros of P n;~a3 iff F

un ðxFn;j Þ þ A½P n;1

Z n ðxn;j Þ F

Z n1 ðxn;j Þ

F

un1 ðxFn;j Þ þ A½P n;2

Z n ðxn;j Þ Z cn2 ðxn;j Þ F

un2 ðxFn;j Þ ¼ 0; j ¼ 1; 2:

ð19Þ

Suppose now that the jth node is a zero of un1 . From (17) it then follows that

  F   Z n xn;j   un2 xFn;j : 0 ¼ C ½P n F c Z n2 xn;j

Due to the previous lemma, a zero of un1 cannot be a zero of un2 , so that C ½P n ¼ 0. As a result, the third condition is not satisfied. For this reason, the first part of the first condition ensures that none of the nodes are zeros of un1 , and that   F   Z n xn;j ½P ½P ½P vn – 1. The expressions for An;1 and An;2 then follow by dividing the left-hand side in (19) with   un1 xFn;j , and solv½P

ing for An;1

Z n1

    F F ½P and An;2 . Further, for b xn;1 ¼ b xn;2 it holds that

F

xn;j

    h    i F F F F a xn;1  a xn;2 ¼ En Z n1 xn;1  Z n1 xn;2 – 0; so that a0 – 0 iff the second part of the first condition is satisfied. Next, for an – an1 we deduce from (17) that P n;~a3 2 Ln n Ln1 (and hence, in de Ln due to a0 – 0) iff the second condition is h i n2 ðxÞ satisfied. Note that bðan Þ is well defined. Indeed, since un1 2 de Ln1 , it follows that u is finite whenever an ¼ an2 . u ðxÞ n1

x¼an

Similar as in the proof of Theorem 10, we then have that 1. conditions 3–5 ensure the existence of a rational Gaussian quadrature formula Gn ff g so that Fff g ¼ Gff g ¼ Gn ff g for F F G þ every f 2 Rn1;n2 , and hence, xn;k ¼ xG n;k and kn;k ¼ kn;k 2 R0 for k ¼ 1; . . . ; n; 2. supposing the rational Gauss–Lobatto quadrature exists, while the fixed nodes are zeros of a 2-qORF of the second kind that does not satisfy the third condition, leads to a contradiction. Finally, note that we only need to verify the second condition if an 2 R0 , because conditions 3 and 4 ensure that all the zeros are real. This ends the proof. h 4.3. Measures supported on an interval So far we considered an arbitrary HPD linear functional F. In many practical situations, however, this linear functional F is defined by

Fff g ¼

Z a

b

f ðxÞdlðxÞ;

10313

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

where l is a positive bounded Borel measure (in short, a measure) with suppðlÞ # ½a; bˆR, and the poles A ¼ fa1 ; a2 ; . . .g are bounded away from ½a; b. In this case, it is not only required that the nodes in an nth PRIQ are all real, but in addition that they are also all inside the interval ½a; b. So, the aim of this subsection is to derive additional conditions to ensure that the zeros of an m-qORF are all in ½a; b. Since for every a and b, the interval ½a; b can be mapped onto the interval I ¼ ½1; 1 by means of a properly chosen polynomial or rational transformation, we can-without lost of generality-restrict ourselves in the remainder to the interval I. It is well known that the zeros of an nth ORF with respect to a measure l on I are all inside ð1; 1Þ whenever an 2 RI . Further, in [23, Thm. 4.1] it is proved that a 1-qORF of the second kind Pn;~a2 2 r Ln with all real zeros has at least ðn  1Þ zeros inside I, and that there exist ða1 =a0 Þ 2 C such that all its zeros are inside I. A similar result holds for m-qORFs in general. (The proof of the following lemma is practically the same as in [23, Thm. 4.1], and hence, we omit it.) Lemma 15. Suppose that the zeros of an m-qORF of the first or second kind in r Ln are all real. Further, assume that the poles fanðmþ1iÞ ; . . . ; ani g, with i ¼ 1 for the first kind or i ¼ 0 for the second kind, are all real and/or appear in complex conjugate pairs, and outside I. Then at least ðn  mÞ of the zeros are inside the interval I. Note that the condition on the poles in the previous lemma is needed to fulfill the necessary condition for the existence of a PRIQ with domain of validity Rni;nðmþ1iÞ . Whenever the condition is not satisfied, we do not need to investigate the location of the zeros of the m-qORF due to the fact that the corresponding PRIQ cannot exist anyway. First, consider the case of rational Gauss–Radau for i ¼ 1 (the case in which i ¼ 0 has already been dealt with in [8]). The condition on the poles in the previous lemma then says that either an ¼ an1 , or fan1 ; an g 2 RI . Since the first situation (i.e., an ¼ an1 ) can be considered as a special case of rational Gauss–Radau for i ¼ 0, we can restrict ourselves to the second sit½Q  uation, with an – an1 . From Lemma 9 it then follows that the zeros are all real iff vn 2 ð1; 1Þ, while they are for sure all in I if a1 ¼ 0. Thus, we now have the following. Theorem 16. Suppose that fan1 ; an g 2 RI , with an – an1 . Then a 1-qORF Q n;~a2 ðxÞ ¼ a0 un ðxÞ þ a1 un1 ðxÞ 2 r Ln has all its zeros inside I iff

v½Qn  2 ðm; MÞ; with

 v½Q n as defined before in Lemma 9, and

m ¼ minfcð1Þ; cð1Þg 2 ð1; 0Þ and M ¼ maxfcð1Þ; cð1Þg 2 ð0; 1Þ; where

  un ðxÞ 1 1 cð~xÞ ¼ E1 : n ðan  an1 Þ un1 ðxÞ x¼~x  v½Q n tending to 1, it follows from the proof of Lemma 8 that one of the zeros of Q n;~ a2 ðxÞ tends to an1 R I. On the ½Q  other hand, for vn tending to 1 (i.e., for a0 tending to 0), it is clear that one of the zeros tends to an R I. Since all the zeros are ½Q  inside ð1; 1Þ for vn ¼ 0, and taking into account that the zeros of Q n;~a2 ðxÞ change continuously when continuously changing ½Q     ~ ½Q ^ ½Q vn , it follows that somewhere in between of the limiting values 1 and 1 for v½Q n there exist values v n and v n such that

Proof. For

Q n;~a~2 ð1Þ ¼ 0 and Q n;~a^2 ð1Þ ¼ 0 respectively. The statement now easily follows. h Finally, consider the case of rational Gauss–Lobatto for i ¼ 0. In this case, the condition on the poles in Lemma 15 says that

an1 2 RI . From Lemma 12 it now follows that the zeros are all real iff v½P n 2 ð1; 1Þ and



 En 2 v½P n j En1 j Tfan2 g  Tfan g , with Y n 2 R. First, we need the following lemma. ½P

½P

a1 a0

En ¼ Y n þ i

Lemma 17. For each vn 2 ð1; 1 there exists a unique constant Y n;n ðvn Þ 2 R, with n 2 f 1g, such that the corresponding 2-qORF P n;~a3 ðxÞ has a zero in x ¼ n, while the remaining ðn  1Þ zeros are all real. Proof. Clearly, a 2-qORF P n;~a3 ðxÞ has a zero in x ¼ n iff

  a1 a2 ¼  aðnÞ þ bðnÞ ; a0 a0 where að~ xÞ and bð~ xÞ are as defined before in Theorem 14. Multiplying both sides with En yields

" #    En 2  a1   bðnÞEn1 ¼: Z n;n v½P ; En ¼  aðnÞEn þ v½P n  n  a0 En1

10314

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

so that

RefZ n;n



"



#    En 2   RefbðnÞEn1 g ¼ Y n;n v½P : n En1 

½P  v½P n g ¼  RefaðnÞEn g þ vn 

Thus, it remains to prove that the remaining zeros are all real. For this we need to prove that

ImfZ n;n



   En 2  Tfan2 g  Tfan g; En1 



½P  v½P n g ¼ vn 

8v½P n 2 ð1; 1:

ð20Þ

The statement in the lemma has already been proved in [8] for the case of a 1-qORF P n;~a2 ðxÞ, which implies that the equality ½P ½P in (20) surely holds true for vn ¼ 0. On the other hand, for vn ¼ 1 it follows from the proof of Lemma 9 that Pn;~a3 ðxÞ has one zero in

x¼

Dn þ aa10 En jEn j2   ; 1  Dn þ aa10 En jEn j2 a1 n1

while the hremaining iðn  1Þ zeros are in ð1; 1Þ due to the fact that they are the zeros of un1 . As a result, ½P Z n;n ð1Þ ¼  1n=nan1 þ Dn jEn j2 , while I mfZ n;n ð1Þg ¼ I mfDn gjEn j2 . Thus, the equality in (20) holds true for vn ¼ 1 too, which ½P implies that the equality holds true for every vn 2 ð1; 1. h With this we now can prove the following. Z n ðxÞ n ðxÞ r Theorem 18. Suppose that an1 2 RI . Then a 2-qORF P n;~a3 ðxÞ ¼ a0 un ðxÞ þ a1 ZZn1 ðzÞ un1 ðxÞ þ a2 Z cn2 ðxÞ un2 ðxÞ 2 Ln has all its zeros inside I iff

a1 En ¼ Y n þ i a0

v½P and n 2 ½X m ; 1Þ where



v½P n j

En 2 ½P j Tfan2 g  Tfan g ; Y n 2 ½mðv½P n Þ; Mðvn Þ; En1

v½P n is defined as before in Lemma 12, aðnÞ  aðnÞ

Xm ¼

a2

< 0;

n 2 f 1g;

ð21Þ

aðnÞ  aðnÞ  2En n a2 n11 n1

with að~ xÞ as defined before in Theorem 14, and ½P ½P mðv½P and Mðv½P n Þ ¼ min Y n;n ðvn Þ n Þ ¼ max Y n;n ðvn Þ; n2f 1g

n2f 1g

with

Y n;n ðv½P n Þ ¼





2 v½P n  1 RefaðnÞEn g  jEn j



n þ RefDn g : 1  n=an1 ½P

Proof. Consider the half plane ð1; 1Þ  ð1; 1Þ, i.e., the collection of points ðvn ; Y n Þ for which the corresponding 2-qORF Pn;~a3 ðxÞ has n real distinct zeros. Due to the previous lemma, there exist functions

" #  2   ½P  En  Y n;n ðv½P Þ ¼  RefaðnÞE g þ v RefbðnÞE g ; n n1 n n  En1  ½P

n 2 f 1g;

½P

such that the points ðvn ; Y n;n ðvn ÞÞ correspond with 2-qORFs for which P n;~a3 ðnÞ ¼ 0. Furthermore, from [6, Thm. 3.4] it follows that there exists a unique X m 2 ð1; 1Þ for which Y n;1 ðX m Þ ¼ Y n;1 ðX m Þ, such that Pn;~a3 ð1Þ ¼ 0 ¼ Pn;~a3 ð1Þ, while the remaining n  2 zeros are in ð1; 1Þ. From Lemma 14 we deduce that

Xm ¼



aðnÞ  aðnÞ En1 : bðnÞ  bðnÞ En

Thus, the graphs of Y n;1 1. S1 :¼ 2. S2 :¼

n n









½P and Y n;1 vn divide the half plane into four sub-planes: v½P n





½P ½P ½P v½P n ; Y n : vn < X m and mðvn Þ 6 Y n 6 M vn ½P

ðvn ; Y n Þ : Y n < m



v½P n

o

o ;

; o ½P ½P ; 3. S3 :¼ ðvn ; Y n Þ : Y n > M vn n     o ½P ½P ½P 4. S4 :¼ 6 Y n 6 M vn . vn ; Y n : vn P X m and m v½P n n



K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

10315

Taking into account that the zeros of P n;~a3 ðxÞ change continuously when continuously changing the parameter a0 ; a1 and a2 , and using the fact that one of the zeros should tend to an1 R I whenever a0 tends to 0, it follows that P n;~a3 ðxÞ has at least one ½P ½P zero outside I whenever ðvn ; Y n Þ 2 S1 [ S2 [ S3 , and has all its zeros in I for ðvn ; Y n Þ 2 S4 . Finally, from the proof of the previous lemma it follows that

     En 2 n  bðnÞEn1 ¼ Z n;n ð1Þ ¼ aðnÞEn þ  þ Dn jEn j2 ;  1  n=an1 En1 so that

     En 2 n 2   En1  bðnÞEn1 ¼ 1  n=an1 þ Dn jEn j  aðnÞEn :



4.4. Construction of rational Gauss-type quadrature rules In the first two subsections, we related the rational Gauss–Radau and Gauss–Lobatto quadrature rules (for respectively i ¼ 1 and i ¼ 0) with certain rational Gaussian quadrature rules by means of the Lemmas 9 and 12 respectively. Supposing the quadrature rule exists for a given (set of) node(s), the aim of this subsection is then to investigate the construction of the quadrature rule (i.e., the computation of the nodes and weights in the quadrature formula). For this, let X be fixed to either Q ½Q  ½P ½Q  ½P or P, and suppose that the matrices Jn ; Dn ; In and Sn , are defined as above in Theorem 3. Further, let vn ; vn ; An;1 and An;1 , be defined as before in Lemmas 9 and 12, and Theorems 10 and 14, and set ½u

½X

cn2 ¼

cn2

; ½X

gK n

½X

½X

an1 ¼ cn2 ;

  ½X ½u  ½X ½u and bn1 ¼ bn1  An;1 cn1 =L½X n ;

ð22Þ

where jgj ¼ 1,

K ½X n

8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi >  < 1 þ v½Q X ¼ Q; n ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > ½P : 1= 1  v X ¼ P; n

( L½X n

and

¼

½Q 

1 þ vn

X ¼ Q;

1

X ¼ P:

With this we can define the modified matrices

0 J½X n

B ¼B @

~ 0n2

Jn1

1

C ½X ~ cn2 C A; ½X bn1

½X ~ 0Tn2 an1

0 @ S½X n ¼

Sn1 ~ 0Tn2 Zc

~ 0n2

½X

an1

n2 ðan2 Þ

0

1 A;

and ½X ½X B½X n ¼ Jn Dn þ In  Sn ;

where ~ 0n2 ¼ ð0 . . . 0ÞT 2 Cn2

and modified column vector



T

½X ~ ½X u : n ðxÞ ¼ u0 ðxÞu1 ðxÞ . . . K n un1 ðxÞ

We then have the following theorem for the construction of the rational Gauss–Radau and Gauss–Lobatto quadrature rules from the previous subsections. F

Theorem 19. Given the notation introduced above, the nodes fxn;k g in the rational Gauss–Radau quadrature rule (15) (X ¼ Q ), respectively in the rational Gauss–Lobatto quadrature rule (18) ðX ¼ PÞ, are eigenvalues of the GEP ½X

½X

F

½X ~ ~ J½X n v n;k ¼ xn;k Bn v n;k ;

with ½X ~ v ½X n;k ¼ gn

( ) n2   2 h i2   2 1=2     X     F F F ½X  ¼ 1; ~ ½X  un1 xn;k  u xn;k ; g½X uj xn;k  þ K n n n j¼0

the corresponding normalized eigenvector. The corresponding weights satisfy

   ½X 2 F kn;k ¼ v n;k  k1k2F ; where

v ½X n;k

denotes the first entry in ~ v ½X . n;k

ð23Þ

10316

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

Proof. The statement directly follows from Lemmas 8 and 9 for X ¼ Q (respectively from Lemmas 11 and 12 for X ¼ P), together with (13) and Theorem 1. h

5. Numerical examples In the numerical examples that follow, we construct rational Gauss-type quadrature formulas for the approximation of

Fff g ¼ J w ðf Þ :¼

Z

f ðxÞwðxÞdx;

I

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where I ¼ ½1; 1 and wðxÞ is the Chebyshev weight function of the first kind wðxÞ ¼ 1= 1  x2 . Let x ¼ JðzÞ denote the Joukowski Transformation x ¼ 12 ðz þ z1 Þ, which maps the open unit disc D ¼ fz 2 C : jzj < 1g onto the cut Riemann sphere CI and the unit circle T ¼ fz 2 C : jzj ¼ 1g onto the interval I. The inverse mapping is denoted by z ¼ J inv ðxÞ and is chosen so that z 2 D if x 2 CI . Further, with the sequence A ¼ fa1 ; a2 ; . . .g  CI we associate a sequence B ¼ fb1 ; b2 ; . . .g  D, so that bk ¼ J inv ðak Þ. The so-called Chebyshev nORFs (with respect to the Chebyshev weight function c wðxÞ and inner product hf ; giw ¼ J w ðfg Þ), with arbitrary complex poles outside I, are then given by (see [10])

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

1  jbk j2 zBck1 ðzÞ 1 uk ðxÞ ¼ þ ; 1  bk z ðz  bk ÞBk1 ðzÞ 2p

1 u0 ðxÞ  pffiffiffiffi ;

p

k P 1;

ð24Þ

where

B0 ðzÞ  1;

Bk ðzÞ ¼

k Y z  bj

1  bj z

j¼1

k P 1:

;

These nORFs are regular and satisfy the three-term recurrence relation (3) with coefficients (see [21, Section 4] and [24, Thm. 3.5])

Ek ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð1  jbk j2 Þð1  jbk1 j2 Þð1  bk bk1 Þ ð1 þ b2k Þð1 þ b2k1 Þ

;

k>1

and

Dk ¼ 

1 þ b2k1 2

2ð1  jbk1 j Þ



ð1  jbk1 j2 Þðbk þ bk2 Þ þ 2Refbk1 gð1  bk bk2 Þ ; ð1  bk bk1 Þð1  bk1 bk2 Þ

k>1

were we put by convention b0 ¼ J inv ða0 Þ ¼ 0. For k ¼ 1, the coefficients are given by

E1 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð1  jb1 j2 Þ 1 þ b21

and D1 ¼ b1 :

Unless stated differently, the computations in the examples were done in double precision using MATLABÒ7.

5

Example 1. For the first example, we consider the sequence of poles A1 ¼ fa; a; a; a; . . .g  RI , with a > 1. We then numerically approximate

J w ðuk ucl Þ ¼ dk;l ;

n  2 6 k 6 n;

n  2 6 l < n;

where the uk ’s are the Chebyshev nORFs (24) with poles in A1 , by means of an nth rational Gauss–Radau quadrature formula

J n ðuk ucl Þ ¼ kn;1 uk ðxn;1 Þucl ðxn;1 Þ þ

n X kn;j uk ðxn;j Þucl ðxn;j Þ; j¼2

with fixed node in xn;1 ¼ 1 or xn;1 ¼ 1. Thus, for n > 1 we have with b ¼ J inv ðaÞ 2 ð0; 1Þ that

( 1þb ½Q An;1 ð1Þ

¼

½Q 1=An;1 ð1Þ

¼

1b

n ¼ odd;

1b 1þb

n ¼ even:

Further, the coefficients (4) are given by

1 þ b2 ½u  ½u c0 ¼ a1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2ð1  b2 Þ 5

½u

½u

ck1 ¼ ak ¼

1 þ b2 2ð1  b2 Þ

MATLAB is a registered trademark of The MathWorks, Inc.

;

k>1

10317

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

and ½u 

½u

b0 ¼ b;

b1 ¼ 

bð1 þ b2 Þ 2ð1  b2 Þ

;

½u

bk ¼ 0;

k > 1:

½Q 

Thus, we have that An;1 ð 1Þ 2 C, while

v

½Q n ð1Þ

¼

8 < :

2b ð1þbÞ2

2b ð1bÞ2

n ¼ odd; n ¼ even;

v

and

½Q  n ð1Þ

¼

8 <

2b ð1bÞ2

:

2b ð1þbÞ2

n ¼ odd; n ¼ even;

so that the four conditions in Theorem 10 are all satisfied. Hence, these rational Gauss–Radau quadrature rules exist, and we can use Theorem 19 to compute their nodes and weights. Tables 1–4 then show the absolute error of the approximation for several values of n and a. To get an idea of the accuracy of the computed nodes and weights, we included the absolute error

n ðju0 j2 Þ of the approximation of Jw ðju0 j2 Þ ¼ 1 too, as well as the difference in absolute value between the fixed node xn;1 and ½comp

the nearest computed node xn;j

. The numerical results clearly show that the quadrature rules are exact for every f 2 Rn;n2 ,

but not for f 2 Rn1;n1 n Rn1;n2 . Example 2. For the second and last example, we consider the function

f ðxÞ ¼

1

p

( sin

)

1

ðx2  JðbÞ2 Þðx2  JðibÞ2 Þ

;

b 2 ð0; 1Þ;

ð25Þ j

which is analogous to the one in Example 5.4 from [24]. This function has essential singularities in x ¼ Jði bÞ; j ¼ 0; . . . ; 3. For b very close to one, this function is extremely oscillating near the center and near the endpoints of the interval I. Since an essential singularity can be viewed as a pole of infinity multiplicity, this suggests taking k1

A2 ¼ fak ¼ Jði

bÞ; k ¼ 1; 2; . . .g:

We then numerical approximate J w ðf Þ by means of an nth rational Gauss-type (Gaussian, Gauss–Radau and/or Gauss–Lobatto) quadrature rule J n ðf Þ. Note that the rational Gaussian quadrature rule only exists for n odd, while the rational Gauss–Radau and Gauss–Lobatto quadrature rules can only exist for n even. The coefficients (4) are now given by

Table 1 Absolute error in the rational Gauss–Radau quadrature formulas with fixed node in 1 for the estimation of J w ðuk ucl Þ.

a

n

ðk; lÞ

Error

n ðju0 j2 Þ

n o ½comp  1j min16j6n jxn;j

2

8

(6, 6) (7, 6) (8, 6) (7, 7) (30, 30) (31, 30) (32, 30) (31, 31) (126, 126) (127, 126) (128, 126) (127, 127)

3.1107e015 6.7671e016 3.7793e015 5.0000e001 6.5060e014 3.4982e014 6.9833e014 5.0000e001 8.6710e014 2.4629e013 1.6776e013 5.0000e001

1.1102e015

4.4409e016

2.2204e016

6.6613e016

8.6597e015

2.2204e016

(6, 6) (7, 6) (8, 6) (7, 7) (30, 30) (31, 30) (32, 30) (31, 31) (126, 126) (127, 126) (128, 126) (127, 127)

2.8869e015 5.9539e015 2.9089e014 9.0909e001 2.5824e013 4.6078e014 3.0015e013 9.0909e001 5.8576e013 5.0504e016 7.8601e013 9.0909e001

6.6613e016

3.3307e016

5.3291e015

4.4409e016

6.2728e014

4.4409e016

32

128

1.1

8

32

128

10318

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

Table 2 Absolute error in the rational Gauss–Radau quadrature formulas with fixed node in 1 for the estimation of J w ðuk ucl Þ.

a

n

1.001

8

32

128

ðk; lÞ

Error

n ðju0 j2 Þ

n o ½comp  1j min16j6n jxn;j

(6, 6) (7, 6) (8, 6) (7, 7) (30, 30) (31, 30) (32, 30) (31, 31) (126, 126) (127, 126) (128, 126) (127, 127)

2.6725e012 6.9373e014 3.7845e012 9.9900e001 1.9464e011 3.3857e013 2.2508e011 9.9900e001 5.3282e010 1.1467e011 5.2787e010 9.9900e001

8.9040e014

8.8818e016

1.8097e013

4.4409e016

5.8442e013

2.2204e015

Table 3 Absolute error in the rational Gauss–Radau quadrature formulas with fixed node in 1 for the estimation of Jw ðuk ucl Þ.

a 2

n 7

31

127

ðk; lÞ

Error

n ðju0 j2 Þ

n o ½comp þ 1j min16j6n jxn;j

(5, 5) (6, 5) (7, 5) (6, 6) (29, 29) (30, 29) (31, 29) (30, 30) (125, 125) (126, 125) (127, 125) (126, 126)

2.1101e015 6.9117e015 9.2753e015 5.0000e001 7.5829e014 3.8533e014 6.8099e014 5.0000e001 3.2041e013 7.0064e014 3.1795e013 5.0000e001

6.6613e016

4.4409e016

3.7748e015

4.4409e016

5.5511e015

6.6613e016

Table 4 Absolute error in the rational Gauss–Radau quadrature formulas with fixed node in 1 for the estimation of Jw ðuk ucl Þ.

a 1.1

n 7

31

127

1.001

7

31

127

ðk; lÞ

Error

n ðju0 j2 Þ

n o ½comp þ 1j min16j6n jxn;j

(5, 5) (6, 5) (7, 5) (6, 6) (29, 29) (30, 29) (31, 29) (30, 30) (125, 125) (126, 125) (127, 125) (126, 126)

1.8874e014 7.4455e015 1.0161e014 9.0909e001 4.8961e013 1.2114e013 5.6687e013 9.0909e001 1.4577e012 2.8959e013 1.6543e012 9.0909e001

3.9968e015

2.2204e016

2.4425e015

8.8818e016

7.9492e014

8.8818e016

(5, 5) (6, 5) (7, 5) (6, 6) (29, 29) (30, 29) (31, 29) (30, 30) (125, 125) (126, 125) (127, 125) (126, 126)

5.5702e012 1.3406e013 4.2642e012 9.9900e001 3.0297e011 6.5270e013 2.2020e011 9.9900e001 3.4639e010 8.8512e012 3.5332e010 9.9900e001

6.8834e014

2.2204e016

1.1113e013

1.1102e016

8.6042e013

1.5543e015

10319

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

Gaussian Gauss−Radau Gauss−Lobatto

0

10

−2

10

−4

Relative error

10

−6

10

−8

10

−10

10

−12

10

−14

10

0

10

20

30

40

50

60

Number of interpolation points Fig. 1. Relative error in the rational Gauss-type quadrature formulas for the estimation of J w ðf Þ, where f is given by (25).

½u

½u

ck1 ¼ ak ¼

8 1þb2 > ; < pffiffiffiffiffiffiffiffiffiffiffiffi 2

k ¼ 1;

> :

k>1

2ð1b Þ

1þb2 ; 2ð1ð1Þk ib2 Þ

and

½u

bk1 ¼

8 b; > > > bð1þb2 Þ½2þið1b2 Þ > > < 2 2

k¼1 k ¼ 2;

2ð1b Þð1ib Þ

> 0; k ¼ 2m þ 1; > > > 2 > ÞÞ : ð1Þm bð1þb2 Þ½ð1þb2 Þþið1b ; k ¼ 2m þ 2; ð1b2 Þð1ib2 Þ2

m ¼ 1; 2; . . .

Further, since the last pole is complex for n even, the nth rational Gauss–Radau quadrature rule is unique (it is the nth rational Gaussian quadrature rule with the last two poles interchanged), so that (see Theorem 10)

A2m;1 ¼ ð1Þm ½Q 

b½ð1  b2 Þ þ ið1 þ b2 Þ 2

2

ð1 þ ib Þð1  b Þ

and

v½Q2m ¼

2b2 ð1  b2 Þ2

;

m ¼ 1; 2; . . .

For the nth rational Gauss–Lobatto quadrature rule, with n ¼ 2m P 4, we consider the case of fixed nodes in xn;1 ¼ xn;2 ¼ 1. We then have that ½P

A2m;1 ð1; 1Þ ¼ 0;

½P

A2m;2 ð1; 1Þ ¼ 1;

and

v½P m ¼ 2; 3; . . . 2m ð1; 1Þ ¼ 1;

½P

It can be verified that the coefficients A2m;j ð1; 1Þ; j ¼ 1; 2, satisfy the five conditions in Theorem 14, so that these rational Gauss–Lobatto quadrature rules exist. Thus, we can again use Theorem 19 to compute the nodes and weights in the rational Gauss-type quadrature rules. For b ¼ 45 we find with the aid of MAPLEÒ 106 that

J w ðf Þ  0:2888276658908954: Fig. 1 then shows the relative error of the approximation in terms of the number of interpolation points n. Note that for n ¼ 4m þ 1; m ¼ 1; 2; . . ., the rational functions in Ln1 ¼ L4m that interpolate the function f ðxÞ in the nth rational Gaussian nodes, have poles in the four different singularities of f ðxÞ, each with the same multiplicity m. This may explain the oscillating behavior of the relative errors for the rational Gaussian quadrature formulas as a function of n, with better results for the ð4m þ 1Þth rational Gaussian quadrature formulas compared with the ð4m þ 3Þth rational Gaussian, and ð4m þ 2Þth and ð4m þ 4Þth rational Gauss–Radau quadrature formulas. Further, note that the graph for the rational Gauss–Lobatto quadrature rule (more or less) coincides with the graph for the rational Gaussian quadrature rule when shifting the latter to the right (the ðn þ 1Þth rational Gauss–Lobatto and nth rational Gaussian quadrature formulas also have the same domain of 6

MAPLE is a registered trademark of Waterloo Maple, Inc.

10320

K. Deckers, A. Bultheel / Applied Mathematics and Computation 218 (2012) 10299–10320

validity). Due to this shift, we obtain that the rational Gauss–Lobatto quadrature rule performs a little bit better compared with the rational Gauss–Radau quadrature rule, although the latter has larger domain of validity. Finally, to get an idea of the accuracy of the computed nodes and weights in the rational Gaussian and Gauss–Radau quadrature formulas, we also computed the nodes and weights with the aid of the MATLABÒ function rcheb, as described in [25]. For each n, the difference in absolute value between the nodes and weights computed by means of the (modified) GEP and the nodes and weights computed by means of the MATLABÒ function rcheb was less than 5:22e  015 and 1:32e  013 respectively. On the other hand, for the rational Gauss–Lobatto quadrature formulas the difference in absolute value between the fixed nodes and the nearest computed nodes was less than 2:11e  015, while the absolute error of the approximation of J w ðju0 j2 Þ ¼ 1 was less than 9:66e  014. References [1] A. Bultheel, R. Cruz-Barroso, M. Van Barel, On Gauss-type quadrature formulas with prescribed nodes anywhere on the real line, Calcolo 47 (1) (2010) 21–48. [2] A. Bultheel, P. González-Vera, E. Hendriksen, O. Njåstad, Orthogonal Rational Functions, Cambridge Monographs on Applied and Computational Mathematics, vol.5, Cambridge University Press, Cambridge, 1999. [3] A. Bultheel, P. González-Vera, E. Hendriksen, O. Njåstad, Orthogonal rational functions and quadrature on the real half line, Journal of Complexity 19 (3) (2003) 212–230. [4] K. Deckers, A. Bultheel, Orthogonal rational functions with complex poles: the Favard theorem, Journal of Mathematical Analysis and Applications 356 (2) (2009) 764–768. [5] K. Deckers, A. Bultheel, Recurrence and asymptotics for orthogonal rational functions on an interval, IMA Journal of Numerical Analysis 29 (1) (2009) 1–23. [6] K. Deckers, A. Bultheel, R. Cruz-Barosso, F. Perdomo-Pı´o, Positive rational interpolatory quadrature formulas on the unit circle and the interval, Applied Numerical Mathematics 60 (12) (2010) 1286–1299. }–Lobatto quadrature on the interval and the unit circle respectively, [7] K. Deckers, A. Bultheel, F. Perdomo-Pı´o, Rational Gauss–Radau and rational Szego Jaen Journal on Approximation 3 (1) (2011) 15–66. [8] K. Deckers, A. Bultheel, J. Van Deun, A generalized eigenvalue problem for quasi-orthogonal rational functions, Numerische Mathematik 117 (3) (2011) 463–506. [9] K. Deckers, J. Van Deun, A. Bultheel, Computing rational Gauss–Chebyshev quadrature formulas with complex poles: the algorithm, Advances in Engineering Software 40 (8) (2009) 707–717. [10] K. Deckers, J. Van Deun, A. Bultheel, Rational Gauss–Chebyshev quadrature formulas for complex poles outside [-1, 1], Mathematics of Computation 77 (262) (2008) 967–983. [11] M.M. Dzˇrbašian, A survey on the theory of orthogonal systems and some open problems, in: P. Nevai (Ed.), Orthogonal Polynomials: Theory and Practice, NATO-ASI Series C: Mathematical and Physical Sciences, 294, Kluwer Academic Publishers, Boston, 1990, pp. 135–146. [12] D. Fasino, L. Gemignani, Structured eigenvalue problems for rational gauss quadrature, Numerical Algorithms 45 (1–4) (2007) 195–204. [13] W. Gautschi, Algorithm 793: GQRAT–Gauss quadrature for rational functions, ACM Transactions on Mathematical Software 25 (2) (1999) 213–239. [14] A. Ihsan Hascelik, Gauss quadrature rules for a generalized Hermite weight function, Applied Mathematics and Computation 180 (1) (2006) 86–96. [15] J.R. Illán González, Gaussian rational quadrature formulas for ill-scaled integrands, Journal of Computational and Applied Mathematics 233 (3) (2009) 745–748. [16] D.P. Laurie, Computation of Gauss-type quadrature formulas, Journal of Computational and Applied Mathematics 127 (1–2) (2001) 201–217. [17] G. Mastroianni, G.V. Milovanovic´, Weighted integration of periodic functions on the real line, Applied Mathematics and Computation 128 (2–3) (2002) 365–378. [18] K. Müller, A. Bultheel, Translation of the Russian paper Orthogonal systems of rational functions on the unit circle by M.M. Dzˇrbašian, Technical Report TW253, Department of Computer Science, K.U. Leuven, February 1997. [19] E.A. Rovba, Orthogonal systems of rational functions on the segment and quadratures of Gauss-type, Mathematica Balkanica (N.S.) 13 (1–2) (1999) 187–198. [20] W. Van Assche, I. Vanherwegen, Quadrature formulas based on rational interpolation, Mathematics of Computation 61 (204) (1993) 765–783. [21] J. Van Deun, Eigenvalue problems to compute almost optimal points for rational interpolation with prescribed poles, Numerical Algorithms 45 (1–4) (2007) 89–99. [22] J. Van Deun, A. Bultheel, Orthogonal rational functions and quadrature on an interval, Journal of Computational and Applied Mathematics 153 (1–2) (2003) 487–495. [23] J. Van Deun, A. Bultheel, Orthogonal rational functions on an interval, Technical Report TW322, Department of Computer Science, K.U. Leuven, March 2001. [24] J. Van Deun, A. Bultheel, P. González-Vera, On computing rational Gauss–Chebyshev quadrature formulas, Mathematics of Computation 75 (253) (2006) 307–326. [25] J. Van Deun, K. Deckers, A. Bultheel, J.A.C. Weideman, Algorithm 882: near best fixed pole rational interpolation with applications in spectral methods, ACM Transactions on Mathematical Software 35 (2) (2008) 14:1–14:21.