View of Weighted quadrature formulas for semi-infinite range integrals

(1)

J. Numer. Anal. Approx. Theory, vol. 44 (2015) no. 1, pp. 69–80 ictp.acad.ro/jnaat

WEIGHTED QUADRATURE FORMULAS FOR SEMI-INFINITE RANGE INTEGRALS^∗

GRADIMIR V. MILOVANOVI ´C

Dedicated to prof. Ion P˘av˘aloiu on the occasion of his 75th anniversary

Abstract. Weighted quadrature formulas on the half line (a,+∞),a >0, for non-exponentially decreasing integrands are developed. Suchn-point quadrature rules are exact for all functions of the form x 7→ x⁻²P(x⁻¹), where P is an arbitrary algebraic polynomial of degree at most 2n−1. In particular, quadrature formulas with respect to the weight functionx7→w(x) =x^βlog^mx(0≤β <1, m∈N0) are considered and several numerical examples are included.

MSC 2010. 65D30, 65D32.

Keywords. Gaussian quadrature rules, nodes, Christoffel numbers, non- exponentially decreasing integrands.

1. INTRODUCTION

In this paper we consider weighted quadrature formulae on the half line (a,+∞),

(1.1)

Z +∞

a

w(x)f(x) dx=

n

X

k=1

Akf(xk) +Rn(f),

whereais a finite real number andx7→w(x) is a given weight function. Such a quadrature formula fora= 0 andw(x) =x^αe^−x,α >−1, is the well known generalized Gauss-Laguerre quadrature rule (cf. [10, p. 325]), which is exact for all algebraic polynomials of degree at most 2n−1, i.e., whenf ∈P2n−1.

Error analysis and convergence of such Gaussian formulas on unbounded intervals (with the classical measures of Laguerre and Hermite) was given in 1928 by Uspensky [17]. Otherwise, the corresponding problems for quadrature rules on finite intervals was studied much earlier by [15], Markov [8], Stieltjes [16], etc. On some new results in this directions see books [5] and [10], including the so-called truncated quadrature rules obtained by ignoring the last part of its nodes (see Mastroianni and Monegato [9]).

