TitleI Chebyshev series method for computing weighted

quadrature formulas
E. Berriochoa Esnaola, A. Cachafeiro López, J. R. Illán-González, and E.
Martı́nez-Brey
Departamento de Matemática Aplicada I, Universidad de Vigo, Spain

Abstract
In this paper we study convergence and computation of interpolatory quadra-
ture formulas with respect to a wide variety of weight functions. The main
goal is to evaluate accurately a definite integral, whose mass is highly concen-
trated near some points. The numerical implementation of this approach is
based on the calculation of Chebyshev series and some integration formulas
which are exact for polynomials. In terms of accuracy, the proposed method
can be compared with rational Gauss quadrature formula.
Keywords: Interpolatory quadrature formula, rational quadrature formula,
ill-scaled integrand, difficult pole, modified weight function, Chebyshev
series
2010 MSC: 65D32, 41A55

1. Introduction
A common problem is the evaluation of a definite integral over a bounded
interval [a, b]. The solution of this problem can be approached in terms of a
product integration rule, specially when the integrand has singularities. The
procedure is to rewrite the integrand as F (x)W (x), where F is the integrand
and W is a non-negative weight function.

I
The work of A.C.L., E.B.E. and J.R.I.G. was supported by a research grant from the
Ministry of Sciences and Innovation, Spain, project code MTM 2008-00341
Email addresses: esnaola@uvigo.es (E. Berriochoa Esnaola),
acachafe@uvigo.es (A. Cachafeiro López), ebrey@uvigo.es (and E. Martı́nez-Brey)

Preprint submitted to Elsevier October 10, 2011

 Suppose that w is a weight function on [a, b], which is closely related to
W in a sense to be clarified. The integral of F W can be approximated by a
weighted quadrature of the form
Z b Xn
F (x)W (x)dx = λn,k (W )F (xn,k ) + En (F ), (1)
a k=1

where Z b
Pn (x)
λn,k (W ) = W (x)dx, (2)
a (x − xn,k )Pn0 (xn,k )
and Pn (x) = nk=1 (x − xn,k ) is the nth (monic) orthogonal polynomial with
Q
respect to w.
In the rest of this article, we will maintain the notation used in (2). The
aim is to indicate explicitly the weight function to which we refer.
To improve the accuracy of the results, a plausible strategy consists in
modifying conveniently the initial factorization, say, F W = f GW , where
G ≥ 0, and the main property of f is to have no singularities.
Roughly speaking, a real-valued function is said to be a difficult function,
if it has integrable singularities, or it ranges over very large and very small
values (poor scaling). A typical case is one in which the function is mero-
morphic on a region V ⊃ [a, b]. In this context, the adjective difficult also
applies to the poles of the integrand that are very close to [a, b].
To tackle the integration of a difficult function, the technique of selecting
G as a proper rational function has been used intensively, i.e. f = qF/p
and G = p/q, where p and q are polynomials whose roots match the difficult
points of F (cf. [5]–[7],[9]). Monegato [9] considers G = 1/q, and calculates
λn,k (W/q) in terms of λn,k (W ). That approach requires the calculation of
Rb
the integrals Vi = a xi W (x)/q(x)dx, i = 0, 1. Monegato’s method is based
on the partial fraction decomposition of 1/q, and depends on the ability to
calculate Vi . The main objective of [9] is to obtain accurate results when the
integrand has difficult poles, and this is achieved by applying a non-linear
algorithm proposed by Gautschi. A variant of Monegato’s work is given in [5]
where it is considered {Gn }, defined as Gn = pn /qn , n ∈ N, where {pn } and
{qn } are two sequences of polynomials such that pn /qn = O(1). Moreover,
the roots of pn coincide with some zeros of the integrand, and the difficult
poles of the integrand are neutralized by applying an appropriate change of
variables. The examples in [5] demonstrate that the zeros of the integrand
can play a role, but the numerical results are similar to those in [9] when

2
pn ≡ 1. A common feature to both papers [5, 9] is that the convergence of
the quadrature rules is guaranteed by the following condition (cf. [10]).
Z b
W (x)2 /w(x) dx < ∞. (3)
a

The purpose of this paper is to present a flexible approach to evaluate

the integral of difficult functions. As a result, they have been derived several
methods whose performance is illustrated mainly by means of numerical ex-
amples. In a sense, some of these methods are simpler than that of Monegato
[9]. Our strategy has been based on constructing interpolatory quadrature
formulas with respect to a class of varying weights.
From here onwards, we will only consider integrals on the interval [−1, 1].
The paper is organized as follows.
In Section 2 we consider a modification of (1–2) given in terms of a se-
quence {Gn W }, which means that the coefficients are given by λn,k (Gm W ).
Then we obtain convergence when the trio ({Gn }, W, w) satisfies a condition
weaker than (3). Section 3 deals with the case G = K/q, where K ∈ L2 (w),
and focuses on the calculation of λn,k (KW/q). Here we assume that G is
approximated by Gm = pm /q, where pm is a partial sum of the Chebyshev √
series expansion of K. The problem reduces to consider w(x) = 1/ 1 − x2
and W = h1 w, where h1 is a polynomial. It is clear that W is a generic
weight that includes the four classical Chebyshev weight functions. It is used
a technique similar to that by Berriochoa et al [1], to prove an exact formula
for λn,k (KW/q). Unlike the method proposed by Clenshaw and Curtis [4],
here we only use the Chebyshev series expansion of a factor K.
The numerical examples of Section 4 allow to appraise the effectiveness
of the proposed approach. To facilitate this, we compare our results with
those obtained by Gautschi [7], Monegato [9] and Deckers et al [2]. Section
5 is devoted to the calculation of the coefficients of the Chebyshev series
expansion, and Section 6 contains some remarks as a conclusion.
Throughout the paper Π denotes the space of all algebraic polynomials,
Πn = {p ∈ Π; deg(p) ≤ n}, and kf k = sup{|f (x)|; x ∈ [−1, 1]}.

