# On Fast and Stable Implementation of Clenshaw-Curtis and Fejér-Type Quadrature Rules

Authors: Shuhuang Xiang; Guo He; Haiyong Wang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Based upon the fast computation of the coefficients of the interpolation polynomials at Chebyshev-type points by FFT, together with the efficient evaluation of the modified moments by forward recursions or by Oliver’s algorithms, this paper presents fast and stable interpolating integration algorithms, by using the coefficients and modified moments, for Clenshaw-Curtis, Fejér’s first- and second-type rules for Jacobi weights or Jacobi weights multiplied by a logarithmic function. Numerical examples illustrate the stability, efficiency, and accuracy of these quadratures.


---

## Body

## 1. Introduction

The interpolation quadrature of the Clenshaw-Curtis rule as well as Fejér-type formulas for(1)I[f]=∫-11f(x)w(x)dx≈IN[f]=∑k=0Nwkf(xk)
has been extensively studied since Fejér [1, 2] in 1933 and Clenshaw and Curtis [3] in 1960, where the nodes {xk} are of Chebyshev types while the weights {wk} are computed by sums of trigonometric functions. When xk=cos⁡((2k+1)π/(2N+2)) (k=0,1,…,N), this quadrature is calledFejér’s first-type rule. This kind of points is called the first kind of Chebyshev points, whileFejér’s second-type rule is corresponding to the Filippi points xk=cos⁡((k+1)π/(N+2)) (k=0,1,…,N) and theClenshaw-Curtis-type quadrature to the Clenshaw-Curtis points (the second kind of Chebyshev points) xk=cos⁡(kπ/N) (k=0,1,…,N). For more details, see Davis and Robinowitz [4], Sloan and Smith [5, 6], Sommariva [7], Trefethen [8], Waldvogel [9], and so forth.In the casew(x)≡1, a connection between the Fejér, Clenshaw-Curtis quadrature rules, and discrete Fourier transforms (DFTs) was given by Gentleman [10, 11], where the Clenshaw-Curtis rule is implemented with (N+1) nodes by means of a discrete cosine transformation (DCT). An independent approach along the same lines, unified algorithms based on DFTs of order N for generating the weights of the two Fejér rules and of the Clenshaw-Curtis rule, was presented in Waldvogel [9]. A streamlinedmatlab code is given as well in [9]. In addition, Clenshaw and Curtis [3], O’Hara and Smith [12], Trefethen [8, 13], Xiang and Bornemann in [14], Xiang et al. [15–17], and so forth showed that the Gauss, Clenshaw-Curtis, and Fejér quadrature rules are about equally accurate.In this paper, we focus the attention on the weight functionsw(x)=(1-x)α(1+x)β and w(x)=ln⁡⁡((1+x)/2)(1-x)α(1+x)β. For these two weight functions, the Clenshaw-Curtis-type quadrature has been extensively studied in a series of papers of Piessens [18, 19] and Piessens and Branders [20–23], by using Chebyshev interpolant QN[f](x)=∑n=0N‍anTn(x) of f(x) at the (N+1) Clenshaw-Curtis points together with the modified moments Mn=∫-11‍w(x)Tn(x)dx [24]:
(2)I[f]=∫-11f(x)w(x)dx≈INC-C[f]=∑n=0Nan∫-11Tn(x)w(x)dx=∑k=0NanMn‍,
where Tn(x) is the Chebyshev polynomial of degree n and an can be efficiently computed by FFT [8, 10, 11] which is widely used for the approximation of highly oscillatory integrals such as [22, 25–30]. The modified moments Mn=∫-11‍w(x)Tn(x)dx satisfy the following recurrence formulas for Jacobi weights or Jacobi weights multiplied by ln⁡⁡((x+1)/2) [20].(i)
Forw(x)=(1-x)α(1+x)β, by using Fasenmyer’s technique, the recurrence formula for the evaluation of the modified moments
(3)Qn(α,β)=∫-11(1-x)α(1+x)βTn(x)dxis(4)(β+α+k+2)Qn+1(α,β)+2(α-β)Qn(α,β)+(β+α-n+2)Qn-1(α,β)=0with(5)Q0(α,β)=2β+α+1Γ(α+1)Γ(β+1)Γ(β+α+2),Q1(α,β)=2β+α+1Γ(α+1)Γ(β+1)Γ(β+α+2)β-αβ+α+2.The forward recursion is numerically stable [20], except in two cases:
(6)α>β,β=-12,12,32,…(7)β>α,α=-12,12,32,….(ii)
Forw(x)=ln⁡((x+1)/2)(1-x)α(1+x)β, for
(8)Gn(α,β)=∫-11ln⁡⁡(x+12)(1-x)α(1+x)βTn(x)dx,the recurrence formula [20] is
(9)(β+α+k+2)Gn+1(α,β)+2(α-β)Gn(α,β)+(β+α-n+2)Gn-1(α,β)=2Qn(α,β)-Qn-1(α,β)-Qn+1(α,β)with(10)G0(α,β)=-2β+α+1Φ(α,β+1),G1(α,β)=-2β+α+1[2Φ(α,β+2)-Φ(α,β+1)],where(11)Φ(α,β)=B(α+1,β)[Ψ(α+β+1)-Ψ(β)].B(x,y) is the Beta function and Ψ(x) is the Psi function (see Abramowitz and Stegun [31]). The forward recursion is as numerically stable as (4) except for (6) or (7) [20].Thus, the modified moments can be fast computed by the forward recursions (4) and (9) except the two cases (6) and (7), and the total costs for IN[f] are O(Nlog⁡N) operations.However, in case (6) or (7), the accuracy of the forward recursion is catastrophic particularly when α-β≫1 or β-α≫1 and n≫1 (see Table 1). In case (6) the relative errors ϵn of the computed values Qn(α,β) obtained by the forward recursion behave approximately as
(12)ϵn~n2(α-β),n⟶∞
and in case (7) as
(13)ϵn~n2(β-α),n⟶∞.Table 1
Computation ofMn(α,β)=∫-11‍(1-x)α(1+x)βTn(x)dx with different n and (α,β) by the forward recursion of (4).
n
5
10
100
Exact value for (20,-0.5)
-
1.734810854604316
e
+
05
4.049003666168904
e
+
03
-
3.083991348593134
e
-
41
(4) for (20,-0.5)
-
1.734810854604308
e
+
05
4.049003666169083
e
+
03
1.787242305340324
e
-
11
Exact value for (100,-0.5)
-
2.471295049468578
e
+
29
1.174275526131223
e
+
29
2.805165440968788
e
-
29
(4) for (100,-0.5)
-
2.471295049468764
e
+
29
1.174275526131312
e
+
29
-
1.380038973213404
e
+
13In this paper, we will consider interpolation approaches for Clenshaw-Curtis rules as well as Fejér’s first- and second-type formulas for(14)I[f]=∫-11f(x)w(x)dx≈∑k=0Nan∫-11Tn(x)w(x)dx‍=∑k=0NanMn‍≔IN[f]
with w(x)=(1-x)α(1+x)β or w(x)=(1-x)α(1+x)βln⁡⁡((1+x)/2), which can be efficiently calculated in O(Nlog⁡N) operations. Computing the modified moments Mn in cases (6) and (7) by Oliver’s algorithm [32] or Lozier’s algorithm [33] with one starting value and one end value, as well as offering the very short codes for the evaluation of the coefficients by FFT, is the topic of this paper.This paper is organized as follows. In Section2.1, we studied the asymptotic expansions of the modified moments Mn. Based on the results of the asymptotic expansions, we discussed Oliver’s algorithm for the modified moments Mn when the recursion is unstable in Section 2.2. Moreover, we gave some test to verify the accuracy and the CPU time of Oliver’s algorithm. In Section 2.3, the concisematlab codes for evaluation of the coefficients an by FFT are presented. The efficiency and accuracy of the three quadratures are illustrated in Section 3.


