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

Xiang, Shuhuang; He, Guo; Wang, Haiyong

doi:https://doi.org/10.1155/2014/436164

Abstract and Applied Analysis

On this page

Abstract Introduction References Copyright Related Articles

Special Issue

Recent Theory and Applications on Numerical Algorithms and Special Functions

View this Special Issue

Research Article | Open Access

Volume 2014 | Article ID 436164 | https://doi.org/10.1155/2014/436164

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

Shuhuang Xiang,¹Guo He,¹and Haiyong Wang²

Academic Editor: Robert A. Van Gorder

Received27 Jun 2014

Accepted05 Aug 2014

Published27 Aug 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.

1. Introduction

The interpolation quadrature of the Clenshaw-Curtis rule as well as Fejér-type formulas for has been extensively studied since Fejér [1, 2] in 1933 and Clenshaw and Curtis [3] in 1960, where the nodes are of Chebyshev types while the weights are computed by sums of trigonometric functions. When (), this quadrature is called Fejér’s first-type rule. This kind of points is called the first kind of Chebyshev points, while Fejér’s second-type rule is corresponding to the Filippi points () and the Clenshaw-Curtis-type quadrature to the Clenshaw-Curtis points (the second kind of Chebyshev points) (). For more details, see Davis and Robinowitz [4], Sloan and Smith [5, 6], Sommariva [7], Trefethen [8], Waldvogel [9], and so forth.

In the case , 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 nodes by means of a discrete cosine transformation (DCT). An independent approach along the same lines, unified algorithms based on DFTs of order for generating the weights of the two Fejér rules and of the Clenshaw-Curtis rule, was presented in Waldvogel [9]. A streamlined MATLAB 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 functions and 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 of at the Clenshaw-Curtis points together with the modified moments [24]: where is the Chebyshev polynomial of degree and 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 satisfy the following recurrence formulas for Jacobi weights or Jacobi weights multiplied by [20].(i)For , by using Fasenmyer’s technique, the recurrence formula for the evaluation of the modified moments is with The forward recursion is numerically stable [20], except in two cases: (ii)For , for the recurrence formula [20] is with where is the Beta function and 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 are operations.

However, in case (6) or (7), the accuracy of the forward recursion is catastrophic particularly when or and (see Table 1). In case (6) the relative errors of the computed values obtained by the forward recursion behave approximately as and in case (7) as

In this paper, we will consider interpolation approaches for Clenshaw-Curtis rules as well as Fejér’s first- and second-type formulas for with or , which can be efficiently calculated in operations. Computing the modified moments 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 Section 2.1, we studied the asymptotic expansions of the modified moments . Based on the results of the asymptotic expansions, we discussed Oliver’s algorithm for the modified moments 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 concise MATLAB codes for evaluation of the coefficients 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 can be efficiently and stably computed then can be so. Note that the asymptotic behaviour of two linearly independent solutions of (4) satisfies (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]). If , , and is times continuously differentiable for , then where

Lemma 2. If , , and is times continuously differentiable for , then where

Proof. The proof can be directly derived from the proof of Lemma 1 and the proof of the Theorem 5 in Erdélyi [36].

Theorem 3. If , then where

Proof. For , taking in (3), we have where Consequently, the desired result can be derived by applying (16) to (22).
The above outcome can be extended to the case of , , since can be written as where denotes the largest integer less than or equal to .

Theorem 4. If , then where

Proof. Letting in (8) and using , we have 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 for and 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 , which is equivalent to (20) except a minor clerical error. In [18], should be modified by . The first term of (25) is different from the first asymptotic term for proposed in [18]; that is, . 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

Letwhere “” denotes the transpose; then the modified moments can be solved by where is computed by hypergeometric function [20] when : If , 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 operations.

In addition, for the weight 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 and one end . Consequently, the modified moments can be solved by where The end value can be calculated by its asymptotic formula (25) with the first four-term truncation.

Tables 2 and 3 show the accuracy of Oliver’s algorithm for and , and , respectively. Table 4 displays the accuracy of Oliver’s algorithm for and , and for . Table 5 displays the accuracy of Oliver’s algorithm for and , and for . 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. The MATLAB codes on Oliver’s algorithms and all the MATLAB codes in this paper can be downloaded from [38].

2.3. Fast Computation of the Coefficients by FFT

The coefficients for the interpolation polynomial at the three kinds of points can be efficiently computed by FFT. For the Clenshaw-Curtis points, the elegant MATLAB code on the is from [8] (see Algorithm 1).

function 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(2N:−1:N+2); g(N+1)]; % Chebyshev coefficients

For 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 following MATLAB functions.

For the first kind of Chebyshev points, we presented a MATLAB code for computing the coefficients by FFT in Algorithm 2.

function a = fejer1(f,N)
x = cos(pi(2(0:N)’+1)/(2*N+2)); % coefficients for Fejér’s first rule
fx = 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(21ipi(0:2n+1)/(4*n+4)).g’); % FFT
a = hx(1:n+1);a(1)=0.5*a(1); % Chebyshev coefficients