2. Convergence of modified quadrature rules

Definition 1. Let W and w be two weight functions on [−1, 1], both with
infinitely many points of increase. Let F be an integrable function with res-
pect to W and assume that F = f G, where f is continuous, and G is a

3
nonnegative function such that G ∈ L1 (W ). Let {Gn } be a sequence of
positive and continuous functions on (−1, 1), such that Gn ∈ L1 (W ), n ∈ N,
and Gn W dx → GW dx, as n → ∞, in the weak ∗ topology of measures. We
say that
Z 1 n
X
F (x)W (x)dx = λn,k (Gm W )f (xn,k ) + En,m (f ) (4)
−1 k=1

is a modified interpolatory quadrature formula (M IQF ) of order (n, m), with

respect to the pair ({Gn }, w), if
Z b
Pn (x)
λn,k (Gm W ) = Gm (x)W (x)dx, (5)
a (x − xn,k )Pn0 (xn,k )
Qn
where Pn (x) = k=1 (x − xn,k ) is the nth (monic) orthogonal polynomial with
respect to w.

As usual, the term En,m (f ) stands for the quadrature error.

Proposition 1. Suppose that (4) is a M IQF with respect to ({Gn }, w), and
let n be fixed. Then
Z 1 n
X
f (x)Gm (x)W (x)dx = λn,k (Gm W )f (xn,k ) (6)
−1 k=1

holds for all f ∈ Πn−1 , and all m ∈ N.

Proof. Just note that λn,k (Gm W ) is the kth coefficient of the interpolatory
rule of polynomial type which corresponds to the weight Gm W .
The construction of (4) depends on the integrand F = f G, which in-
creases the cost of the numerical procedure. In effect, the quadrature coef-
ficients must be calculated each time one changes G. On the other hand,
{Gn } represents a discretization method.
If G = p/q, with p, q ∈ Π, then we say that (4) is a rational formula.
To study the convergence of a M IQF , we assume that m = m(n) is
an increasing sequence of integers. For simplicity, we will write Gn and En ,
instead of Gm(n) and En,m(n) , respectively.
The following result brings us closer to the theory in [10].

4
Proposition 2. Let En be the error functional of a M IQF with respect to
({Gn }, w), with Gn = G, for all n. If condition (3) holds true, then
1. limn En (f ) = 0, and
R1
2. limn nk=1 |λn,k (GW )|f (xn,k ) = −1 f (x)G(x)W (x)dx,
P

for all bounded Riemann integrable functions f ,

Proof. Let M > 0 such that G(x) ≤ M for all x ∈ [−1, 1]. According to [10,
Theorem
R1 1], we only have to check
R that 2
2 2 1
−1
(W (x)G(x)) dx/w(x) ≤ M −1 W (x) dx/w(x) < ∞.
The following result is proved using the technique of [10] (See also [5]).

Pn a M IQF with respect to ({Gn }, w), and

Proposition 3. If (4) represents
2
Gn ∈ L2 (W /w), n ∈ N, then j=1 |λn,j (Gn W )| = O(ρn ), where

1 1/2
Gn (x)2 W (x)2
Z
ρn = dx . (7)
−1 w(x)

Moreover, if supn ρn < ∞, then limn En (f ) = 0 for all f continuous.

Proof. Consider the orthonormal basis of polynomials with respect to w,

and let Sn,m be the mth partial sum of the corresponding Fourier series of
W Gn /w ∈ L2 (w). Taking into account that deg(Pn (x)/(x − xn,j )) = n − 1,
we can write
Z 1
Pn (x)Sn,n−1 (x)
λn,j (Gn W ) = 0
w(x)dx.
−1 Pn (xn,j )(x − xn,j )

Now we use the orthogonality of Pn with respect to the following polynomial

Sn,n−1 (x) − Sn,n−1 (xn,j )

,
x − xn,j

with degree n − 2, to obtain that

λn,j (Gn W ) = Sn,n−1 (xn,j )cn,j , (8)

where cn,j are the coefficients of the Gauss rule associated with w.

5
 1/2 1/2
By writing cn,j = cn,j cn,j and adding the first n coefficients in (8) we can
use the Cauchy-Schwartz inequality to obtain the following estimate

n n
!1/2 n
!1/2
X X X
2
|λn,j (Gn W )| ≤ cn,j Sn,n−1 (xn,j ) cn,j . (9)
j=1 j=1 j=1

2 2
Because deg(Sn,n−1 ) = 2n−2, we can apply to Sn,n−1 the exactness condition
for the Gauss quadrature formula associated with w. The Bessel inequality
allows us to obtain
X n Z 1 1/2
|λn,j (Gn W )| ≤ ρn w(x)dx , (10)
j=1 −1

which proves the firstR part of the proposition.P

1
Now let Ln (f ) = −1 f (x)Gn (x)W (x)dx− nk=1 λn,k (Gn W )f (xn,k ), where
f is continuous on [−1, 1]. It follows that
Z 1 Z 1 1/2
f (x)Gn (x)W (x)dx ≤ ρn kf k w(x)dx .
−1 −1