## 2. Computation of the Modified Moments and the Coefficients of the Interpolation Polynomials

From (4) and (9), we see that if Qn(α,β) can be efficiently and stably computed then Gn(α,β) can be so. Note that the asymptotic behaviour of two linearly independent solutions of (4) satisfies
(15)yn,1~n-2α-2[1+O(n-2)],yn,2~n-2β-2[1+O(n-2)],n-2βi[1+O(n-2)]n⟶∞
(see Denef and Piessens [34] and Piessens and Branders [20]). Consequently, the forward recursions for (4) and (9) are perfectly stable except cases (6) and (7). In these two cases, both the forward recursions and backward recursions for (4) and (9) are numerically unstable. In the following, we will consider Oliver’s algorithms to evaluate modified moments based on their asymptotic formulas of the moments with one starting value and one end value for these two cases.


### 2.1. Asymptotic Expansions of the Modified Moments

Lemma 1 (Erdélyi [35]).
If0<λ, μ≤1, and ϕ(t) is m times continuously differentiable for α≤t≤β, then
(16)∫αβeixt(t-α)λ-1(β-t)μ-1ϕ(t)dt‍=Bm(x)-Am(x)+O(x-m),
where
(17)Am(x)=∑n=0m-1Γ(n+λ)n!eiπ(n+λ-2)/2x-n-λ×eixα[dndtn{(β-t)μ-1ϕ(t)}]t=α,Bm(x)=∑n=0m-1Γ(n+μ)n!eiπ(n-μ)/2x-n-μ×eixβ[dndtn{(t-α)λ-1ϕ(t)}]t=β.Lemma 2.
If0<λ, 0<μ, and ϕ(t) is m times continuously differentiable for α≤t≤β, then
(18)∫αβeixt(t-α)λ-1(β-t)μ-1ln⁡⁡(β-t)ϕ(t)dt=Bm(x)-Am(x)+O(x-m),
where
(19)Am(x)=∑n=0m-1Γ(n+λ)n!eiπ(n+λ-2)/2x-n-λ×eixα[dndtn{(β-t)μ-1ln⁡⁡(β-t)ϕ(t)}]t=α,Bm(x)=∑n=0m-1Γ(n+μ)n!eiπ(n-μ)/2×x-n-μ{Ψ(n+μ)-ln⁡⁡(x)-iπ2}×eixβ[dndtn{(t-α)λ-1ϕ(t)}]t=β.Proof.
The proof can be directly derived from the proof of Lemma1 and the proof of the Theorem 5 in Erdélyi [36].Theorem 3.
If-1<α,-1<β, then
(20)Qn(α,β)=2β-α∑k=0m-1ak(α,β)h(α+k)+(-1)n2α-β∑k=0m-1ak(β,α)h(β+k)+O(n-2min⁡⁡(α,β)-2m),
where
(21)h(α)=cos⁡[π(α+1)]Γ(2α+2)n-2α-2,a0(α,β)=1,a1(α,β)=-α12-β4-16,a2(α,β)=1120+19α1440+α2288+αβ48+β32+β232,a3(α,β)=-15040-β960-107α181440-β2384-α21920-β3384-α310368-7αβ2880-α2β1152-αβ2384.Proof.
For-1<α,β≤-1/2, taking x=cos⁡(θ) in (3), we have
(22)Qn(α,β)=∫0π(1-cos⁡⁡(θ))α+1/2×(1+cos⁡⁡(θ))β+1/2cos⁡⁡(nθ)dθ=∫0πφ(θ)θ2α+1(π-θ)2β+1cos⁡⁡(nθ)dθ,
where
(23)φ(θ)=(1-cos⁡⁡(θ)θ2)α+1/2(1+cos⁡⁡(θ)(π-θ)2)β+1/2.
Consequently, the desired result can be derived by applying (16) to (22).
The above outcome can be extended to the case of-1/2<α, -1/2<β, since Qn(α,β) can be written as
(24)Qn(α,β)=∫0π[φ(θ)θ⌊2α+1⌋+1(π-θ)⌊2β+1⌋+1]×θ2α+1-⌊2α+1⌋-1(π-θ)2β+1-⌊2β+1⌋-1×cos⁡⁡(nθ)dθ,
where ⌊z⌋ denotes the largest integer less than or equal to z.Theorem 4.
If-1<α,-1<β, then
(25)Gn(α,β)=2β-α∑k=0m-1ckh(α+k)+(-1)n2α-β‍×∑k=0m-1h(β+k)(2ak(β,α)ϕ(β+k)+bk)+O(n-2m),
where
(26)ϕ(β)=Ψ(2β+2)-ln⁡⁡(2n)-π2tan⁡(πβ),b0=0,b1=-112,b2=191440+α48+β144,b3=-7α2880-β960-α2384-β23456-107181440-αβ576,c0=0,c1=-14,c2=132+α48+β16,c3=-7α2880-β192-α21152-β2128-1960-αβ192.Proof.
Lettingx=cos⁡(θ) in (8) and using ln⁡⁡((1+cos⁡(θ))/2)=ln⁡⁡((1+cos⁡(θ))/2(π-θ)2)+2ln⁡⁡(π-θ), we have
(27)Gn(α,β)=∫0πln⁡⁡(1+cos⁡⁡(θ)2(π-θ)2)×φ(θ)θ2α+1(π-θ)2β+1cos⁡⁡(nθ)dθ+2∫0πφ(θ)ln⁡⁡(π-θ)θ2α+1‍×(π-θ)2β+1cos⁡⁡(nθ)dθ.
Applying Lemmas 1 and 2 to the two integrals on the right hand side of (27), respectively, leads to the desired result.Remark 5.
Piessens and Branders gave the first term of asymptotic expansion forQn(α,β) and Gn(α,β) in [20] using the asymptotic theory of Fourier coefficients (Lighthill [37]). Furthermore, Piessens [18] presented the asymptotic expansion with explicit formulas for the first three terms for Qn(α,β), which is equivalent to (20) except a minor clerical error. In [18], a1(α,β)=-((α+2β+2)/3)2β-α-2 should be modified by a1(α,β)=-((α+3β+2)/3)2β-α-2. The first term of (25) is different from the first asymptotic term for Gn(α,β) proposed in [18]; that is, 2β-α∑k=01‍ckh(α+k)+(-1)n2α-β  ∑k=00‍h(β+k)  [2ak(β,α)ϕ(β+k)+bk]. In order to ensure the accuracy of the modified moments evaluated by Oliver’s method, the end value of modified moments must be computed very precisely. In this paper, we use the first four-term truncation.