For the Filippi points, we presented a matlab code for computing the coefficients by FFT in Algorithm 3.

function a = fejer2(f,N) % coefficients for Fejé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(2N+4:−1:N+4); g(N+3)];
b(N+1:−2:1) = b(N+1:-2:1)−2b(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 coefficients

In addition, the coefficients of the interpolation polynomials 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 is closely related to the discrete Fourier transform but using purely real numbers and takes operations for The discrete sine transform denoted by and its inverse by both take operations for

Notice that the coefficients for the interpolation polynomial at are represented by and for the interpolation polynomial at satisfies Then both can be efficiently calculated by DCT and IDST, respectively. The MATLAB 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.

function a = fejer1dct(f,N) % coefficients for Fejé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 coefficients

function a = fejer2idst(f,N) % coefficients for Fejé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 coefficients

For the first kind of Chebyshev points, a MATLAB code for computing the coefficients by DCT is presented in Algorithm 4.

As far as the Filippi points, a MATLAB code for computing the coefficients by IDST is presented in Algorithm 5.

Remark 6. For the modified moments on the second kind Chebyshev polynomial , we have here we use the simple equation (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 functions , and by the algorithms presented in this paper, which are compared with the Gauss-Jacobi quadrature used in CHEBFUN 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This paper is supported by the National Natural Science Foundation of China, no. 11371376.

References

L. Fejér, “On the infinite sequences arising in the theories of harmonic analysis, of interpolation, and of mechanical quadratures,” Bulletin of the American Mathematical Society, vol. 39, no. 8, pp. 521–534, 1933.
View at: Publisher Site | Google Scholar | MathSciNet
L. Fejér, “Mechanische Quadraturen mit positiven Cotesschen Zahlen,” Mathematische Zeitschrift, vol. 37, no. 1, pp. 287–309, 1933.
View at: Publisher Site | Google Scholar | MathSciNet
C. W. Clenshaw and A. R. Curtis, “A method for numerical integration on an automatic computer,” Numerische Mathematik, vol. 2, pp. 197–205, 1960.
View at: Publisher Site | Google Scholar | MathSciNet
P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, Academic Press, New York, NY, USA, 2nd edition, 1984.
View at: MathSciNet
I. H. Sloan and W. E. Smith, “Product-integration with the Clenshaw-CURtis and related points,” Numerische Mathematik, vol. 30, no. 4, pp. 415–428, 1978.
View at: Publisher Site | Google Scholar | MathSciNet
I. H. Sloan and W. E. Smith, “Product integration with the Clenshaw-CURtis points: implementation and error estimates,” Numerische Mathematik, vol. 34, no. 4, pp. 387–401, 1980.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
A. Sommariva, “Fast construction of Fejér and Clenshaw-Curtis rules for general weight functions,” Computers & Mathematics with Applications, vol. 65, no. 4, pp. 682–693, 2013.
View at: Publisher Site | Google Scholar | MathSciNet
L. N. Trefethen, “Is Gauss quadrature better than Clenshaw-Curtis?” SIAM Review, vol. 50, no. 1, pp. 67–87, 2008.
View at: Publisher Site | Google Scholar | MathSciNet
J. Waldvogel, “Fast construction of the Fejér and Clenshaw-Curtis quadrature rules,” BIT, vol. 46, no. 1, pp. 195–202, 2006.
View at: Publisher Site | Google Scholar | MathSciNet
W. M. Gentleman, “Implementing Clenshaw-Curtis quadrature,” Communications of the ACM, vol. 15, no. 5, pp. 337–346, 1972.
View at: Google Scholar
W. Morven, “Algorithm 424: Clenshaw-Curtis quadrature [D1],” Communications of the ACM, vol. 15, no. 5, pp. 353–355, 1972.
View at: Google Scholar
H. O'Hara and F. J. Smith, “Error estimation in the Clenshaw-CURtis quadrature formula,” The Computer Journal, vol. 11, pp. 213–219, 1968.
View at: Publisher Site | Google Scholar | MathSciNet
L. N. Trefethen, Approximation Theory and Approximation Practice, SIAM, New York, NY, USA, 2013.
View at: MathSciNet
S. Xiang and F. Bornemann, “On the convergence rates of Gauss and Clenshaw-CURtis quadrature for functions of limited regularity,” SIAM Journal on Numerical Analysis, vol. 50, no. 5, pp. 2581–2587, 2012.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
S. Xiang, X. Chen, and H. Wang, “Error bounds for approximation in Chebyshev points,” Numerische Mathematik, vol. 116, no. 3, pp. 463–491, 2010.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
S. Xiang, “On convergence rates of Fejér and Gauss-Chebyshev quadrature rules,” Journal of Mathematical Analysis and Applications, vol. 405, no. 2, pp. 687–699, 2013.
View at: Publisher Site | Google Scholar | MathSciNet
S. Xiang, “On the optimal rates of convergence for quadratures derived from Chebyshev points,” http://arxiv.org/abs/1308.4322.
View at: Google Scholar
R. Piessens, “Modified Clenshaw-Curtis integration and applications to numerical computation of integral transforms,” in Numerical Integration, P. Keast and G. Fairweather, Eds., pp. 35–51, 1987.
View at: Google Scholar | MathSciNet
R. Piessens, “Computing integral transforms and solving integral equations using Chebyshev polynomial approximations,” Journal of Computational and Applied Mathematics, vol. 121, no. 1-2, pp. 113–124, 2000.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
R. Piessens and M. Branders, “The evaluation and application of some modified moments,” BIT Numerical Mathematics, vol. 13, no. 4, pp. 443–450, 1973.
View at: Publisher Site | Google Scholar | MathSciNet
R. Piessens and M. Branders, “Modified Clenshaw-CURtis method for the computation of Bessel function integrals,” BIT Numerical Mathematics, vol. 23, no. 3, pp. 370–381, 1983.
View at: Publisher Site | Google Scholar | MathSciNet
R. Piessens and M. Branders, “Computation of Fourier transform integrals using Chebyshev series expansions,” Computing, vol. 32, no. 2, pp. 177–186, 1984.
View at: Publisher Site | Google Scholar | MathSciNet
R. Piessens and M. Branders, “On the computation of Fourier transforms of singular functions,” Journal of Computational and Applied Mathematics, vol. 43, no. 1-2, pp. 159–169, 1992.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
W. Gautschi, “Computational aspects of three-term recurrence relations,” SIAM Review, vol. 9, pp. 24–82, 1967.
View at: Publisher Site | Google Scholar | MathSciNet
R. Chen and C. An, “On evaluation of Bessel transforms with oscillatory and algebraic singular integrands,” Journal of Computational and Applied Mathematics, vol. 264, pp. 71–81, 2014.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
R. Chen and C. An, “On the evaluation of infinite integrals involving Bessel functions,” Applied Mathematics and Computation, vol. 235, pp. 212–220, 2014.
View at: Publisher Site | Google Scholar | MathSciNet
H. Kang and S. Xiang, “Efficient integration for a class of highly oscillatory integrals,” Applied Mathematics and Computation, vol. 218, no. 7, pp. 3553–3564, 2011.
View at: Publisher Site | Google Scholar | MathSciNet
H. Kang and S. Xiang, “Efficient quadrature of highly oscillatory integrals with algebraic singularities,” Journal of Computational and Applied Mathematics, vol. 237, no. 1, pp. 576–588, 2013.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
H. Wang and S. Xiang, “Uniform approximations to Cauchy principal value integrals of oscillatory functions,” Applied Mathematics and Computation, vol. 215, no. 5, pp. 1886–1894, 2009.
View at: Publisher Site | Google Scholar | MathSciNet
S. Xiang, Y. J. Cho, H. Wang, and H. Brunner, “Clenshaw-Curtis-Filon-type methods for highly oscillatory Bessel transforms and applications,” IMA Journal of Numerical Analysis, vol. 31, no. 4, pp. 1281–1314, 2011.
View at: Publisher Site | Google Scholar | MathSciNet
M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, USA, 1964.
J. Oliver, “The numerical solution of linear recurrence relations,” Numerische Mathematik, vol. 11, pp. 349–360, 1968.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
D. W. Lozier, “Numerical solution of linear difference equations,” Tech. Rep. NBSIR 80-1976, Math. Analysis, Nat. Bureau of Standards, 1980.
View at: Google Scholar
J. Denef and R. Piessens, “The asymptotic behaviour of solutions of difference equations of Poincaré's type,” Bulletin de la Société Mathématique de Belgique, vol. 26, no. 2, pp. 133–347, 1974.
View at: Google Scholar | MathSciNet
A. Erdélyi, “Asymptotic representations of Fourier integrals and the method of stationary phase,” Journal of the Society for Industrial and Applied Mathematics, vol. 3, pp. 17–27, 1955.
View at: Google Scholar | MathSciNet
A. Erdélyi, “Asymptotic expansions of integrals involving logarithmic singularities,” Journal of the Society for Industrial and Applied Mathematics, vol. 1, pp. 38–47, 1956.
View at: Google Scholar
M. J. Lighthill, Fourier Analysis and Generalized Functions, Cambridge University Press, Cambridge, UK, 1959.
http://math.csu.edu.cn/office/teacherpage.aspx?namenumber=56.
L. N. Trefethen et al., “Chebfun Version 4.2,” The Chebfun Development Team, 2011, http://www.maths.ox.ac.uk/chebfun/.
View at: Google Scholar

Copyright

Copyright © 2014 Shuhuang Xiang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

3039

Downloads

1679

Citations