Now suppose that supn ρn < ∞. From (10), it follows the asymptotical
stability of the quadrature formula. Moreover, {Ln } is uniformly bounded.
From the exactness condition (Proposition 1) one obtains that limn Ln (f ) = 0
for all f continuous on [−1, 1].
R1
Using that |En (f )| ≤ −1 f (G − Gn )W dx + |Ln (f )|, the second part of
the proposition is also proved.
If {Gn } satisfies Gn ≤ C0 , n ∈ N, then the asymptotical stability of the
corresponding M IQF can be obtained by assuming condition (3) (cf. [5]).
Let W (x) = 1/(1 − x2 )δ and 2 η
√ h(x) = (1 − x ) with 3/4 ≤ δ 2< 1 and
δ − 3/4 < η. Let w(x) = 1/ 1 − x2 . It implies that h ∈ L2 (W /w) and
(3) is not fulfilled. If Gn ≤ h, n ∈ N, then supn ρn < ∞, where ρn is given
by (7). We conclude that the assumptions of Proposition 3 are weaker than
those of [10, Theorem 1] and [5, Lemma 1].
As expected, the rate of convergence of the quadrature error depends
heavily on the approximation properties of both functions f and G.

6
Proposition 4. Let En be the nth error functional of a M IQF with respect
1/n
to ({Gn }, w). Let us suppose that lim supn ρn ≤ 1, where ρn is given by (7).
If we also have that lim supn kG − Gn k1/n < 1, and f is analytic on [−1, 1],
then lim supn |En (f )|1/n < 1.
Proof. For every n, let p∗n−1 ∈ Πn−1 be the polynomial of best approximation
to f on [−1, 1], in the uniform norm. We obtain that

|En (f )| ≤ kf k kG − Gn kM1 + kf − p∗n−1 kρn (1 + M2 ) ,

R1 R 1/2
1
where M1 = −1 W (x)dx and M2 = w(x)dx
−1
.
Now, we only have to apply the theorem of Bernstein-Walsh which asserts
that lim supn kf − p∗n−1 k1/n < 1.

3. Calculation of quadrature coefficients

For every p ∈ Π we define pe(z) = z n p((z + 1/z)/2), where n = deg(p).
If p(x) = an xn + an−1 xn−1 + · · · + a0 , with an 6= 0, then deg(e p) = 2n and
n
pe(0) = an /2 6= 0. Besides, pe(ζ) = 0 if and only if pe(1/ζ) = 0, ζ 6= 0.
If p1 ∈ Πn1 and p2 ∈ Πn2 , then pg 1 p2 = pe1 pe2 .
In the sequel, we
√ give a more specific meaning to the previous notation.
Set w(x) = 1/ 1 − x2 and W1 = w/π. Now we restrict our attention
to those weight functions which have the form W = h1 W1 , where h1 is an
algebraic polynomial
R1 non-negative
R 1 on [−1, 1]. Notice that condition (3) holds
true, in effect, −1 W dx/w = −1 h1 wdx/π 2 < ∞.
2 2

Let K be a positive and continuous function on [−1, 1], whose Chebyshev

coefficients are defined as usual, i.e.
Z 1 Z 1
A0 = K(x)T0 (x)W1 (x)dx, Aj = 2 K(x)Tj (x)W1 (x)dx, j = 1, 2, ...
−1 −1

where Tj is the jth Chebyshev orthogonal polynomial of the first kind.

Set G = K/q where q is a polynomial positive on [−1, 1]. The follo-
Pn results are given for M IQF s with respect to ({Gn }, w), where Gn =
wing
j=0 Aj Tj /q, n ∈ N.
Lemma 1 is a consequence of [8, Theorem 5.2].
Lemma 1. If K ∈ L2 (w), then Gn W dx → GW dx, as n → ∞, in the weak
∗
topology of measures.

7
Theorem 1. Let q ∈ Πm \Πm−1 be such that q(x) > 0 for all x ∈ [−1,
Pn1]. Set
h1 /q = p1 + p2 /q, where pi ∈ Πni , i = 1, 2 and n2 < m. Let P = i=0 bi Ti ,
where bi ∈ R, i = 0, ..., n. If K ∈ L2 (w), then
Z 1
2 P (x)K(x)h1 (x)W1 (x)/q(x)dx = (11)
−1
nX
1 +n   X ∞ s  !
pe1 (z)Vn (z) X Vn (z)R(z)
< Aj Res ,0 + Aj Res , zi
j=0
z n1 +1−j j=0 i=0
z n2 +1−m−j

where z1 , · · · , zs are the s distinct

Pn zeros of qe in the unit disk, (s ≤ m), z0 = 0,
R(z) = pe2 (z)/e q (z), Vn (z) = i=0 bi (z + z −i ), and Res(F, z) is the residue of
i

F at z.
Proof. If q ∈ Πm possesses s distinct zeros in R\[−1, 1], then, using the
properties of the transform q → qe, one obtains that pe2 (z)/e q (z) has s distinct
poles zi 6= 0 in the unit disk.
h1 (cos(θ))
Let H(θ) = . The following equalities hold true.
q(cos(θ))
Z 1  Z π ∞
P Kh1 W1 1 X
(x)dx = P (cos(θ))H(θ) Aj cos(jθ)dθ =
−1 q 2π −π j=0
Z π n ∞
1 X X Aj
bs cos(sθ)H(θ) eijθ dθ+
2π −π s=0 j=0
2
Z π n ∞
1 X X Aj
bs cos(sθ)H(θ) eijθ dθ =
2π −π s=0 j=0
2
∞
!
pe1 (z) z m pe2 (z)
Z  X
1 1
< Vn (z) + n2 Aj z j−1 dz .
2 2πi |z|=1 z n1 z qe(z) j=0

