HOME

TheInfoList



OR:

Clenshaw–Curtis quadrature and Fejér quadrature are methods for
numerical integration In analysis, numerical integration comprises a broad family of algorithms for calculating the numerical value of a definite integral. The term numerical quadrature (often abbreviated to quadrature) is more or less a synonym for "numerical integr ...
, or "quadrature", that are based on an expansion of the integrand in terms of
Chebyshev polynomials The Chebyshev polynomials are two sequences of orthogonal polynomials related to the cosine and sine functions, notated as T_n(x) and U_n(x). They can be defined in several equivalent ways, one of which starts with trigonometric functions: ...
. Equivalently, they employ a
change of variables In mathematics, a change of variables is a basic technique used to simplify problems in which the original variables are replaced with functions of other variables. The intent is that when expressed in new variables, the problem may become si ...
x = \cos \theta and use a
discrete cosine transform A discrete cosine transform (DCT) expresses a finite sequence of data points in terms of a sum of cosine functions oscillating at different frequency, frequencies. The DCT, first proposed by Nasir Ahmed (engineer), Nasir Ahmed in 1972, is a widely ...
(DCT) approximation for the cosine series. Besides having fast-converging accuracy comparable to
Gaussian quadrature In numerical analysis, an -point Gaussian quadrature rule, named after Carl Friedrich Gauss, is a quadrature rule constructed to yield an exact result for polynomials of degree or less by a suitable choice of the nodes and weights for . Th ...
rules, Clenshaw–Curtis quadrature naturally leads to nested quadrature rules (where different accuracy orders share points), which is important for both adaptive quadrature and multidimensional quadrature ( cubature). Briefly, the function f(x) to be integrated is evaluated at the N extrema or roots of a Chebyshev polynomial and these values are used to construct a polynomial approximation for the function. This polynomial is then integrated exactly. In practice, the integration weights for the value of the function at each node are precomputed, and this computation can be performed in O(N \log N) time by means of
fast Fourier transform A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT). A Fourier transform converts a signal from its original domain (often time or space) to a representation in ...
-related algorithms for the DCT.W. Morven Gentleman, "Implementing Clenshaw-Curtis quadrature I: Methodology and experience," ''Communications of the ACM'' 15(5), p. 337-342 (1972).Jörg Waldvogel,
Fast construction of the Fejér and Clenshaw-Curtis quadrature rules
" ''BIT Numerical Mathematics'' 46 (1), p. 195-202 (2006).


General method

