In practice, the integration weights for the value of the function at each node are precomputed, and this computation can be performed in
time by means of fast Fourier transform-related algorithms for the DCT.
[1][2] A simple way of understanding the algorithm is to realize that Clenshaw–Curtis quadrature (proposed by those authors in 1960)[3] amounts to integrating via a change of variable x = cos(θ).
(except the endpoints are weighted by 1/2, to avoid double-counting, equivalent to the trapezoidal rule or the Euler–Maclaurin formula).
[4][5] That is, we approximate the cosine-series integral by the type-I discrete cosine transform (DCT):
is needed, the formula simplifies further into a type-I DCT of order N/2, assuming N is an even number:
From this formula, it is clear that the Clenshaw–Curtis quadrature rule is symmetric, in that it weights f(x) and f(−x) equally.
A cosine series converges very rapidly for functions that are even, periodic, and sufficiently smooth.
That is, Fejér only used the interior extrema of the Chebyshev polynomials, i.e. the true stationary points.
at a different set of equally spaced points, halfway between the extrema:
Despite the fact that Fejér discovered these techniques before Clenshaw and Curtis, the name "Clenshaw–Curtis quadrature" has become standard.
In practice, several authors have observed that Clenshaw–Curtis can have accuracy comparable to that of Gaussian quadrature for the same number of points.
In fact, recent theoretical results[7] argue that both Gaussian and Clenshaw–Curtis quadrature have error bounded by
[8] As a practical matter, high-order numeric integration is rarely performed by simply evaluating a quadrature formula for very large
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.
can be taken into account a priori, the integration error can be made to depend only on the accuracy in approximating
For example, special methods have been developed to apply Clenshaw–Curtis quadrature to integrands of the form
quadrature methods are problematic because of the high accuracy required to resolve the contribution of rapid oscillations.
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 1/√x) at some point.
Note that Gaussian quadrature can also be adapted for various weight functions, but the technique is somewhat different.
[11] High accuracy, even exponential convergence for smooth integrands, can be retained as long as
to transform an infinite or semi-infinite interval into a finite one, as described in Numerical integration.
, 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[11]), to transform the semi-infinite integral into:
In this case, we have used the fact that the remapped integrand f(L cot θ)/sin2(θ) 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.
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 f(L cotθ)/sin2(θ) is indeed periodic with period π, it is not necessarily smooth at the endpoints if all the derivatives do not vanish there [e.g. the function f(x) = tanh(x3)/x3 decays as 1/x3 but has a jump discontinuity in the slope of the remapped function at θ=0 and π].
[12] Because of the endpoint singularities in this coordinate remapping, however, one uses Fejér's first quadrature rule [which does not evaluate f(−1)] unless g(∞) is finite.
In practice, it is inconvenient to perform a DCT of the sampled function values f(cos θ) for each new integrand.
where D is the matrix form of the (N/2+1)-point type-I DCT from above, with entries (for zero-based indices):
(Note, however, that these weight factors are altered if one changes the DCT matrix D to use a different normalization convention.