∗This paper was supported by the Serbian Ministry of Education, Science and Techno- logical Development (No. #OI 174015).

Serbian Academy of Sciences and Arts, Beograd, Serbia & State University of Novi Pazar, Serbia, e-mail: [email protected].

(2)

Very recently Gautschi [6] has constructed a special logarithmically weighted quadrature formula on (0,+∞), when x7→ (x−1−logx)e^−x. Also, Xu and Milovanovi´c [18] have developed generalized Gaussian quadrature rules of the form (1.1), with x 7→ w(x) = e^−x on (0,+∞), which are exact on the set of basis functions {1,logx, x, xlogx, . . . , xⁿ⁻¹, xⁿ⁻¹logx}. In the other words, these rules are exact for each f(x) = p(x) +q(x) logx, where p, q ∈ Pn−1, so that they can calculate integrals with a sufficient accuracy, regardless of whether their integrands contain a logarithmic singularity, or they do not.

For a similar approach for integrals on the finite intervals see [12] and [14].

On the other side, a large number of integrals of the form ^R_a^+∞F(x) dx which appear in applications do not have exponentially decreasing integrands F(x), and in such cases Gauss-Laguerre quadrature rules are notoriously poor (see Evans [2]). As a starting simple example, Evans [2] has consideredF(x) = 1/(x²+ 0.25), where the convergence of the corresponding integral depends on the 1/x² term for largex. He has proposed a quadrature method based on the set of basis functions {1/x^k} and demonstrated its effectiveness on a series of numerical examples.

In this paper we develop a general approach for constructing a class of n- point generalized quadrature rules (1.1) of Gaussian type on (a,+∞), a >0, which are exact for all functions of the form x 7→ x⁻²P(x⁻¹), where P is an arbitrary algebraic polynomial of degree at most 2n−1. In particular, we consider quadrature formulas with respect to the weight functionx7→w(x) = x^βlog^mx(0≤β <1,m∈N0), which reduces to the constant weight forβ= 0 and m = 0. In order to show the efficiency of the obtained quadrature rules we present a few numerical examples.

2. GENERALIZED WEIGHTED GAUSSIAN RULES

Supposea >0, as well as that the weight function x7→w(x) on (a,+∞) is such that

(2.1) 0<

Z +∞

a

w(x)

x² dx <+∞.

Following Evans [2], we develop a general approach for constructing generalized Gaussian quadrature formulas of the form (1.1). In the cases of integrals on (α,+∞), whenα < a, we simply take

Z +∞

α

w(x)f(x) dx= Z a

α

w(x)f(x) dx+ Z +∞

a

w(x)f(x) dx

and apply to first integral on the right hand side some of rules for calculating integrals on the finite intervals. Also, we mention here that a faster convergence of the corresponding quadrature process can be achieved by taking a greater value of a.

(3)

Thus, the basic idea is to construct a quadrature formula of the form (1.1), which is exact for all functions of the form

x7→ 1 x²P_m1

x

, m= 0,1, . . . ,2n−1,

whereP_m(t) are arbitrary selected algebraic polynomials intof degreem, i.e., (2.2)

Z +∞

a

w(x) 1 x²P_m1

x dx=

n

X

k=1

A_k x²_kP_m 1

x_k

, m= 0,1, . . . ,2n−1.

Remark. Because of linearity, it is easy to see that this system of 2nnon- linear equations inx_k andA_k,k= 1, . . . , n, is equivalent to the corresponding system with monomials, i.e., when Pm(x) =x^m,m= 0,1, . . . ,2n−1.

On the other side we consider the Gauss-Christoffel quadrature formula with respect to the weight functiont7→w(1/t) on (0,1/a), i.e.,

(2.3)

Z 1/a 0

w1 t

g(t) dt=

n

X

k=1

Bkg(τk) +R^G_n(g),

whereτ_kandB_kare its nodes and Christoffel numbers, respectively, andR^G_n(g) is the corresponding remainder term. According to (2.1), such quadrature formulas exist uniquely, because the all momentsµ_k=^R₀^1/aw(1/t)t^kdt,k≥0, exist and µ₀ >0.

It is known that the nodesτ_k in (2.3) are eigenvalues of the following sym- metric tridiagonal Jacobi matrix (cf. [10, pp. 325–328])

(2.4) Jn=







α₀ √

β₁ O

√β₁ α₁ √ β₂

√β₂ α₂ . ..

. .. . .. ^pβn−1

O ^pβn−1 αn−1





 ,

whereα_k and β_k are coefficients in the three-term recurrence relation πk+1(t) = (t−αk)πk(t)−βkπk−1(t), k= 0,1, . . . , (2.5)

π₀(t) = 1, π−1(t) = 0,

for the (monic) polynomials {π_k}_k∈_N₀ orthogonal with respect to the inner product

(2.6) (p, q) =

Z 1/a 0

w1 t

p(t)q(t) dt.

In fact,πn(t) = (t−τ1)· · ·(t−τn).

The weight coefficientsB_k in (2.3) are given by B_k =β₀v_k,1² , k= 1, . . . , n,

(4)

where v_k,1 is the first component of the eigenvector v_k (= [v_k,1 . . . v_k,n]^T) corresponding to the eigenvalue τk and normalized such that v^T_kv_k = 1, and β₀ =µ₀ =^R₀^1/aw(1/t) dt.

The most popular method for solving this eigenvalue problem is the Golub- Welsch procedure, obtained by a simplification of QR algorithm [7]. This procedure is implemented in several packages including the most known ORTPOL given by Gautschi [4].

As we can see from (2.4), for constructing Gauss–Christoffel quadratures (2.3) for any number of nodes less than or equal to n, we need the first n recursion coefficientsα_k andβ_k,k= 0,1, . . . , n−1, in (2.5).

In general, the recursion coefficients are known explicitly only for some narrow classes of orthogonal polynomialsc(e.g. for the classical orthogonal polynomials). In the case of the so-called strongly non-classical polynomi- als, these recursion coefficients must be constructed numerically (cf. [3], [5], [10, pp. 159–166]). However, recent progress in symbolic computation and variable-precision arithmetic today makes it possible to generate the recursive coefficients in (2.5) directly by using the original Chebyshev method of moments. Respectively symbolic/variable-precision software for orthogonal polynomials and Gaussian (and similar) quadratures is available. OurMath- ematica package OrthogonalPolynomials (see [1] and [13]), is download- able from the web site http://www.mi.sanu.ac.rs/˜gvm/. Also, there is Gautschi’s software inMatlab(packages OPQand SOPQ).

Now, we can give our main result:

Theorem2.1. Let x7→w(x) be a weight function on (a,+∞),a >0, such that the condition (2.1) holds. Assume also that τ_k and B_k, k = 1, . . . , n, are nodes and Christoffel numbers of the Gaussian quadrature formula (2.3), respectively. Then there exists the generalized Gaussian quadrature formula

(2.7)

Z +∞

a

w(x)f(x) dx=

n

X

k=1

A_kf(x_k) +R_n(f), with

(2.8) xk= 1

τ_k, Ak= Bk

τ_k² >0, k= 1, . . . , n,

which is exact for all functions of the form f(x) = x⁻²P(x⁻¹), where P ∈ P2n−1.

The remainder term in this quadrature rule can be expressed in the following form R_n(f) =R^G_n(g), where g(t) =t⁻²f(t⁻¹).

Proof. We start with the system of 2n nonlinear equations (2.2), whose solution determines the parameters of the quadrature formula (2.7). Our aim is to prove that this solution uniquely exists.

First, we take the sequence of orthogonal polynomials {π_m}²ⁿ⁻¹_m=0 in the system (2.2) and then by a simple change of variablesx= 1/tin the integral

(5)

on the left hand side we obtain Z +∞

a

w(x) 1 x²πm

1 x

dx = Z 1/a

0

w1 t

πm(t) dt

= (π₀, π_m)

= µ0δ0,m,

where the inner product is defined by (2.6) and δk,mis Kronecker’s delta.

Evidently, this leads to the system of equations (2.9)

n

X

k=1

A_k 1 x²_kπm

1 x_k

=µ₀δ0,m, m= 0,1, . . . ,2n−1,

but, by an application of the Gaussian rule (2.3), it gives also another system of equations

(2.10)

n

X

k=1

B_kπm(τ_k) =µ0δ0,m, m= 0,1, . . . ,2n−1,

because R_n^G(g) = 0 for each g ∈ P2n−1. The last system has the unique solution, and it represents the parameters τk and Bk, k = 1, . . . , n, of the Gaussian quadrature (2.3).

Since the systems of equations (2.9) and (2.10) are equivalent, the statement

of this theorem follows directly.

3. SPECIAL CASES AND NUMERICAL EXAMPLES

In this section we consider special cases of quadrature formulas with respect to the weight function x7→ w(x) =x^βlog^mx, where 0≤β < 1 and m∈N0. For β = 0 and m = 0, it reduces to the constant weight w(x) = 1. In order to show the efficiency of the obtained quadrature formulas we present a few numerical examples.

We start this section with the weight functionx7→w(x) =x^β, 0≤β <1.

The condition (2.1) is satisfied, because Z +∞

a

w(x)

x² dx= a^β−1 1−β.

Here we consider only the caseβ = 0, i.e., whenw(x) = 1. Since Z +∞

a

f(x) dx=a Z +∞

1

f(ax) dx,

we see that for this important case the following statement holds.

Proposition 3.1. Let (3.1)

Z +∞

a

f(x) dx=

n

X

k=1

A_k(a)f(x_k(a)) +R_n(f;a), a >0,

(6)

be a generalized Gaussian quadrature (2.7) (with the constant weight function w(x) = 1). Then

A_k(a) =aA_k(1) and x_k(a) =ax_k(1), k= 1, . . . , n.

This means that it is enough to know only quadrature parameters fora= 1.

These parameters can be obtained directly using (2.8) and Gauss-Legendre parameters τ_k and B_k for transformed interval (0,1).

Recursive coefficients in (2.5), in this case for translated monic Legendre polynomials, are

α_k= 1

2, k≥0, β₀ = 1, β_k= k²

4(4k²−1), k≥1.

Otherwise, it can be obtained using ourMathematicaPackageOrthogonal- Polynomialsin symbolic form (see [1] and [13]). For example, if we need the first forty recurrence coefficients, then we start with the first eighty moments µ_k = 1/(k+ 1), k = 0,1, . . . ,79, and then we use the standard Chebyshev algorithm (cf. [10, 160–162]:

<< orthogonalPolynomials‘

mom=Table[1/(k+1), {k,0,79}];

{al,be} = aChebyshevAlgorithm[mom, Algorithm -> Symbolic]

These recursive coefficients enable us to construct quadrature formulas (2.3) for any number of nodes up to 40.

However, in this Legendre case (translated to (0,1)) we can directly use aGaussianNodesWeightsroutine to construct nodes and weights in the Gauss- ian quadrature formula (2.3), as well as ones in the quadrature formula (2.7):

transLeg[n_] := (aGaussianNodesWeights[n, {aLegendre}, WorkingPrecision -> 70, Precision -> 65] + {1,0})/2;

parQF = Table[transLeg[n], {n,2,40,2}];

For[m = 1, m < 21, m++,

parQF[[m]][[2]] = parQF[[m]][[2]]/parQF[[m]][[1]]ˆ2;

parQF[[m]][[1]] = 1/parQF[[m]][[1]];]

Thus, in this way for a = 1, we obtain quadrature parameters x_k and A_k for each n= 2(2)40.

Example 3.2. In order to show the efficiency of our quadrature rule (2.7) we apply it to the integral

J(a;c) = Z +∞

a

1

(x−2)²+c² dx= 1 2c

π−2 arctan a−2 c

,

for different values of a >0 and c > 0. In Figure 3.1 we present graphics of the function

x7→f(x;c) = 1 (x−2)²+c²

(7)

forc= ¹₈,¹₄,¹₂, and 1, as well as the corresponding graphics of the exact values of this integralJ(a;c) (right).

In order to test the quadrature formula (3.1), we apply it to J(a; 1) for a= ¹₂, 1, 2, 3, 4, and 8, when n= 2(2)40.

� � � � �

��

�

� � � � � �

�

��

�(���)

Fig. 3.1. The functionx7→f(x;c) (left) and the integrala7→J(a;c) (right) forc= ¹₈ (blue line),c=¹₄ (black line),c=¹₂ (brown line), andc= 1 (red line).

Relative errors in the quadrature sums Qn(f(·;c);a) =

n

X

k=1

Ak(a)f(xk(a);c), defined by

err_n(f(·;c);a) =Q_n(f(·;c);a)−J(a;c) J(a;c)

,

are displayed in Figure 3.2 in a log-scale. Numerical results show that the

● ● ● ● ● ● ● ● ● ● ● ●

● ● ● ● ● ● ● ●

■ ■ ■ ■ ■

■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■

◆ ◆ ◆

◆ ◆ ◆ ◆

◆ ◆ ◆

◆ ◆ ◆ ◆

◆ ◆ ◆

▲ ▲

▼ ▼

▼

▼ ▼

▼

▼ ▼

▼

○

○ ○

○

10 20 ○30 40

10^-60 10^-40 10^-20

1 10 20 30 40

● a=1/2

■ a=1

◆ a=2

▲ a=3

▼ a=4

○ a=8

Fig. 3.2. Relative errors errn(f(·; 1);a) in quadrature sumsQn(f(·; 1);a).

(8)

convergence is much faster if the parameter a is larger. For example, if a= 2, then for n = 10(10)40, the relative errors are 1.71×10⁻⁷, 1.83×10⁻¹⁴, 1.91×10⁻²¹, 1.94×10⁻²⁸, respectively, while the corresponding errors for a= 4 are 5.52×10⁻¹⁵, 1.21×10⁻²⁹, 1.40×10⁻⁴⁴, 1.44×10⁻⁵⁹.

Otherwise, this integrand f(x;c) has poles at the points 2±ic, which are approaching the real line when c tends to zero. In this case, for small values of a (near 2 or less than 2), the convergence of the quadrature process slows down considerably, because of a strong influence of these singularities.

This effect can be seen from Table 3.1, where quadrature approximations and corresponding relative errors are presented for (a, c) = (1,¹₄), (²¹₁₀,10⁻⁶), and (4,10⁻⁶). In order to save space, in last case only relative errors are given.

Digits in error are underlined, and numbers in parenthesis indicate the decimal exponents.

Notice that the integral J(a; 0) for a≤2 does not exist.

Finally, the last column shows that for (a, c) = (4,10⁻⁶), the convergence of the quadrature rule (3.1) is very fast.

n (a, c) = (1,1/4) (a, c) = (21/10,10⁻⁶) (a, c) = (4,10⁻⁶) 2 2.83088 7.56(−1) 4.21706255691703 5.78(−1) 5.92(−3) 4 5.38719 5.35(−1) 8.01223217799471 1.99(−1) 9.70(−6) 6 7.41379 3.60(−1) 9.47887835712778 5.21(−2) 1.24(−8) 8 8.88711 2.33(−1) 9.88043864297441 1.20(−2) 1.42(−11) 10 9.89102 1.46(−1) 9.97447558340612 2.55(−3) 1.53(−14) 20 11.45438 1.14(−2) 9.99999276505451 7.23(−7) 1.47(−29) 30 11.57808 7.23(−4) 9.99999999813998 1.53(−10) 1.08(−44) 40 11.58606 3.41(−5) 9.99999999966638 2.86(−14) 6.99(−60)

Table 3.1. Quadrature sumsQn(f(·;c);a) and their relative errors errn(f(·;c);a) for integralsJ(a;c).

In the sequel we consider quadrature rules with respect to the weight function x 7→ w(x) = x^βlogx, 0 ≤ β < 1. Here we suppose that a ≥ 1. The condition (2.1) is satisfied, because

(3.2) 0<

Z +∞

a

w(x)

x² dx= a^β−1

(1−β)²[1 + (1−β) loga]. In this case, the moments

µ_k= Z 1/a

0

w(1/t)t^kdt= Z 1/a

0

t^k−βlog1 t dt can be expressed in the form

(3.3) µ_k= a^β−k−1[(k+ 1−β) loga+ 1]

(k+ 1−β)² , k≥0.

Taking the first one hundred moments (mom) (e.g. for a= 1 and β = 1/4) and using Mathematica Package OrthogonalPolynomials, we can get the

(9)

first fifty recurrence coefficients α_k and β_k (denoted by {al1,be1}) in the three-term recurrence relation (2.5) in a symbolic form

mom=Table[(aˆ(-1+b-k)(1+(1-b+k)Log[a]))/(1-b+k)ˆ2, {k,0,99}];

mom1=mom/. {a->1, b->1/4}

{al1,be1} = aChebyshevAlgorithm[mom, Algorithm -> Symbolic]

For example, first four coefficients are α₀= 9

49, α₁= 209897

452025, α₂= 6582284926939

13538179995075, α₃= 7618613698603068100869609 15464687102113919816429449 and

β₀=16

9 , β₁= 11808

290521, β₂= 213147564896

3717280400625, β₃= 421267942813254097088 6997413354065613077481.

Example 3.3. As a test example we consider the function x7→f(x) = 1

(x+ 1)² and integral (see [2])

I(f;a) = Z +∞

a

x^1/4logx (x+ 1)² dx

= 1

36a⁷⁴(a+ 1)

−9(a+ 1)Φ

−1 a,2,7

4

+4a

3(a+ 1)(loga+ 4)2F1

3 4,1;7

4;−1 a

+ 4a+ 9aloga+ 4

,

where Φ and ₂F₁ are the Lerch transcendent and Gauss hypergeometric function, defined by

Φ(z, s, a) =

+∞

X

k=0

z^k

(k+a)^s and ₂F1(a, b;c;z) =

+∞

X

k=0

(a)_k(b)_k (c)_k

z^k k!, respectively, and (a)_k=a(a+ 1). . .(a+k−1) is the Pochhammer symbol.

We consider this integral for two values of the lower bound: a = 1 and a= e, i.e.,

I(f; 1) = 1.359 743 280 976 008 95. . . and I(f; e) = 1.228 976 186 680 372 55. . . .

Applying Gauss-Laguerre rule toI(f; 1) (translated from (1,+∞) to (1,+∞)) gives poor results. Relative errors in the corresponding Gauss-Laguerre quadrature sums are presented in Table 3.2. As we can see only two two decimal digits are true in quadrature sum with 2048 nodes!

Now, we apply our quadrature formula (2.7), with parameters given by (2.8), to I(f; 1) and I(f; e), with only n= 2(2)12 nodes. The relative errors

(10)

n= 2 n= 8 n= 32 n= 128 n= 512 n= 2048 6.72(−1) 3.60(−1) 1.64(−1) 7.00(−2) 2.90(−2) 1.18(−2)

Table 3.2. Relative errors in Gauss-Laguerre quadrature sums withn= 2,8,32,128,512 and 2048 nodes.

in the quadrature sums Qn(f;a) =^Pⁿ_k=1Akf(xk), err_n(f;a) =Q_n(f;a)−I(f;a)

I(f;a) , are presented in Table 3.3.

a n= 2 n= 4 n= 6 n= 8 n= 10 n= 12

1 2.94(−3) 4.24(−6) 5.15(−9) 5.72(−12) 4.74(−13) 7.07(−13) e 2.40(−4) 1.64(−8) 8.91(−13) 8.83(−14) 5.31(−14) 3.80(−14) e² 7.18(−6) 1.28(−11) 3.10(−14)

Table 3.3. Relative errors errn(f;a) in quadrature sumsQn(f;a) for different number of nodesnand three values ofa(= 1,e, and e²).

As we can see, the convergence is faster whenais bigger. In the third line of the same table we also present the corresponding relative errors whena= e² and n= 2,4, and 6.

Finally, we mention that this approach can be applied also in the case of the weight functions

w(x) =wm(x) =x^βlog^mx, 0≤β <1, m= 2,3, . . . , on the interval (a,+∞), witha≥1.

The condition (2.1) forUm=^R_a^+∞x⁻²wm(x) dxis also satisfied, because U_m = 1

1−β

hmUm−1+a^β−1log^maⁱ, m= 2,3, . . . , whereU₁ is given in (3.2). The corresponding moments

µ^[m]_k = Z 1/a

0

t^k−βlog^m1

tdt, k= 0,1, . . . , can be expressed recursively in terms of the momentsµ^[m−1]_k ,

µ^[m]_k = 1 k+ 1−β

mµ^[m−1]_k +a^β−k−1log^ma, m= 2,3, . . . ,

where the momentsµ^[1]_k (≡µ_k) are given by (3.3). For example, form= 2 we get

µ^[2]_k =

a^β−k−1^h(k+ 1−β)²log²a+ 2(k+ 1−β) loga+ 2ⁱ

(k+ 1−β)³ , k≥0.

(11)

The corresponding recursive coefficients, for example for β = 0 and a= 1, are

α₀=1

8, α₁= 115

296, α₂= 28200187

62721512, α₃= 28003451041760695 59414538084233528, . . . and

β0= 2, β1= 37

1728, β2= 211897

4620375, β3= 945381680572419

17600932734728000, . . . .

Example 3.4. We consider the function x 7→ 1/(1 +x²) and the corresponding weighted integral over (a,+∞)

I(f;a) = Z +∞

a

log²x 1 +x²dx, for two values ofa (= 1 and = e), for which

I(f; 1) = 1.937 892 292 518 738 760 967 269 691 69. . . and

I(f; e) = 1.809 886 879 397 869 426 020 164 472 46. . . .

Applying our quadrature formula (2.7), with parameters given by (2.8), to I(f; 1) andI(f; e), with n= 2(2)12 nodes, we get the corresponding quadrature approximations Qn(f;a) with the relative errors errn(f;a) presented in Table 3.4.

a n= 2 n= 4 n= 6 n= 8 n= 10 n= 12

1 1.66(−4) 1.31(−6) 1.98(−10) 5.73(−12) 2.08(−15) 2.56(−17) e 5.33(−5) 5.04(−10) 1.86(−13) 2.05(−17) 1.22(−21) 3.30(−26)

Table 3.4. Relative errors errn(f;a) in quadrature sumsQn(f;a) for different number of nodesnand two values ofa(= 1 and = e).

Here also we can note a faster convergence when ahas a larger value.

REFERENCES

[1] A.S. Cvetkovi´candG.V. Milovanovi´c,The Mathematica Package “OrthogonalPoly- nomials”, Facta Univ. Ser. Math. Inform.,19, pp. 17–36, 2004.

[2] G.A. Evans,Some new thoughts on Gauss-Laguerre quadrature, Int. J. Comput. Math.

82, pp. 721–730, 2005.

[3] W. Gautschi,On generating orthogonal polynomials, SIAM J. Sci. Statist. Comput., 3, pp. 289–317, 1982.

[4] W. Gautschi,Algorithm 726: ORTHPOL – A package of routines for generating or- thogonal polynomials and Gauss-type quadrature rules, ACM Trans. Math. Software,20, pp. 21–62, 1994.

[5] W. Gautschi, Orthogonal Polynomials: Computation and Approximation, Clarendon Press, Oxford, 2004.

[6] W. Gautschi,Gauss quadrature routines for two classes of logarithmic weight functions, Numer. Algor.55, pp. 265–277, 2010.

(12)

[7] G. GolubandJ.H. Welsch,Calculation of Gauss quadrature rules, Math. Comp.23, pp. 221–230, 1969.

[8] A. Markov, Sur la méthode de Gauss pour le calcul approché des intégrales, Math.

Ann.,25, pp. 427–432, 1885.

[9] G. Mastroianni and G. Monegato, Truncated quadrature rules over (0,∞) and Nystr¨om-type methods, SIAM J. Numer. Anal.41, pp. 1870–1892, 2003.

[10] G. MastroianniandG.V. Milovanovi´c,Interpolation Processes – Basic Theory and Applications, Springer-Verlag, Berlin – Heidelberg – New York, 2008.

[11] G.V. Milovanovi´c,Construction and applications of Gaussian quadratures with non- classical and exotic weight function, Stud. Univ. Babe¸s-Bolyai Math.,60, pp. 211–233, 2015.

[12] G.V. Milovanovi´c, Generalized Gaussian quadratures for integrals with logarithmic singularity, FILOMAT (to appear).

[13] G.V. Milovanovi´candA.S. Cvetkovi´c,Special classes of orthogonal polynomials and corresponding quadratures of Gaussian type, Math. Balkanica,26, pp. 169–184, 2012.

[14] G.V. Milovanović, T.S. IgićandD. Turnić,Generalized quadrature rules of Gauss- ian type for numerical evaluation of singular integrals, J. Comput. Appl. Math.278, pp.

306–25, 2015.

[15] C. Posse,Sur les quadratures, Nouv. Ann. Math. (2)14, pp. 49–62, 1875.

[16] T.J. Stieltjes, Quelques recherches sur la théorie des quadratures dites mécaniques, Ann. Sci. Éc. Norm. Paris, Sér. 2,1, pp. 409–426, 1884.

[17] J.V. Uspensky,On the convergence of quadrature formulas related to an infinite in- terval, Trans. Amer. Math. Soc.,30, pp. 542–559, 1928.

[18] Zhenhua Xuand G.V. Milovanovi´c,Efficient method for the computation of oscil- latory Bessel transform and Bessel Hilbert transform(submitted).

Received by the editors: November 6, 2015.