### 2.2. Oliver’s Algorithm

Let(28)AN∶=(2(α-β)α+β+2+1α+β+2-22(α-β)α+β+2+2⋱⋱⋱α+β+2-(N-1)2(α-β)α+β+2+(N-1)α+β+2-N2(α-β)),(29)bN∶=(-(β+α+1)Q0(α,β)0⋯0-(α+β+2+N)QN+1)T,where “·T” denotes the transpose; then the modified moments can be solved by
(30)ANM=bN,M=(Q1(α,β),Q2(α,β),…,QN(α,β))T,
where QN+1(α,β) is computed by hypergeometric function [20] when N≤2000:
(31)QN+1(α,β)=2α+β+1Γ(α+1)Γ(β+1)Γ(α+β+2)×3F2([N+1,-N-1,α+1],[12,α+β+2],1).
If N>2000, QN+1(α,β) is computed by using the asymptotic expression (20) with the first four-term truncation. Particularly, Oliver’s algorithm can be fast implemented by applying LU factorization (chasing method) with O(N) operations.In addition, for the weightw(x)=ln⁡⁡((x+1)/2)(1-x)α(1+x)β in case of (6) or (7), the computation of the modified moments can be fixed up by Oliver’s algorithm similar to (29) with one starting G0(α,β)=-2α+β+1Φ(α,β+1) and one end GN+1(α,β). Consequently, the modified moments can be solved by
(32)ANG=CN,G=(G1(α,β),G2(α,β),…,GN(α,β))T,
where
(33)CN=(c1c2⋯cN-1cN)T,c1=-(β+α+1)G0(α,β)+2Q1(α,β)-Q0(α,β)-Q2(α,β),ck=2Qk(α,β)-Qk-1(α,β)-Qk+1(α,β),k=2,…,N-1,cN=-(β+α+N+2)GN+1(α,β)+2QN(α,β)-QN-1(α,β)-QN+1(α,β).
The end value GN+1(α,β) can be calculated by its asymptotic formula (25) with the first four-term truncation.Tables2 and 3 show the accuracy of Oliver’s algorithm for α=100 and β=-0.5, α=-0.5 and β=100, respectively. Table 4 displays the accuracy of Oliver’s algorithm for α=0.6 and β=-0.5, α=10 and β=-0.5 for Qn(α,β). Table 5 displays the accuracy of Oliver’s algorithm for α=-0.4999 and β=-0.5, α=0.9999 and β=-0.5 for Gn(α,β). Table 6 shows the CPU time for implementation of the two Oliver’s algorithms. From these tables, we see that Oliver’s algorithms have high efficiency and are precise. Thematlab codes on Oliver’s algorithms and all thematlab codes in this paper can be downloaded from [38].Table 2
Computation ofQn(α,β)=∫-11‍(1-x)α(1+x)βTn(x)dx with (α,β)=(100,-0.5) and different n by Oliver's algorithm.
n
100
500
1000
Exact value for (100,-0.5)
2.805165440968788
e
-
29
-
2.283851909785347
e
-
198
-
1.247890461118514
e
-
259
Oliver's method for (100,-0.5)
2.805165440968861
e
-
29
-
2.283851909785405
e
-
198
-
1.247890461118544
e
-
259Table 3
Computation ofGn(α,β)=∫-11‍(1-x)α(1+x)βln⁡((1+x)/2)Tn(x)dx with (α,β)=(-0.5,100) and different n by Oliver's algorithm compared with that computed by the forward recursion (9).
n
100
500
1000
Exact value for (-0.5,100)
1.089944378602585
e
-
28
7.222157005510106
e
-
198
5.715301877322031
e
-
259
Oliver's method for (-0.5,100)
1.089944378602206
e
-
28
7.222157005510654
e
-
198
5.715301877322483
e
-
259
(9) for (-0.5,100)
-
5.331299059334499
e
+
14
-
1.061058894110758
e
+
14
-
5.304494050667818
e
+
13Table 4
Computation ofQn(α,β)=∫-11‍(1-x)α(1+x)βTn(x)dx with different n and (α,β) by Oliver's algorithm.
n
2000
4000
8000
Exact value for (0.6,-0.5)
9.551684021848334
e
-
12
1.039402748103725
e
-
12
1.131065744497495
e
-
13
Oliver's method for (0.6,-0.5)
9.551684021848822
e
-
12
1.039402748103918
e
-
12
1.131065744497332
e
-
13
Exact value for (10,-0.5)
-
8.412345942129556
e
-
57
-
2.005493070382270
e
-
63
-
4.781368848995069
e
-
70
Oliver's method for (10,-0.5)
-
8.412345942129623
e
-
57
-
2.005493070382302
e
-
63
-
4.781368848995179
e
-
70Table 5
Computation ofGn(α,β)=∫-11‍(1-x)α(1+x)βln⁡((1+x)/2)Tn(x)dx with different n and (α,β) by Oliver's algorithm.
n
10
100
500
Exact value for (-0.4999,-0.5)
-
0.314181354550401
-
0.031418104511487
-
0.006283620842004
Oliver's method for (-0.4999,-0.5)
-
0.314181354550428
-
0.031418104511490
-
0.006283620842004
Exact value for (0.9999,-0.5)
-
0.895286620533541
-
0.088858164406923
-
0.017770353274330
Oliver's method for (0.9999,-0.5)
-
0.895286620533558
-
0.088858164406925
-
0.017770353274330Table 6
The CPU time for calculation of the first(N+1) modified moments by the Oliver's methods for α=-0.5 and β=100.
Modified moments
N
=
1
0
3
N
=
1
0
4
N
=
1
0
5
N
=
1
0
6
{
Q
n
(
α
,
β
)
}
n
=
0
N
0.003954 s
0.011553 s
0.121553 s
1.200070 s
{
G
n
(
α
,
β
)
}
n
=
0
N
0.006087 s
0.028416 s
0.296032 s
2.933928 s


