Quadratura de Clenshaw-Curtis: diferència entre les revisions

:<math>\int_0^\pi f(\cos \theta) \sin(\theta)\, d\theta \approx a_0 + \sum_{k=1}^{N/2-1} \frac{2 a_{2k}}{1 - (2k)^2} + \frac{a_{N}}{1 - N^2}.</math>

 

The reason that this is connected to the Chebyshev polynomials <math>T_k(x)</math> is that, by definition, <math>T_k(\cos\theta) = \cos(k\theta)</math>, and so the cosine series above is really an approximation of <math>f(x)</math> by Chebyshev polynomials:

 

The fact that such [[Chebyshev approximation]] is just a cosine series under a change of variables is responsible for the rapid convergence of the approximation as more terms <math>T_k(x)</math> are included. A cosine series converges very rapidly for functions that are [[even and odd functions|even]], periodic, and sufficiently smooth. This is true here, since <math>f(\cos\theta)</math> is even and periodic in <math>\theta</math> by construction, and is ''k''-times differentiable everywhere if <math>f(x)</math> is ''k''-times differentiable on <math>[-1,1]</math>. (In contrast, directly applying a cosine-series expansion to <math>f(x)</math> instead of <math>f(\cos \theta)</math> will usually ''not'' converge rapidly because the slope of the even-periodic extension would generally be discontinuous.)

 

[[Lipót Fejér|Fejér]] proposed two quadrature rules very similar to Clenshaw–Curtis quadrature, but much earlier (in 1933).<ref name=Fejer33>Leopold Fejér, "[http://projecteuclid.org/euclid.bams/1183496842 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, "[http://www.digizeitschriften.de/resolveppn/GDZPPN002374498 Mechanische Quadraturen mit positiven Cotesschen Zahlen], ''Mathematische Zeitschrift'' '''37''', 287 (1933).</ref>

 

Despite the fact that Fejér discovered these techniques before Clenshaw and Curtis, the name "Clenshaw–Curtis quadrature" has become standard.

 

The classic method of [[Gaussian quadrature]] evaluates the integrand at <math>N+1</math> points and is constructed to ''exactly'' integrate polynomials up to [[degree of a polynomial|degree]] <math>2N+1</math>. In contrast, Clenshaw–Curtis quadrature, above, evaluates the integrand at <math>N+1</math> points and exactly integrates polynomials only up to degree <math>N</math>. It may seem, therefore, that Clenshaw–Curtis is intrinsically worse than Gaussian quadrature, but in reality this does not seem to be the case.

 

One often cited advantage of Clenshaw–Curtis quadrature is that the quadrature weights can be evaluated in <math>O(N \log N)</math> time by [[fast Fourier transform]] algorithms (or their analogues for the DCT), whereas the Gaussian quadrature weights require <math>O(N^2)</math> time to compute. As a practical matter, however, high-order numeric integration is rarely performed by simply evaluating a quadrature formula for very large <math>N</math>. 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 formula]]s or similar methods. Nested rules are also important for [[sparse grid]]s in multidimensional quadrature, and Clenshaw–Curtis quadrature is a popular method in this context.<ref>Erich Novak and Klaus Ritter, "High dimensional integration of smooth functions over cubes," ''Numerische Mathematik'' vol. '''75''', pp. 79–97 (1996).</ref>

 

More generally, one can pose the problem of integrating an arbitrary <math>f(x)</math> against a fixed ''weight function'' <math>w(x)</math> that is known ahead of time:

 

Note that [[Gaussian quadrature]] 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 <math>w(x)</math>, corresponding to the extrema or roots of a Chebyshev polynomial. In Gaussian quadrature, different weight functions lead to different [[orthogonal polynomials]], and thus different roots where the integrand is evaluated.

 

It is also possible to use Clenshaw–Curtis quadrature to compute integrals of the form <math>\int_0^\infty f(x) dx</math> and <math>\int_{-\infty}^\infty f(x) dx</math>, using a coordinate-remapping technique.<ref name=Boyd87>John P. Boyd, "Exponentially convergent Fourier–Chebshev <nowiki>[</nowiki>''sic''<nowiki>]</nowiki> quadrature schemes on bounded and infinite intervals," ''J. Scientific Computing'' '''2''' (2), p. 99-109 (1987).</ref> High accuracy, even exponential convergence for smooth integrands, can be retained as long as <math>f(x)</math> decays sufficiently quickly as |''x''| approaches infinity.

 

Another coordinate-remapping approach was suggested for integrals of the form <math>\int_0^\infty e^{-x} g(x) dx</math>, in which case one can use the transformation <math>x = -\ln[(1 + \cos\theta)/2]</math> to transform the integral into the form <math>\int_0^\pi f(\cos\theta)\sin\theta \,d\theta</math> where <math>f(u) = g(-\ln[(1+u)/2])/2</math>, at which point one can proceed identically to Clenshaw–Curtis quadrature for ''f'' as above.<ref name=Basu77>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).</ref> Because of the endpoint singularities in this coordinate remapping, however, one uses Fejér's first quadrature rule [which does not evaluate ''f''(&minus;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. Instead, one normally precomputes quadrature weights <math>w_n</math> (for ''n'' from 0 to ''N''/2, assuming that ''N'' is even) so that

 

The weights <math>w_n</math> are positive and their sum is equal to one.<ref>J. P. Imhof, "On the Method for Numerical Integration of Clenshaw and Curtis", ''Numerische Mathematik'' '''5''', p. 138-141 (1963).</ref>

 

* [[Euler–Maclaurin formula]]

* [[Gauss–Kronrod quadrature formula]]