From Lemma 1 we obtain that

Z 1  ∞  Z
P Kh1 W1 1 X 1 Vn (z)pe1
(x)dx = < Aj dz+
−1 q 2 j=0
2πi |z|=1 z n1 +1−j

Vn (z)pe2 (z)z m−n2 −1+j

Z 
1
dz .
2πi |z|=1 qe(z)
The proof finishes after applying the residue theorem.

8
 Notice that Vn (z) = σ(z)/z n , where σ ∈ Π2n .
Next we prove three corollaries of Theorem 1.
Corollary 1. Let P = ni=0 bi Ti , where bi ∈ R, i = 0, ..., n. If K ∈ L2 (w)
P
then the following formula takes place.
Z 1 d+n
!!
1 X Vn (z)he1 (z)
P (x)(Kh1 W1 )(x)dx = < Aj Res −j+d+1
,0 , (12)
−1 2 j=0
z

where d = deg(h1 ) and Vn (z) = ni=0 bi (z i + z −i ).

P
In particular,
Z 1 n
!
1 X
P (x)(KW1 )(x)dx = A0 b0 + Aj bj . (13)
−1 2 j=0

Proof. To prove (12), we only have to put H(θ) = h1 (cos(θ)) in Theorem 1.

Equation (13) follows easily by direct calculation of the residues.
Pn
Corollary 2. Let P = i=0 bi Ti , where bi ∈ R, i = 0, ..., n. Let q ∈
Πm \Πm−1 be such that q(x) > 0 for all x ∈ [−1, 1]. If z1 , · · · , zs are the s
distinct zeros of qe in the unit disk, then
Z 1
P (x)h1 (x)W1 (x)/q(x)dx = (14)
−1
  X s  !
1 pe1 (z)Vn (z) Vn (z)pe2 (z)
< Res ,0 + Res , zi ,
2 z n1 +1 i=0
q
e (z)z n2 +1−m

where ni = deg(pi ), i = 1, 2, z0 = 0 and Vn (z) = ni=0 bi (z i + z −i ).

P

Proof. The proof is reduced to put K ≡ 1 in Theorem 1.

Corollary 3. Let Pq ∈ Πm \Πm−1 , m ≥ 1, and q(x) > 0 for all x ∈ [−1, 1].
Suppose that P = ni=0 bi Ti , where bi ∈ R, i = 0, ..., n. If K ∈ L2 (w), then
Z 1 ∞ s  !
X X Vn (z)
P (x)dx/q(x) = < Aj Res 1−m−j
, zk , (15)
−1 j=0 k=0
q
e (z)z

where A2j = −2/(4j 2 − 1), A2j−1 = 0, j = 1, 2, · · · , and A0 = 1. Moreover,

1 , · · · , zs are the s distinct zeros of q
zP e in the unit disk, z0 = 0, Vn (z) =
n i −i
i=0 bi (z + z ) and Res(F, z) is the residue of F at z.

9
Proof.
√ (15) follows from Theorem 1 when p1 ≡ 0, p2 ≡ h1 ≡ 1, K(x) =
π 1 − x2 , and every Aj is replaced by their exact expression (cf. [8]).
If deg(h1 ) is too large, then we recommend the use of equation (12) of
Corollary 1, to facilitate the calculation of coefficients Aj .
The presence of q and K at the same time might appear to be over-
abundant. Actually, we can assign all the benign poles to K, while those
that have proven to be difficult must be zeros of q. Corollary 3 includes
the presence of both K and q, and also provides a closed formula for the
coefficients Aj . This corollary allows us to tackle the problem of computing
R1
−1
F (x)W (x)dx when W ≡ 1 and F = f /q. This F is meromorphic and
the zeros√of q should be difficult poles. In principle, an adverse fact is that
K(x) = 1 − x2 has singularities at ±1, which makes its Chebyshev series
slowly converging. However, the numerical results are acceptably accurate,
as demonstrated by some examples. If we are using (12), (13) or (14), then
no matter what is the speed of convergence of the Chebyshev series, because
these formulas are given in terms of finite sums.
Any numerical implementation of the previous formulas is based on the
following issues: i) the nodes coincide with the roots of Chebyshev polyno-
mials, and ii) the coefficients λn,k (GW1 ), 1 ≤ k ≤ n, have to be calculated
for every difficult factor G. Moreover, the most favorable choice is one in
which G does not depend on n, so that the set {Am }N m=0 is calculated only
once with N > n, using F F T algorithm whose order is O(N log(N )).
If P = ln−1,j ∈ Πn−1 is the jth fundamental polynomial of Lagrange
interpolation, then the equations (11)–(15) correspond to the quadrature
coefficient λn,j (GW ), where G = K/q and W = h1 W1 . Moreover, ln−1,j can
be simulated efficiently by a routine based on its representation as a linear
combination of Chebyshev polynomials of the first kind, that is,
1
ln−1,j = (b0,n,j T0 + b1,n,j T1 + · · · + bn−1,n,j Tn−1 ) , (16)
n
where b0,n,j = 1, and bk,n,j = 2 cos((n − k)π(2j + 1)/(2n)), k = 1, ..., n − 1.
Corollary 1 provides the simplest formulas of this section. Equation (12)
can be expressed in terms of some derivatives of hed at z = 0, and equation
(13) only requires the knowledge of bj = bk,n,j , j = 1, ..., n.
The total number of floating point operations to calculate each bj does
not depend on n. Thus, when using (13), the process of calculating the n