### 2.3. Fast Computation of the Coefficients by FFT

The coefficientsan for the interpolation polynomial QN[f](x)=∑n=0N‍anTn(x) at the three kinds of points can be efficiently computed by FFT. For the Clenshaw-Curtis points, the elegantmatlab code on the an is from [8] (see Algorithm 1).Algorithm 1:
f
unction  a  =  clenshaw_curtis(f,N)                    %  (N+1)-coefficients for  C-C quadrature
x  =  cos(pi*(0:N)’/N);                   %  C-C  points
fx  =  feval(f,x)/(2*N);                       %  f  evaluated at these points
g  =  fft(fx([1:N+1  N:−1:2]));                 %  FFT
a  =  [g(1);  g(2:N)+g(2*N:−1:N+2);  g(N+1)];            %  Chebyshev coefficientsFor the other two classes of points, by the sums of trigonometric functions at these two point-set, it is not difficult to get the FFT implementation. Here, we will not give details but just offer the followingmatlab functions.For the first kind of Chebyshev points, we presented amatlab code for computing the coefficients an by FFT in Algorithm 2.Algorithm 2:
f
unction a = fejer1
(f,N)
x =
c
os(pi
*
(2
*
(0:N)’+1)/(2*N+2));                % coefficients for Fejér’s first rule
f
x  = feval(f,x)/(N+1);                       % first kind of Chebyshev points
g = fft(fx([1:N+1  N+1:−1:1]));               % f  evaluated at these points
hx = real(exp(2
*
1i
*
pi
*
(0:2
*
n+1)/(4*n+4))
.
*
g’);   % FFT
a  = hx(1:n+1);a(1)=0.5*a(1);             % Chebyshev coefficientsFor the Filippi points, we presented a matlab code for computing the coefficientsan by FFT in Algorithm 3.Algorithm 3:
f
unction a = fejer2(f,N)                % coefficients forFejér’s  second rule
x = cos(pi
*
(0:N+2)’/(N+2));
fx = feval(f,x)/(2N+4);                % f evaluated at these points
g = fft(fx([1:N+3  N+2:−1:2]));            %  FFT
b = [g(1);  g(2:N+2)+g(2*N+4:−1:N+4);  g(N+3)];
b(N+1:
−
2:1) = b(N+1:-2:1)
−
2
*
b(N+3);
b(N:
−
2:1) = b(N:-2:1)
−
b(N+2);
b(1) = b(1)+mod(N+1,2)
*
b(N+3)+mod(N,2)
*
b(N+2)/2;
a = b(1:N+1);                     %  Chebyshev coefficientsIn addition, the coefficients of the interpolation polynomialsQN[f](x) at these Chebyshev-type points may be also efficiently evaluated by DCT and inverse discrete sine transform (IDST), respectively. The discrete cosine transform denoted by Y=dct(X) is closely related to the discrete Fourier transform but using purely real numbers and takes O(Nlog⁡N) operations for
(34)mmmY(k)=σ(k)∑s=1NX(s)cos⁡⁡((k-1)π(2s-1)2N)withσ(1)=1N,σ(k)=2N,for2≤k≤N.
The discrete sine transform denoted by Y=dst(X) and its inverse by X=idst(Y) both take O(Nlog⁡N) operations for
(35)Y(k)=∑s=1NX(s)sin⁡(kπsN+1).Notice that the coefficientsaj for the interpolation polynomial QN[f](x)=∑n=1N+1‍′an-1Tn-1(x) at {cos⁡((2k-1)π/(2N+2))}k=1N+1 are represented by
(36)an-1=2N+1∑s=1N+1‍f(cos⁡((2s-1)π2N+2))×cos⁡((2s-1)(n-1)π2N+2),n=1,2,…,N+1,
and a~n for the interpolation polynomial QN[f](x)=∑n=1N+1‍a~n-1Un-1(x) at {cos⁡(kπ/(N+2))}k=1N+1 satisfies
(37)f(cos⁡⁡(kπN+2))sin⁡(kπN+2)=∑s=1N+1a~s-1sin⁡(skπN+2),k=1,2,…,N+1.
Then both can be efficiently calculated by DCT and IDST, respectively. Thematlab codes on the DCT and IDST are very short and just need three rows. Now, we give the codes in Algorithm 4 and Algorithm 5.Algorithm 4:
f
unction a = fejer1dct(f,N)            % coefficients forFejér’s  second rule
x = cos(pi
*
(2
*
(0:N)’+1)/(2*N+2));         % The first kind of Chebyshev points
fx = feval(f,x);                   %f  evaluated at these points
a = dct(fx)
*
sqrt(2/(N+1));a(1)=a(1)/sqrt(2);    % Chebyshev coefficientsAlgorithm 5:
f
unction a = fejer2idst(f,N)            % coefficients forFejér’s  second rule
x  = cos(pi*(1:N+1)’/(N+2));             % Filippi points
fx = feval(f,x).
*
sin(pi
*
(1:N+1)’/(N+2));      % f evaluated at  these points
a = idst(fx);                     % Chebyshev coefficientsFor the first kind of Chebyshev points, amatlab code for computing the coefficients an by DCT is presented in Algorithm 4.As far as the Filippi points, amatlab code for computing the coefficients a~n by IDST is presented in Algorithm 5.Remark 6.
For the modified moments on the second kind Chebyshev polynomialM~n∶=∫-11‍w(x)Un(x)dx, we have
(38)M~0=M0,M~1=2M1,M~n=2Mn+M~n-2,n=2,3,…;
here we use the simple equation Un+2=2Tn+2+Un (see [31, pp. 778]).