A simple way of understanding the algorithm is to realize that Clenshaw–Curtis quadrature (proposed by those authors in 1960)C. W. Clenshaw and A. R. Curtis
A method for numerical integration on an automatic computer
''Numerische Mathematik'' 2, 197 (1960).
amounts to integrating via a change of variable . The algorithm is normally expressed for integration of a function over the interval ��1,1(any other interval can be obtained by appropriate rescaling). For this integral, we can write: \int_^1 f(x)\,dx = \int_0^\pi f(\cos \theta) \sin(\theta)\, d\theta . That is, we have transformed the problem from integrating f(x) to one of integrating f(\cos \theta) \sin \theta. This can be performed if we know the cosine series for f(\cos \theta): f(\cos \theta) = \frac + \sum_^\infty a_k \cos (k\theta) in which case the integral becomes: \int_0^\pi f(\cos \theta) \sin(\theta)\, d\theta = a_0 + \sum_^\infty \frac . Of course, in order to calculate the cosine series coefficients a_k = \frac \int_0^\pi f(\cos \theta) \cos(k \theta)\, d\theta,\quad k=0,1,2,\dots, one must again perform a numeric integration, so at first this may not seem to have simplified the problem. Unlike computation of arbitrary integrals, however, Fourier-series integrations for
periodic functions A periodic function, also called a periodic waveform (or simply periodic wave), is a function that repeats its values at regular intervals or periods. The repeatable part of the function or waveform is called a ''cycle''. For example, the tr ...
(like f(\cos\theta), by construction), up to the
Nyquist frequency In signal processing, the Nyquist frequency (or folding frequency), named after Harry Nyquist, is a characteristic of a Sampling (signal processing), sampler, which converts a continuous function or signal into a discrete sequence. For a given S ...
k=N, are accurately computed by the N+1 equally spaced and equally weighted points \theta_n = n \pi / N for n = 0,\ldots,N (except the endpoints are weighted by 1/2, to avoid double-counting, equivalent to the
trapezoidal rule In calculus, the trapezoidal rule (or trapezium rule in British English) is a technique for numerical integration, i.e., approximating the definite integral: \int_a^b f(x) \, dx. The trapezoidal rule works by approximating the region under the ...
or the
Euler–Maclaurin formula In mathematics, the Euler–Maclaurin formula is a formula for the difference between an integral and a closely related sum. It can be used to approximate integrals by finite sums, or conversely to evaluate finite sums and infinite series using ...
). That is, we approximate the cosine-series integral by the type-I
discrete cosine transform A discrete cosine transform (DCT) expresses a finite sequence of data points in terms of a sum of cosine functions oscillating at different frequency, frequencies. The DCT, first proposed by Nasir Ahmed (engineer), Nasir Ahmed in 1972, is a widely ...
(DCT): a_k \approx \frac \left \frac + \frac (-1)^k + \sum_^ f(\cos[n\pi/N \cos(n k \pi/N) \right">\pi_N.html" ;"title="\frac + \frac (-1)^k + \sum_^ f(\cos[n\pi/N">\frac + \frac (-1)^k + \sum_^ f(\cos[n\pi/N \cos(n k \pi/N) \right/math> for k = 0,\ldots,N and then use the formula above for the integral in terms of these a_k. Because only a_ is needed, the formula simplifies further into a type-I DCT of order , assuming ''N'' is an even number: a_ \approx \frac \left[ \frac + f(0) (-1)^k + \sum_^ \left\ \cos\left(\frac\right) \right] From this formula, it is clear that the Clenshaw–Curtis quadrature rule is symmetric, in that it weights ''f''(''x'') and ''f''(−''x'') equally. Because of
aliasing In signal processing and related disciplines, aliasing is a phenomenon that a reconstructed signal from samples of the original signal contains low frequency components that are not present in the original one. This is caused when, in the ori ...
, one only computes the coefficients a_ up to , since discrete sampling of the function makes the frequency of indistinguishable from that of ''N''–2''k''. Equivalently, the a_ are the amplitudes of the unique bandlimited trigonometric interpolation polynomial passing through the points where is evaluated, and we approximate the integral by the integral of this interpolation polynomial. There is some subtlety in how one treats the a_ coefficient in the integral, however—to avoid double-counting with its alias it is included with weight 1/2 in the final approximate integral (as can also be seen by examining the interpolating polynomial): \int_0^\pi f(\cos \theta) \sin(\theta)\, d\theta \approx a_0 + \sum_^ \frac + \frac.


Connection to Chebyshev polynomials

The reason that this is connected to the Chebyshev polynomials T_k(x) is that, by definition, T_k(\cos\theta) = \cos(k\theta), and so the cosine series above is really an approximation of f(x) by Chebyshev polynomials: f(x) = \frac T_0(x) + \sum_^\infty a_k T_k(x), and thus we are "really" integrating f(x) by integrating its approximate expansion in terms of Chebyshev polynomials. The evaluation points x_n = \cos(n\pi/N) correspond to the extrema of the Chebyshev polynomial T_N(x). The fact that such
Chebyshev approximation In mathematics, approximation theory is concerned with how functions can best be approximated with simpler functions, and with quantitatively characterizing the errors introduced thereby. What is meant by ''best'' and ''simpler'' will depend ...
is just a cosine series under a change of variables is responsible for the rapid convergence of the approximation as more terms T_k (x) are included. A cosine series converges very rapidly for functions that are even, periodic, and sufficiently smooth. This is true here, since f(\cos\theta) is even and periodic in \theta by construction, and is ''k''-times differentiable everywhere if f(x) is ''k''-times differentiable on 1,1/math>. (In contrast, directly applying a cosine-series expansion to f(x) instead of f(\cos \theta) will usually ''not'' converge rapidly because the slope of the even-periodic extension would generally be discontinuous.)


Fejér quadrature

Fejér proposed two quadrature rules very similar to Clenshaw–Curtis quadrature, but much earlier (in 1933).Leopold 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'' 39 (1933), pp. 521–534. Leopold Fejér,
Mechanische Quadraturen mit positiven Cotesschen Zahlen
''Mathematische Zeitschrift'' 37, 287 (1933).
Of these two, Fejér's "second" quadrature rule is nearly identical to Clenshaw–Curtis. The only difference is that the endpoints f(-1) and f(1) are set to zero. That is, Fejér only used the ''interior'' extrema of the Chebyshev polynomials, i.e. the true stationary points. Fejér's "first" quadrature rule evaluates the a_k by evaluating f(\cos\theta) at a different set of equally spaced points, halfway between the extrema: \theta_n = (n + 0.5) \pi / N for 0 \leq n < N. These are the ''roots'' of T_N(\cos\theta), and are known as the
Chebyshev nodes In numerical analysis, Chebyshev nodes (also called Chebyshev points or a Chebyshev grid) are a set of specific algebraic numbers used as nodes for polynomial interpolation and numerical integration. They are the Projection (linear algebra), pr ...
. (These equally spaced midpoints are the only other choice of quadrature points that preserve both the even symmetry of the cosine transform and the translational symmetry of the periodic Fourier series.) This leads to a formula: a_k \approx \frac \sum_^ f(\cos n+0.5)\pi/N \cos n+0.5) k \pi/N which is precisely the type-II DCT. However, Fejér's first quadrature rule is not nested: the evaluation points for do not coincide with any of the evaluation points for , unlike Clenshaw–Curtis quadrature or Fejér's second rule. Despite the fact that Fejér discovered these techniques before Clenshaw and Curtis, the name "Clenshaw–Curtis quadrature" has become standard.


Comparison to Gaussian quadrature

The classic method of
Gaussian quadrature In numerical analysis, an -point Gaussian quadrature rule, named after Carl Friedrich Gauss, is a quadrature rule constructed to yield an exact result for polynomials of degree or less by a suitable choice of the nodes and weights for . Th ...
evaluates the integrand at N+1 points and is constructed to ''exactly'' integrate polynomials up to degree 2N+1. In contrast, Clenshaw–Curtis quadrature, above, evaluates the integrand at N+1 points and exactly integrates polynomials only up to degree N+1 (assuming N is an even number). It may seem, therefore, that Clenshaw–Curtis is intrinsically worse than Gaussian quadrature, but in reality this does not seem to be the case. In practice, several authors have observed that Clenshaw–Curtis can have accuracy comparable to that of Gaussian quadrature for the same number of points. This is possible because most numeric integrands are not polynomials (especially since polynomials can be integrated analytically), and approximation of many functions in terms of Chebyshev polynomials converges rapidly (see
Chebyshev approximation In mathematics, approximation theory is concerned with how functions can best be approximated with simpler functions, and with quantitatively characterizing the errors introduced thereby. What is meant by ''best'' and ''simpler'' will depend ...
). In fact, recent theoretical results argue that both Gaussian and Clenshaw–Curtis quadrature have error bounded by O( N/k) for a ''k''-times differentiable integrand. One often cited advantage of Clenshaw–Curtis quadrature is that the quadrature weights can be evaluated in O(N \log N) time by
fast Fourier transform A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT). A Fourier transform converts a signal from its original domain (often time or space) to a representation in ...
algorithms (or their analogues for the DCT), whereas most algorithms for Gaussian quadrature weights required O(N^2) time to compute. However, recent algorithms have attained O(N) complexity for Gauss–Legendre quadrature. As a practical matter, high-order numeric integration is rarely performed by simply evaluating a quadrature formula for very large N. Instead, one usually employs an adaptive quadrature scheme that first evaluates the integral to low order, and then successively refines the accuracy by increasing the number of sample points, possibly only in regions where the integral is inaccurate. To evaluate the accuracy of the quadrature, one compares the answer with that of a quadrature rule of even lower order. Ideally, this lower-order quadrature rule evaluates the integrand at a ''subset'' of the original ''N'' points, to minimize the integrand evaluations. This is called a nested quadrature rule, and here Clenshaw–Curtis has the advantage that the rule for order ''N'' uses a subset of the points from order 2''N''. In contrast, Gaussian quadrature rules are not naturally nested, and so one must employ Gauss–Kronrod quadrature formulas or similar methods. Nested rules are also important for sparse grids in multidimensional quadrature, and Clenshaw–Curtis quadrature is a popular method in this context.


Integration with weight functions

More generally, one can pose the problem of integrating an arbitrary f(x) against a fixed ''weight function'' w(x) that is known ahead of time: \int_^1 f(x) w(x)\,dx = \int_0^\pi f(\cos \theta) w(\cos\theta) \sin(\theta)\, d\theta . The most common case is w(x) = 1, as above, but in certain applications a different weight function is desirable. The basic reason is that, since w(x) can be taken into account ''a priori'', the integration error can be made to depend only on the accuracy in approximating f(x), regardless of how badly behaved the weight function might be. Clenshaw–Curtis quadrature can be generalized to this case as follows. As before, it works by finding the cosine-series expansion of f(\cos \theta) via a DCT, and then integrating each term in the cosine series. Now, however, these integrals are of the form W_k = \int_0^\pi w(\cos \theta) \cos(k \theta) \sin(\theta)\, d\theta . For most w(x), this integral cannot be computed analytically, unlike before. Since the same weight function is generally used for many integrands f(x), however, one can afford to compute these W_k numerically to high accuracy beforehand. Moreover, since w(x) is generally specified analytically, one can sometimes employ specialized methods to compute W_k. For example, special methods have been developed to apply Clenshaw–Curtis quadrature to integrands of the form f(x) w(x) with a weight function w(x) that is highly oscillatory, e.g. a
sinusoid A sine wave, sinusoidal wave, or sinusoid (symbol: ∿) is a periodic wave whose waveform (shape) is the trigonometric sine function. In mechanics, as a linear motion over time, this is '' simple harmonic motion''; as rotation, it correspond ...
or
Bessel function Bessel functions, named after Friedrich Bessel who was the first to systematically study them in 1824, are canonical solutions of Bessel's differential equation x^2 \frac + x \frac + \left(x^2 - \alpha^2 \right)y = 0 for an arbitrary complex ...
(see, e.g., Evans & Webster, 1999G. A. Evans and J. R. Webster, "A comparison of some methods for the evaluation of highly oscillatory integrals," ''Journal of Computational and Applied Mathematics'', vol. 112, p. 55-69 (1999).). This is useful for high-accuracy
Fourier series A Fourier series () is an Series expansion, expansion of a periodic function into a sum of trigonometric functions. The Fourier series is an example of a trigonometric series. By expressing a function as a sum of sines and cosines, many problems ...
and Fourier–Bessel series computation, where simple w(x) = 1 quadrature methods are problematic because of the high accuracy required to resolve the contribution of rapid oscillations. Here, the rapid-oscillation part of the integrand is taken into account via specialized methods for W_k, whereas the unknown function f(x) is usually better behaved. Another case where weight functions are especially useful is if the integrand is unknown but has a known singularity of some form, e.g. a known discontinuity or integrable divergence (such as ) at some point. In this case the singularity can be pulled into the weight function w(x) and its analytical properties can be used to compute W_k accurately beforehand. Note that
Gaussian quadrature In numerical analysis, an -point Gaussian quadrature rule, named after Carl Friedrich Gauss, is a quadrature rule constructed to yield an exact result for polynomials of degree or less by a suitable choice of the nodes and weights for . Th ...
can also be adapted for various weight functions, but the technique is somewhat different. In Clenshaw–Curtis quadrature, the integrand is always evaluated at the same set of points regardless of w(x), corresponding to the extrema or roots of a Chebyshev polynomial. In Gaussian quadrature, different weight functions lead to different
orthogonal polynomials In mathematics, an orthogonal polynomial sequence is a family of polynomials such that any two different polynomials in the sequence are orthogonal In mathematics, orthogonality (mathematics), orthogonality is the generalization of the geom ...
, and thus different roots where the integrand is evaluated.


Integration on infinite and semi-infinite intervals

It is also possible to use Clenshaw–Curtis quadrature to compute integrals of the form \int_0^\infty f(x)\,dx and \int_^\infty f(x)\,dx, using a coordinate-remapping technique.John P. Boyd, "Exponentially convergent Fourier–Chebshev /nowiki>''sic''/nowiki> quadrature schemes on bounded and infinite intervals," ''J. Scientific Computing'' 2 (2), p. 99-109 (1987). High accuracy, even exponential convergence for smooth integrands, can be retained as long as f(x) decays sufficiently quickly as , ''x'', approaches infinity. One possibility is to use a generic coordinate transformation such as \int_^f(x)\,dx = \int_^ f\left(\frac\right)\frac\,dt, to transform an infinite or semi-infinite interval into a finite one, as described in
Numerical integration In analysis, numerical integration comprises a broad family of algorithms for calculating the numerical value of a definite integral. The term numerical quadrature (often abbreviated to quadrature) is more or less a synonym for "numerical integr ...
. There are also additional techniques that have been developed specifically for Clenshaw–Curtis quadrature. For example, one can use the coordinate remapping x = L \cot^2(\theta/2), where ''L'' is a user-specified constant (one could simply use ''L''=1; an optimal choice of ''L'' can speed convergence, but is problem-dependent), to transform the semi-infinite integral into: \int_0^\infty f(x)\,dx = 2L \int_0^\pi \frac \sin(\theta)\,d\theta . The factor multiplying sin(''θ''), ''f''(...)/(...)2, can then be expanded in a cosine series (approximately, using the discrete cosine transform) and integrated term-by-term, exactly as was done for ''f''(cos ''θ'') above. To eliminate the singularity at ''θ''=0 in this integrand, one merely requires that ''f''(''x'') go to zero sufficiently fast as ''x'' approaches infinity, and in particular ''f''(''x'') must decay at least as fast as 1/''x''3/2. For a doubly infinite interval of integration, one can use the coordinate remapping x = L\cot(\theta) (where ''L'' is a user-specified constant as above) to transform the integral into: \int_^\infty f(x)\,dx = L \int_0^\pi \frac\,d\theta \approx \frac \sum_^ \frac. In this case, we have used the fact that the remapped integrand is already periodic and so can be directly integrated with high (even exponential) accuracy using the trapezoidal rule (assuming ''f'' is sufficiently smooth and rapidly decaying); there is no need to compute the cosine series as an intermediate step. Note that the quadrature rule does not include the endpoints, where we have assumed that the integrand goes to zero. The formula above requires that ''f''(''x'') decay faster than 1/''x''2 as ''x'' goes to ±∞. (If ''f'' decays exactly as 1/''x''2, then the integrand goes to a finite value at the endpoints and these limits must be included as endpoint terms in the trapezoidal rule.). However, if ''f'' decays only polynomially quickly, then it may be necessary to use a further step of Clenshaw–Curtis quadrature to obtain exponential accuracy of the remapped integral instead of the trapezoidal rule, depending on more details of the limiting properties of ''f'': the problem is that, although is indeed periodic with period Ï€, it is not necessarily smooth at the endpoints if all the derivatives do not vanish there .g. the function ''f''(''x'') = tanh(''x''3)/''x''3 decays as 1/''x''3 but has a jump discontinuity in the slope of the remapped function at θ=0 and Ï€ Another coordinate-remapping approach was suggested for integrals of the form \int_0^\infty e^ g(x)\,dx, in which case one can use the transformation x = -\ln 1 + \cos\theta)/2/math> to transform the integral into the form \int_0^\pi f(\cos\theta)\sin\theta \,d\theta where f(u) = g(-\ln 1+u)/2/2, at which point one can proceed identically to Clenshaw–Curtis quadrature for ''f'' as above.Nirmal Kumar Basu and Madhav Chandra Kundu, "Some methods of numerical integration over a semi-infinite interval," ''Applications of Mathematics'' 22 (4), p. 237-243 (1977). Because of the endpoint singularities in this coordinate remapping, however, one uses Fejér's first quadrature rule hich does not evaluate ''f''(−1)unless ''g''(∞) is finite.


Precomputing the quadrature weights

In practice, it is inconvenient to perform a DCT of the sampled function values ''f''(cos Î¸) for each new integrand. Instead, one normally precomputes quadrature weights w_n (for ''n'' from 0 to ''N''/2, assuming that ''N'' is even) so that \int_^1 f(x)\,dx \approx \sum_^ w_n \left\ . These weights w_n are also computed by a DCT, as is easily seen by expressing the computation in terms of
matrix Matrix (: matrices or matrixes) or MATRIX may refer to: Science and mathematics * Matrix (mathematics), a rectangular array of numbers, symbols or expressions * Matrix (logic), part of a formula in prenex normal form * Matrix (biology), the m ...
algebra Algebra is a branch of mathematics that deals with abstract systems, known as algebraic structures, and the manipulation of expressions within those systems. It is a generalization of arithmetic that introduces variables and algebraic ope ...
. In particular, we computed the cosine series coefficients a_ via an expression of the form: c = \begin a_0 \\ a_2 \\ a_4 \\ \vdots \\ a_N \end = D \begin y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_ \end = Dy, where ''D'' is the matrix form of the (''N''/2+1)-point type-I DCT from above, with entries (for zero-based indices): D_ = \frac \cos\left(\frac\right) \times \begin 1/2 & n=0,N/2 \\ 1 & \mathrm \end and y_n is y_n = f(\cos \pi/N + f(-\cos \pi/N . As discussed above, because of
aliasing In signal processing and related disciplines, aliasing is a phenomenon that a reconstructed signal from samples of the original signal contains low frequency components that are not present in the original one. This is caused when, in the ori ...
, there is no point in computing coefficients beyond a_N, so ''D'' is an (N/2+1)\times(N/2+1) matrix. In terms of these coefficients ''c'', the integral is approximately: \int_^1 f(x)\,dx \approx a_0 + \sum_^ \frac + \frac = d^T c, from above, where ''c'' is the vector of coefficients a_ above and ''d'' is the vector of integrals for each Fourier coefficient: d = \begin 1 \\ 2/(1-4) \\ 2/(1-16) \\ \vdots \\ 2 / (1- -22) \\ 1 / (1-N^2) \end. (Note, however, that these weight factors are altered if one changes the DCT matrix ''D'' to use a different normalization convention. For example, it is common to define the type-I DCT with additional factors of 2 or factors in the first and last rows or columns, which leads to corresponding alterations in the ''d'' entries.) The d^T c summation can be re-arranged to: \int_^1 f(x)\,dx \approx d^T c = d^T D y = (D^T d)^T y = w^T y where is the vector of the desired weights w_n above, with:w = D^T d. Since the
transpose In linear algebra, the transpose of a Matrix (mathematics), matrix is an operator which flips a matrix over its diagonal; that is, it switches the row and column indices of the matrix by producing another matrix, often denoted by (among other ...
d matrix D^T is also a DCT (e.g., the transpose of a type-I DCT is a type-I DCT, possibly with a slightly different normalization depending on the conventions that are employed), the quadrature weights ''w'' can be precomputed in ''O''(''N'' log ''N'') time for a given ''N'' using fast DCT algorithms. The weights w_n are positive and their sum is equal to one.J. P. Imhof, "On the Method for Numerical Integration of Clenshaw and Curtis", ''Numerische Mathematik'' 5, p. 138-141 (1963).


See also

*
Euler–Maclaurin formula In mathematics, the Euler–Maclaurin formula is a formula for the difference between an integral and a closely related sum. It can be used to approximate integrals by finite sums, or conversely to evaluate finite sums and infinite series using ...
* Gauss–Kronrod quadrature formula


References

{{DEFAULTSORT:Clenshaw-Curtis Quadrature Numerical integration