10
coefficients λn,j (GW ) has a computational cost proportional to n2 , provided
that the work of calculating {Aj }nj=0 is considered separately.
The use of (14)–(15) involves the calculation of residues of rational func-
tions at some points of the unit disk. Let p and q be polynomials with
q(z0 ) 6= 0. The following algorithm can be used to estimate the residue of
rk (z) = p(z)/q(z)/(z − z0 )k at z = z0 , for k = 1, 2, · · · , ν.

0
 deg(q) ≥ 1, calculate V = q , A = q(z0 ) and B = V (z0 ),
If
For k = 1, 2, ..., ν, do
 Res(rk , z0 ) = pk−1 (z0 )/A, where p0 = p, and
pk−1 = (qp0k−2 − (k − 1)V pk−2 )/(k − 1)/Ak−1 , k = 2, 3, ..., ν.
If q ≡ 1, then Res(rk , z0 ) = pk−1 (z0 ), where pk−1 = p0k−2 /(k − 1).

In cases that concern us, both deg(p) P and ν increase with n. In effect,
suppose that n > m − n2 and z0 = 0. Let Res (rν , zi ) be as in (14). Then
deg(p) = 2(n + n2 ), ν = n2 + n + 1 − m, deg(q) = 2m and
deg(pk−1 ) = deg(p) + (k − 1)(deg(q) − 1), k = 1, · · · , ν.
If z0 6= 0, then ν is a constant that is usually not large. Note that |q(z0 )|
can be very small because the mapping z → (z + 1/z)/2 transforms every
difficult pole into two points that are very close to each other.
In general terms, let us suppose that deg(p) = O(n). We use algorithms
of order O(n) to carried out the calculation of p0k (z), pk (z0 ) and pk q, where q
is fixed. It follows that the previous procedure to compute {Res(rk , z0 )}νk=1
is of order O(n2 ), provided that ν = O(n). Finally, if we apply this algorithm
to (14)–(15), then the calculation of the coefficients of the n-point quadrature
requires O(n3 ) operations. Moreover, if we approximate λn,k (1/q) by the N th
partial sum of the series expansion in (15), then the computational effort is
N times greater than that we make to evaluate λn,k (h1 W1 /q), by using (14).
The range 64 ≤ N ≤ 512 is suggested by Example 5, though this should vary
according to the features of the integrand.
For the purposes of comparison, the reader is referred to [6, 1997] and
the bibliography therein. More recently, Deckers et al [3, 2009] presented
an algorithm of order O(n) to compute rational Gauss quadrature formulas,
provided that the set of poles is finite and does not depend on n. It is also
proved in [3] that the order is O(n2 ) when all poles are different.

11
4. Examples
From the previous corollaries we can derive several integration methods.
The following examples show that these procedures behave differently, de-
pending on how we select f and G.
Z δ√
1 + 0.5θx k −x
Example 1. Iδ (k, η, θ) = x e dx.
0 e + e−x
−η

Fermi-Dirac integrals are defined as F(k, η, θ) = limδ→∞ Iδ (k, η, θ). Fur-

ther, Iδ (k, η, θ) can be expressed as follows.
Z 1√
k+1 1 + 0.5θδx k −δx
Iδ (k, η, θ) = δ x e dx = (17)
0 e−η + e−δx
Z 1 p
k+1 1 + 0.5θδ(1 − x2 ) −δ(1−x2 ) (1 − x2 )k+1/2
δ e |x| √ dx.
−1 e−η + e−δ(1−x2 ) 1 − x2

Table 1: The integral of Example 1 is evaluated by rational Gauss quadrature (EG ),

and by an interpolatory rule (E) based on formula (12), with θ = 10−4 , η = −1, k = 0.5
and m = deg(q).

n m EG n m E
1 2 5.2e − 02 1 ∞ 9.6e − 04
2 4 1.3e − 03 2 ∞ 4.6e − 04
3 6 2.7e − 05 3 ∞ 1.0e − 07
4 8 2.2e − 07 4 ∞ 4.4e − 08
5 10 4.0e − 07 5 ∞ 1.7e − 11
7 14 1.9e − 12 7 ∞ 4.6e − 15

Concerning the integrals of Fermi-Dirac, the only values of k which are

interesting for physics are k = 1/2, 3/2, 5/2.√ In this case, the Chebyshev
weight function is W (x) = (1 − x2 )k+1/2 /(π 1 − x2 ).
For comparison purposes, the values we have adopted for η and θ, corre-
spond to examples previously studied by Gautschi [7]. Next we describe the
process of calculating F(k, η, θ), using routines based upon Corollary 1.
The very first step is that the integration interval must be truncated con-
veniently. Instead of F(1/2, −1, 10−4 ), we try to evaluate I40 (1/2, −1, 10−4 ),
which approximates the former up to 16 decimal figures. In turn, the new

12
 R1
integral on the interval [0, 40] can be transformed into −1 f (x)G(x)W (x)dx,
where
p E(x)
f (x) = π403/2 1 + 0.002(1 − x2 ), G(x) = , (18)
q∞ (x)
2 2
with E(x) = e−40(1−x ) |x|, q∞ (x) = e + e−40(1−x ) , and h1 (x) = (1 − x2 ).
The zeros of q∞ are much closer to [−1, 1] than those of the denominator
(e + e−x ), a fact which appears to be an added disadvantage.
Formula (12) allows us to obtain accurate results when calculating the
integral I40 (−1, 1/2, 10−4 ), up to the quadrature order n = 512. For purposes
of comparison we have tabulated some of them in Table 1 under the label
E. Column EG enumerates the relative errors reported by Gautschi [7] when
rational Gauss quadrature is applied to Fermi-Dirac integrals.
Z δ√
1 + 0.5θx k −x
Example 2. Jδ (k, η, θ) = x e dx.
0 e − e−x
−η