## 2.1. Asymptotic Expansions of the Modified Moments

Lemma 1 (Erdélyi [35]).
If0<λ, μ≤1, and ϕ(t) is m times continuously differentiable for α≤t≤β, then
(16)∫αβeixt(t-α)λ-1(β-t)μ-1ϕ(t)dt‍=Bm(x)-Am(x)+O(x-m),
where
(17)Am(x)=∑n=0m-1Γ(n+λ)n!eiπ(n+λ-2)/2x-n-λ×eixα[dndtn{(β-t)μ-1ϕ(t)}]t=α,Bm(x)=∑n=0m-1Γ(n+μ)n!eiπ(n-μ)/2x-n-μ×eixβ[dndtn{(t-α)λ-1ϕ(t)}]t=β.Lemma 2.
If0<λ, 0<μ, and ϕ(t) is m times continuously differentiable for α≤t≤β, then
(18)∫αβeixt(t-α)λ-1(β-t)μ-1ln⁡⁡(β-t)ϕ(t)dt=Bm(x)-Am(x)+O(x-m),
where
(19)Am(x)=∑n=0m-1Γ(n+λ)n!eiπ(n+λ-2)/2x-n-λ×eixα[dndtn{(β-t)μ-1ln⁡⁡(β-t)ϕ(t)}]t=α,Bm(x)=∑n=0m-1Γ(n+μ)n!eiπ(n-μ)/2×x-n-μ{Ψ(n+μ)-ln⁡⁡(x)-iπ2}×eixβ[dndtn{(t-α)λ-1ϕ(t)}]t=β.Proof.
The proof can be directly derived from the proof of Lemma1 and the proof of the Theorem 5 in Erdélyi [36].Theorem 3.
If-1<α,-1<β, then
(20)Qn(α,β)=2β-α∑k=0m-1ak(α,β)h(α+k)+(-1)n2α-β∑k=0m-1ak(β,α)h(β+k)+O(n-2min⁡⁡(α,β)-2m),
where
(21)h(α)=cos⁡[π(α+1)]Γ(2α+2)n-2α-2,a0(α,β)=1,a1(α,β)=-α12-β4-16,a2(α,β)=1120+19α1440+α2288+αβ48+β32+β232,a3(α,β)=-15040-β960-107α181440-β2384-α21920-β3384-α310368-7αβ2880-α2β1152-αβ2384.Proof.
For-1<α,β≤-1/2, taking x=cos⁡(θ) in (3), we have
(22)Qn(α,β)=∫0π(1-cos⁡⁡(θ))α+1/2×(1+cos⁡⁡(θ))β+1/2cos⁡⁡(nθ)dθ=∫0πφ(θ)θ2α+1(π-θ)2β+1cos⁡⁡(nθ)dθ,
where
(23)φ(θ)=(1-cos⁡⁡(θ)θ2)α+1/2(1+cos⁡⁡(θ)(π-θ)2)β+1/2.
Consequently, the desired result can be derived by applying (16) to (22).
The above outcome can be extended to the case of-1/2<α, -1/2<β, since Qn(α,β) can be written as
(24)Qn(α,β)=∫0π[φ(θ)θ⌊2α+1⌋+1(π-θ)⌊2β+1⌋+1]×θ2α+1-⌊2α+1⌋-1(π-θ)2β+1-⌊2β+1⌋-1×cos⁡⁡(nθ)dθ,
where ⌊z⌋ denotes the largest integer less than or equal to z.Theorem 4.
If-1<α,-1<β, then
(25)Gn(α,β)=2β-α∑k=0m-1ckh(α+k)+(-1)n2α-β‍×∑k=0m-1h(β+k)(2ak(β,α)ϕ(β+k)+bk)+O(n-2m),
where
(26)ϕ(β)=Ψ(2β+2)-ln⁡⁡(2n)-π2tan⁡(πβ),b0=0,b1=-112,b2=191440+α48+β144,b3=-7α2880-β960-α2384-β23456-107181440-αβ576,c0=0,c1=-14,c2=132+α48+β16,c3=-7α2880-β192-α21152-β2128-1960-αβ192.Proof.
Lettingx=cos⁡(θ) in (8) and using ln⁡⁡((1+cos⁡(θ))/2)=ln⁡⁡((1+cos⁡(θ))/2(π-θ)2)+2ln⁡⁡(π-θ), we have
(27)Gn(α,β)=∫0πln⁡⁡(1+cos⁡⁡(θ)2(π-θ)2)×φ(θ)θ2α+1(π-θ)2β+1cos⁡⁡(nθ)dθ+2∫0πφ(θ)ln⁡⁡(π-θ)θ2α+1‍×(π-θ)2β+1cos⁡⁡(nθ)dθ.
Applying Lemmas 1 and 2 to the two integrals on the right hand side of (27), respectively, leads to the desired result.Remark 5.
Piessens and Branders gave the first term of asymptotic expansion forQn(α,β) and Gn(α,β) in [20] using the asymptotic theory of Fourier coefficients (Lighthill [37]). Furthermore, Piessens [18] presented the asymptotic expansion with explicit formulas for the first three terms for Qn(α,β), which is equivalent to (20) except a minor clerical error. In [18], a1(α,β)=-((α+2β+2)/3)2β-α-2 should be modified by a1(α,β)=-((α+3β+2)/3)2β-α-2. The first term of (25) is different from the first asymptotic term for Gn(α,β) proposed in [18]; that is, 2β-α∑k=01‍ckh(α+k)+(-1)n2α-β  ∑k=00‍h(β+k)  [2ak(β,α)ϕ(β+k)+bk]. In order to ensure the accuracy of the modified moments evaluated by Oliver’s method, the end value of modified moments must be computed very precisely. In this paper, we use the first four-term truncation.


