Recent Theory and Applications on Numerical Algorithms and Special FunctionsView this Special Issue
Research Article | Open Access
On Fast and Stable Implementation of Clenshaw-Curtis and Fejér-Type Quadrature Rules
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.
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  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 , Sloan and Smith [5, 6], Sommariva , Trefethen , Waldvogel , 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 . A streamlined MATLAB code is given as well in . In addition, Clenshaw and Curtis , O’Hara and Smith , Trefethen [8, 13], Xiang and Bornemann in , 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 : 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 .(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 , except in two cases: (ii)For , for the recurrence formula  is with where is the Beta function and is the Psi function (see Abramowitz and Stegun ). The forward recursion is as numerically stable as (4) except for (6) or (7) .
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  or Lozier’s algorithm  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  and Piessens and Branders ). 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 ). If , , and is times continuously differentiable for , then where
Lemma 2. If , , and is times continuously differentiable for , then where
Theorem 3. If , then where
Proof. For , taking in (3), we have
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
Remark 5. Piessens and Branders gave the first term of asymptotic expansion for and in  using the asymptotic theory of Fourier coefficients (Lighthill ). Furthermore, Piessens  presented the asymptotic expansion with explicit formulas for the first three terms for , which is equivalent to (20) except a minor clerical error. In , should be modified by . The first term of (25) is different from the first asymptotic term for proposed in ; 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  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 .