Bose-Einstein integrals are defined as B(k, η, θ) = limδ→∞ Jδ (k, η, θ).

Using a procedure similar to that applied to Fermi-Dirac integrals, we
have estimated the values of Jδ (−1, 0.5, 10−4 ). The corresponding errors are
shown in Table 2, to be compared with the results obtained by Gautschi [7].

Table 2: The integral of Example 2 is evaluated by rational Gauss quadrature (EG ),

and by an interpolatory rule (E) based on formula (12), with θ = 10−4 , η = −1, k = 0.5
and deg(q) = m.

n m EG n m E
1 1 2.5e − 01 1 ∞ 9.6e − 04
2 3 7.8e − 03 2 ∞ 4.6e − 04
3 5 1.7e − 04 3 ∞ 1.0e − 07
4 7 2.7e − 06 4 ∞ 4.6e − 08
5 9 6.5e − 07 5 ∞ 1.9e − 11
6 11 2.5e − 08 6 ∞ 7.3e − 12

1
x/ω √
Z
Example 3. 1 − x2 dx, ω > 1.
−1 sin(πx/ω)
If |ω − 1| is small and zn = ±nω, n = 1, 2, ..., then some points zn could
be classified as difficult poles.

13
 Table 3 depicts the relative errors when it is used a procedure based on
Corollary 2, with h1 (x) = (1 − x2 ) and q(x) = (1 − x2 /1.0012 ). Note that
when n is large, the numerical method is fairly stable despite the presence
of difficult poles. The entire process consumes less than a second when
n = 2048.
Table 3: Errors when the integral of Example 3 is evaluated by an interpolatory rule
based on Corollary 2, with ω = 1.001

n E n E
2 2.0e − 02 64 3.3e − 16
4 1.4e − 03 128 4.9e − 16
8 9.1e − 06 256 4.9e − 16
16 1.8e − 10 1024 9.8e − 16
32 1.6e − 16 2048 2.3e − 15

1
πx/ω √
Z
Example 4. 1 − x2 dx, |ω| > 1.
−1 sinh(πx/ω)
Column E1 of Table
√ 4 gives the relative errors whenQn (13)2is used to com-
pute λ2n+1,k (Gn (x)/ 1 − x ), with Gn (x) = (1 − x ) k=1 (x /(k ω ) + 1)−1 .
2 2 2 2

Column E2 lists the relative errors when all the poles have been moved to
the weight function, i.e. G(x) = (1 − x2 )(π 2 x/ω)/ sinh(πx/ω). The results
are compared with those obtained by using Gauss-Chebyshev quadrature for-
mula of rational type (cf. [2, Example 5.2]). In this example, all the poles of
the integrand are benign.

Table 4: Example 4 is evaluated by rational Gauss-Chebyshev quadrature (GCh), and

by an interpolatory rule (E) based on (13), with ω = 1.001.

n GCh E1 E2 n GCh E1 E2
3 1.1e − 03 1.3e − 02 1.4e − 14 17 3.8e − 16 3.7e − 15 1.3e − 14
5 4.5e − 07 4.7e − 05 1.4e − 14 33 1.9e − 16 4.0e − 15 1.3e − 14
9 5.8e − 16 7.1e − 11 1.3e − 14

Z 1
x/ω
Example 5. dx, ω > 1.
−1 sin(πx/ω)

14
 The difficult poles of Q
the integrand are exactly the same as those of
Example 3. Let q(x) = m 2 2 2
j=1 (j ω − x ), with ω > 1. In order to ap-
ply Corollary 3 we rewrite the integrand as √f (x)K(x)W1 (x)/q(x), where
f (x) = (πt/ω)q(x)/ sin(πt/ω) and K(x) = π 1 − x2 . Column EG of Ta-
ble 5 shows the errors when this integral is evaluated by a rational Gauss
quadrature according to Gautschi’s approach [7]. The numbers below EN are
the relative errors produced by the interpolatory rule (4),
√ when the coeffi-
cients are obtained from (15), and the Chebyshev series of 1 − x2 is replaced
by its N th partial sum. The errors corresponding to N = 64, 128, 256, are
displayed to show that is not worth to choose N > 128.

Table 5: Errors when the integral of Example 5 is evaluated by rational Gauss quadrature
(EG ) and by an interpolatory rule (EN ) based on Corollary 3, with ω = 1.001.

n EG E64 E128 E256

deg(q) = 2 2 8.4e − 02 1.8e − 01 1.9e − 01 1.8e − 01
3 1.3e − 03 2.7e − 02 2.7e − 02 2.7e − 02
4 1.7e − 04 1.1e − 02 1.1e − 02 1.1e − 02
5 2.0e − 05 1.5e − 03 1.5e − 03 1.5e − 03
6 9.3e − 06 6.4e − 04 6.6e − 04 6.6e − 04
7 7.8e − 06 5.8e − 05 8.1e − 05 8.1e − 05
10 2.0e − 05 1.8e − 06 1.9e − 06
12 2.2e − 05 7.1e − 08 9.2e − 08

deg(q) = 4 2 2.9e − 03 1.1e − 01 1.1e − 01 1.1e − 01

3 3.5e − 04 9.6e − 03 9.6e − 03 9.6e − 03
4 1.9e − 05 3.9e − 03 3.9e − 03 3.9e − 03
5 6.6e − 06 2.5e − 04 2.7e − 04 2.7e − 04
6 6.1e − 06 9.2e − 05 1.1e − 04 1.1e − 04
7 1.5e − 05 6.9e − 06 6.9e − 06
10 2.2e − 05 4.8e − 08 6.9e − 08
12 2.2e − 05 1.9e − 08 1.5e − 09