## 2.2. Oliver’s Algorithm

Let(28)AN∶=(2(α-β)α+β+2+1α+β+2-22(α-β)α+β+2+2⋱⋱⋱α+β+2-(N-1)2(α-β)α+β+2+(N-1)α+β+2-N2(α-β)),(29)bN∶=(-(β+α+1)Q0(α,β)0⋯0-(α+β+2+N)QN+1)T,where “·T” denotes the transpose; then the modified moments can be solved by
(30)ANM=bN,M=(Q1(α,β),Q2(α,β),…,QN(α,β))T,
where QN+1(α,β) is computed by hypergeometric function [20] when N≤2000:
(31)QN+1(α,β)=2α+β+1Γ(α+1)Γ(β+1)Γ(α+β+2)×3F2([N+1,-N-1,α+1],[12,α+β+2],1).
If N>2000, QN+1(α,β) is computed by using the asymptotic expression (20) with the first four-term truncation. Particularly, Oliver’s algorithm can be fast implemented by applying LU factorization (chasing method) with O(N) operations.In addition, for the weightw(x)=ln⁡⁡((x+1)/2)(1-x)α(1+x)β in case of (6) or (7), the computation of the modified moments can be fixed up by Oliver’s algorithm similar to (29) with one starting G0(α,β)=-2α+β+1Φ(α,β+1) and one end GN+1(α,β). Consequently, the modified moments can be solved by
(32)ANG=CN,G=(G1(α,β),G2(α,β),…,GN(α,β))T,
where
(33)CN=(c1c2⋯cN-1cN)T,c1=-(β+α+1)G0(α,β)+2Q1(α,β)-Q0(α,β)-Q2(α,β),ck=2Qk(α,β)-Qk-1(α,β)-Qk+1(α,β),k=2,…,N-1,cN=-(β+α+N+2)GN+1(α,β)+2QN(α,β)-QN-1(α,β)-QN+1(α,β).
The end value GN+1(α,β) can be calculated by its asymptotic formula (25) with the first four-term truncation.Tables2 and 3 show the accuracy of Oliver’s algorithm for α=100 and β=-0.5, α=-0.5 and β=100, respectively. Table 4 displays the accuracy of Oliver’s algorithm for α=0.6 and β=-0.5, α=10 and β=-0.5 for Qn(α,β). Table 5 displays the accuracy of Oliver’s algorithm for α=-0.4999 and β=-0.5, α=0.9999 and β=-0.5 for Gn(α,β). Table 6 shows the CPU time for implementation of the two Oliver’s algorithms. From these tables, we see that Oliver’s algorithms have high efficiency and are precise. Thematlab codes on Oliver’s algorithms and all thematlab codes in this paper can be downloaded from [38].Table 2
Computation ofQn(α,β)=∫-11‍(1-x)α(1+x)βTn(x)dx with (α,β)=(100,-0.5) and different n by Oliver's algorithm.
n
100
500
1000
Exact value for (100,-0.5)
2.805165440968788
e
-
29
-
2.283851909785347
e
-
198
-
1.247890461118514
e
-
259
Oliver's method for (100,-0.5)
2.805165440968861
e
-
29
-
2.283851909785405
e
-
198
-
1.247890461118544
e
-
259Table 3
Computation ofGn(α,β)=∫-11‍(1-x)α(1+x)βln⁡((1+x)/2)Tn(x)dx with (α,β)=(-0.5,100) and different n by Oliver's algorithm compared with that computed by the forward recursion (9).
n
100
500
1000
Exact value for (-0.5,100)
1.089944378602585
e
-
28
7.222157005510106
e
-
198
5.715301877322031
e
-
259
Oliver's method for (-0.5,100)
1.089944378602206
e
-
28
7.222157005510654
e
-
198
5.715301877322483
e
-
259
(9) for (-0.5,100)
-
5.331299059334499
e
+
14
-
1.061058894110758
e
+
14
-
5.304494050667818
e
+
13Table 4
Computation ofQn(α,β)=∫-11‍(1-x)α(1+x)βTn(x)dx with different n and (α,β) by Oliver's algorithm.
n
2000
4000
8000
Exact value for (0.6,-0.5)
9.551684021848334
e
-
12
1.039402748103725
e
-
12
1.131065744497495
e
-
13
Oliver's method for (0.6,-0.5)
9.551684021848822
e
-
12
1.039402748103918
e
-
12
1.131065744497332
e
-
13
Exact value for (10,-0.5)
-
8.412345942129556
e
-
57
-
2.005493070382270
e
-
63
-
4.781368848995069
e
-
70
Oliver's method for (10,-0.5)
-
8.412345942129623
e
-
57
-
2.005493070382302
e
-
63
-
4.781368848995179
e
-
70Table 5
Computation ofGn(α,β)=∫-11‍(1-x)α(1+x)βln⁡((1+x)/2)Tn(x)dx with different n and (α,β) by Oliver's algorithm.
n
10
100
500
Exact value for (-0.4999,-0.5)
-
0.314181354550401
-
0.031418104511487
-
0.006283620842004
Oliver's method for (-0.4999,-0.5)
-
0.314181354550428
-
0.031418104511490
-
0.006283620842004
Exact value for (0.9999,-0.5)
-
0.895286620533541
-
0.088858164406923
-
0.017770353274330
Oliver's method for (0.9999,-0.5)
-
0.895286620533558
-
0.088858164406925
-
0.017770353274330Table 6
The CPU time for calculation of the first(N+1) modified moments by the Oliver's methods for α=-0.5 and β=100.
Modified moments
N
=
1
0
3
N
=
1
0
4
N
=
1
0
5
N
=
1
0
6
{
Q
n
(
α
,
β
)
}
n
=
0
N
0.003954 s
0.011553 s
0.121553 s
1.200070 s
{
G
n
(
α
,
β
)
}
n
=
0
N
0.006087 s
0.028416 s
0.296032 s
2.933928 s


## 2.3. Fast Computation of the Coefficients by FFT

The coefficientsan for the interpolation polynomial QN[f](x)=∑n=0N‍anTn(x) at the three kinds of points can be efficiently computed by FFT. For the Clenshaw-Curtis points, the elegantmatlab code on the an is from [8] (see Algorithm 1).Algorithm 1:
f
unction  a  =  clenshaw_curtis(f,N)                    %  (N+1)-coefficients for  C-C quadrature
x  =  cos(pi*(0:N)’/N);                   %  C-C  points
fx  =  feval(f,x)/(2*N);                       %  f  evaluated at these points
g  =  fft(fx([1:N+1  N:−1:2]));                 %  FFT
a  =  [g(1);  g(2:N)+g(2*N:−1:N+2);  g(N+1)];            %  Chebyshev coefficientsFor the other two classes of points, by the sums of trigonometric functions at these two point-set, it is not difficult to get the FFT implementation. Here, we will not give details but just offer the followingmatlab functions.For the first kind of Chebyshev points, we presented amatlab code for computing the coefficients an by FFT in Algorithm 2.Algorithm 2:
f
unction a = fejer1
(f,N)
x =
c
os(pi
*
(2
*
(0:N)’+1)/(2*N+2));                % coefficients for Fejér’s first rule
f
x  = feval(f,x)/(N+1);                       % first kind of Chebyshev points
g = fft(fx([1:N+1  N+1:−1:1]));               % f  evaluated at these points
hx = real(exp(2
*
1i
*
pi
*
(0:2
*
n+1)/(4*n+4))
.
*
g’);   % FFT
a  = hx(1:n+1);a(1)=0.5*a(1);             % Chebyshev coefficientsFor the Filippi points, we presented a matlab code for computing the coefficientsan by FFT in Algorithm 3.Algorithm 3:
f
unction a = fejer2(f,N)                % coefficients forFejér’s  second rule
x = cos(pi
*
(0:N+2)’/(N+2));
fx = feval(f,x)/(2N+4);                % f evaluated at these points
g = fft(fx([1:N+3  N+2:−1:2]));            %  FFT
b = [g(1);  g(2:N+2)+g(2*N+4:−1:N+4);  g(N+3)];
b(N+1:
−
2:1) = b(N+1:-2:1)
−
2
*
b(N+3);
b(N:
−
2:1) = b(N:-2:1)
−
b(N+2);
b(1) = b(1)+mod(N+1,2)
*
b(N+3)+mod(N,2)
*
b(N+2)/2;
a = b(1:N+1);                     %  Chebyshev coefficientsIn addition, the coefficients of the interpolation polynomialsQN[f](x) at these Chebyshev-type points may be also efficiently evaluated by DCT and inverse discrete sine transform (IDST), respectively. The discrete cosine transform denoted by Y=dct(X) is closely related to the discrete Fourier transform but using purely real numbers and takes O(Nlog⁡N) operations for
(34)mmmY(k)=σ(k)∑s=1NX(s)cos⁡⁡((k-1)π(2s-1)2N)withσ(1)=1N,σ(k)=2N,for2≤k≤N.
The discrete sine transform denoted by Y=dst(X) and its inverse by X=idst(Y) both take O(Nlog⁡N) operations for
(35)Y(k)=∑s=1NX(s)sin⁡(kπsN+1).Notice that the coefficientsaj for the interpolation polynomial QN[f](x)=∑n=1N+1‍′an-1Tn-1(x) at {cos⁡((2k-1)π/(2N+2))}k=1N+1 are represented by
(36)an-1=2N+1∑s=1N+1‍f(cos⁡((2s-1)π2N+2))×cos⁡((2s-1)(n-1)π2N+2),n=1,2,…,N+1,
and a~n for the interpolation polynomial QN[f](x)=∑n=1N+1‍a~n-1Un-1(x) at {cos⁡(kπ/(N+2))}k=1N+1 satisfies
(37)f(cos⁡⁡(kπN+2))sin⁡(kπN+2)=∑s=1N+1a~s-1sin⁡(skπN+2),k=1,2,…,N+1.
Then both can be efficiently calculated by DCT and IDST, respectively. Thematlab codes on the DCT and IDST are very short and just need three rows. Now, we give the codes in Algorithm 4 and Algorithm 5.Algorithm 4:
f
unction a = fejer1dct(f,N)            % coefficients forFejér’s  second rule
x = cos(pi
*
(2
*
(0:N)’+1)/(2*N+2));         % The first kind of Chebyshev points
fx = feval(f,x);                   %f  evaluated at these points
a = dct(fx)
*
sqrt(2/(N+1));a(1)=a(1)/sqrt(2);    % Chebyshev coefficientsAlgorithm 5:
f
unction a = fejer2idst(f,N)            % coefficients forFejér’s  second rule
x  = cos(pi*(1:N+1)’/(N+2));             % Filippi points
fx = feval(f,x).
*
sin(pi
*
(1:N+1)’/(N+2));      % f evaluated at  these points
a = idst(fx);                     % Chebyshev coefficientsFor the first kind of Chebyshev points, amatlab code for computing the coefficients an by DCT is presented in Algorithm 4.As far as the Filippi points, amatlab code for computing the coefficients a~n by IDST is presented in Algorithm 5.Remark 6.
For the modified moments on the second kind Chebyshev polynomialM~n∶=∫-11‍w(x)Un(x)dx, we have
(38)M~0=M0,M~1=2M1,M~n=2Mn+M~n-2,n=2,3,…;
here we use the simple equation Un+2=2Tn+2+Un (see [31, pp. 778]).