1
ex
Z
Example 6. dx.
−1 x2 + ε2
This integrand has difficult poles at ±εi, when |ε| is √
small. We rewrite
the integrand as f GW1 , where f (x) = e and G(x) = π 1 − x2 /(x2 + ε2 ).
x

Column E1 of Table 6 shows the accuracy obtained when this integral is

15
evaluated by formula (4), and its coefficients are estimated
√ using equation
(12) of Corollary 1, with K(x) = G(x). If K(x) = π 1 − x2 , then we
apply Corollary 3 to Example 6, with N = 512, to obtain the errors which
are shown under E3 . To compare with Monegato’s method, the only three
results reported in [9] appear in column EM . The results in Table 6 are
relatively poor for n ≥ 10, because Chebyshev nodes of the first kind seem
to be inadequate for integrals whose mass is highly concentrated near x = 0.
This also explains why the errors tend to be smaller when n is odd.
Table 6: Errors when the integral of Example 6 is evaluated by an interpolatory rule
based on Corollary 1 (E1 ) and Corollary 3 (E3 ), with ε = 0.01.

n E1 E3 EM n E1 E3 EM
1 3.2e − 03 3.2e − 01 7 3.7e − 09 3.7e − 09
2 2.6e − 01 2.6e − 01 8 1.8e − 07 1.8e − 07 1.3e − 07
3 1.1e − 04 1.1e − 04 9 9.3e − 12 3.1e − 11
4 5.2e − 03 5.2e − 03 3.5e − 03 10 5.0e − 10 4.8e − 10
5 8.3e − 07 8.3e − 07 11 6.5e − 13 2.1e − 11
6 4.2e − 05 4.2e − 05 16 6.3e − 13 2.1e − 11 2.0e − 16

In order to provide a variant of the previous method, we use the following

transformation of the integrand to move its poles to a region near the ends
of the interval.
Z 1 Z 1 √
ex cosh( 1 − x2 )|x|dx
2 2
dx = √ (19)
−1 x + ε
2 2 2
−1 (1 − x + ε ) 1 − x