## 3. Numerical Examples

In this section, we illustrate the accuracy and efficiency of the Clenshaw-Curtis, Fejér’s first- and second-type rules for the functionstan|x|, 1/(1+16x2) and |x-0.5|0.6 by the algorithms presented in this paper, which are compared with the Gauss-Jacobi quadrature used [x,w]=jacpts(n,α,β) inchebfun v4.2 [39] (see Figure 1). The first column computed by Gauss-Jacobi quadrature in Figure 1 takes 85.7386 seconds and the others totally take 2.9211 seconds in a Lenovo computer with Intel Core 3.20 GHz and 3.47 GB RAM. Figure 2 shows the convergence errors by the three quadratures, which takes 7.336958 seconds.Figure 1
The relative errors compared with Gauss quadrature for∫-11‍(1-x)α(1+x)βf(x)dx evaluated by the Clenshaw-Curtis, Fejér’s first- and second-type rules with N nodes, respectively: f(x)=tan|x| or 1/(1+16x2) and N=10:1000.Figure 2
The absolute errors for∫-11‍(1-x)α(1+x)βln⁡⁡((1+x)/2)f(x)dx evaluated by the Clenshaw-Curtis, Fejér’s first- and second-type rules with N nodes, respectively: f(x)=tan⁡|x| or |x-0.5|0.6 and N=10:1000.



---

*Source: 436164-2014-08-27.xml*