Now, let us examine Example 6 using (19). The following table depicts
three different ways
√ of choosing the integrand f and the difficult factor G,
with W (x) = 1/(π 1 − x2 ) in all cases.
Case f √ G
I π cosh( 1 − x2 ) √ |x|/(1 − x2 + 0.012 )
2 2 2
II π|x| √ 1 − x )/(1 − x2 + 0.012 )
cosh(
III π cosh( 1 − x2 )|x|/(1 − x + 0.01 )
Table 7 lists the relative errors for the three cases, when the procedure is
based upon equation (13). The results for case II are affected by the presence
of |x| in the integrand. To increase accuracy, we have calculated 5 × 107
coefficients Aj . Case III of Table 7 indicates that, for n ≤ 15, the procedure
based on (19) can provide higher precision than Monegato’s method.

16
 Table 7: The integral of Example 6 is evaluated using (19), with ε = 0.01.

n I II III n I II III
2 2.6e − 01 2.9e − 01 6.0e − 15 64 6.9e − 15 6.0e − 06 8.0e − 15
4 5.2e − 03 3.5e − 02 6.0e − 15 128 1.2e − 14 5.8e − 07 1.0e − 14
8 1.8e − 07 4.4e − 03 6.4e − 15 256 1.0e − 14 7.3e − 08 1.0e − 14
16 6.9e − 15 5.3e − 04 7.1e − 15 512 1.4e − 14 1.6e − 08 1.2e − 14
32 5.1e − 15 5.9e − 05 5.3e − 15 1024 1.2e − 14 3.9e − 09 1.0e − 14

5. Note on the coefficients Aj

Both formulas (12) and (13) have shown to be a good reference to develop
routines in any programming language. They depend on P the first d + n + 1
coefficients of the Chebyshev series expansion K(x) = ∞ j=0 Aj Tj (x), which
in turn could depend on other parameters. For example, the previous schema
2
for Fermi-Dirac integrals refers to q∞ (x) = e−η +e−40(1−x ) , hence Aj = Aj (η).
A physical fact is that the parameter η depends on temperature, so the
coefficients Aj are often calculated and the procedure becomes more costly.
In such case it seems to be more convenient to put q∞ ≡ 1, namely
p
1 + 0.002(1 − x2 ) 2
f (x) = π403/2 −η −40(1−x 2) , G(x) = e−40(1−x ) |x|. (20)
e +e
Table 8 enumerates relative errors for different values of η when both the
integrand and the weight function are selected according to (20). Despite
of the convergence being slower, it is quite easy to reach a precision of 14
decimal figures. The advantage of this variant is that the coefficients are
calculated only once.

Table 8: Errors when the integral of Example 1 is evaluated for several values of η, with
θ = 10−4 , k = 0.5, q∞ ≡ 1.

η n=4 n = 16 n = 32 n = 64 n = 80 n = 128
−10−2 3.0e − 01 1.1e − 02 8.5e − 05 1.3e − 10 1.8e − 13 5.0e − 16
−1 1.2e − 01 1.0e − 02 2.2e − 05 1.4e − 12 1.3e − 14 9.6e − 16
−5 2.4e − 03 3.0e − 04 1.5e − 07 2.0e − 15 1.6e − 14 7.3e − 16
−10 1.6e − 05 2.0e − 06 1.3e − 09 6.9e − 15 2.1e − 14 4.7e − 15

17
 The coefficient Aj is the integral of K with respect to the weight 2Tj W1 ,
so we can construct a low cost procedure from the following formula.
n
X
Aj ≈ λn,k (2Tj W1 )K(xn,k ), j = 1, 2, ..., (21)
k=1
Pn
and A0 ≈ k=1 λn,k (W1 )K(xn,k ).
Let Pn−1,K be the Lagrange polynomial of degree n−1,
Pn−1which interpolates
K at the Chebyshev points {xn,k }. Then Pn−1,K (x) = j=0 dn,j Tj (x), where

n−1  !
2 X iπj(2k + 1)
dn,j = < K(xn,k ) exp , j ≥ 1, (22)
n k=0
2n

and dn,0 = n−1

P
k=0 K(xn,k )/n.
We use orthogonality and Proposition 2 to prove that (21) is exact in the
limit, i.e. limn dn,j = Aj .
From the above we deduce the following formulas for Aj .
(1)
If K is even, then A2j ≈ 2<(Aj,n ), where
n−1
e−jπi/n X
 
(1) 2kjπi
Aj,n = K(xn,k ) exp − , (23)
n k=0 n

and A0 must be halved.

(2)
If K is odd, then A2j+1 ≈ 2<(Aj,n ), where
n−1
e−jπi/n X
 
(2) −iπ(2k+1)/(2n)

2kjπi
Aj,n = K(xn,k )e exp − . (24)
n k=0 n

Both expressions (23) and (24) can be computed efficiently by means of any
F F T algorithm. The larger the value of n, the better the approximation
of the coefficients Aj . Nevertheless, the convergence of both formulas (23)
and (24), is not uniform with respect to j, as n → ∞. The size of n in
formulas (23) and (24) depends on the code we use to simulate K, as well as
the mathematical complexity of the latter.
In the experimental stage of this work, we have used Matlabr tools to
carry on the calculations with a precision of 15 decimal figures. The programs
have been run using a 64-bit CPU, with frequency of 3.3 GHz.

18
6. Concluding remarks
Formulation (4)–(5) and the theory in Section 2 attempt to go beyond the
well known rational approach, although the discretization process is mainly
established in terms of rational functions, as shown in Section 3. The non-
rational focus we want to introduce is given by Corollary 1. Using the latter
we have shown that, in some cases, a modified interpolatory rule can pro-
vide superior accuracy to that obtained by rational Gauss formula. Besides,
Corollary 1 allows to remove one of the worst drawbacks of the rational meth-
ods, namely, the need of calculating all the difficult poles of the integrand.
Likewise, the calculation of the coefficients can be done in a way much simpler
for some special cases.
An objection to the results of Section 3 is that some of them contemplate
the calculation of residues of rational functions at some poles, which is an
ill-conditioned problem when the denominator is near a polynomial with
multiple roots. That is why we have suggested a very simple algorithm in
Section 4. On the other hand, equation (12) of Corollary 1 is established
solely in terms of residues at the point z = 0, that are fairly easy to calculate
in terms of the derivatives of a polynomial. In principle, all cases can be
reduced to formula (13), which requires little coding effort.
We have implemented our approach for a type of integral which includes
the four kinds of Chebyshev weight functions. Any definite integral on a
bounded interval can be transformed into another integral in which appears
one of the four Chebyshev weight functions. The following equation differs
from (19) and corroborates the above statement.
Z b √
b − a 1 Z((b − a) 1 − x2 + a)|x|dx
Z
Z(x)dx = √ . (25)
a 2 −1 1 − x2

However, no formula is appropriate for all cases, and (25) is no exception.

Acknowledgments
The authors gratefully acknowledge the thorough work of the referees,
whose remarks have improved significantly the presentation of this article.

19
References
[1] E. Berriochoa, A. Cachafeiro, J.M. Garcı́a-Amor, Characterizing the
measures on the unit circle with exact quadrature formulas in the space
of polynomials, Comput. Math. Appl. 58, 2009, pp. 1370-1382.

[2] K. Deckers, J. Van Deun and A. Bultheel, Rational Gauss-Chebyshev

quadrature formulas for complex poles outside [−1, 1], Math. Comp.
77(262), 2008, pp. 967–983.

[3] K. Deckers, J. Van Deun and A. Bultheel, Computing rational Gauss-

Chebyshev quadrature formulas with complex poles: the algorithm, Ad-
vances in Engineering Software, 40(8), 2009, pp. 707–717.

[4] C. W. Clenshaw and A. R. Curtis, A method for numerical integration

on an automatic computer, Numer. Math. 2, 1960, pp. 197–205.

[5] J. R. Illán González, Interpolatory quadrature formulas for meromorphic

integrands, Jaen J. Approx. 1(2), 2009, pp. 175–193.

[6] W. Gautschi, Moments in Quadrature Problems, Computers Math. Ap-

plic. Vol. 33, No. 1/2, 1997, pp. 105–118.

[7] W. Gautschi, Algorithm 793: GQRAT-Gauss quadrature for rational

functions, ACM Trans. Math. Software, 25, 1999, pp. 213–239.

[8] J. C. Mason and D. C. Handscomb, Chebyshev polynomials, Chapman &

Hall/CRC, 2003.

[9] G. Monegato, Quadrature formulas for functions with poles near the in-
terval of integration, Math. Comp. 47 (1986), pp. 301–312.

[10] I. H. Sloan and W. E. Smith, Properties of interpolatory product inte-

gration rules, SIAM J. Numer. Anal. Vol 19, No 2, 1982, pp. 427–442.

20