Gauss and Clenshaw-Curtis quadrature

Nick Trefethen, September 2010

Suppose you have a function f on an interval:

x = chebfun('x');
f = @(x) x.*sin(2*exp(2*sin(2*exp(2*x))));
fc = chebfun(f);
LW = 'linewidth'; FS = 'fontsize'; MS = 'markersize';
figure, plot(fc,LW,1.2)
title('Function f',FS,16)

In Chebfun you would normally compute the integral like this:

format long
Ichebfun = sum(fc)

Ichebfun =
   0.336732834781727

Chebfun's method is Clenshaw-Curtis quadrature, i.e., the integration of the polynomial representing f by interpolation or piecewise interpolation in Chebyshev points. Here is the number of quadrature points:

Npts = length(fc)

Npts =
   651

If we wanted, we could also perform the integration by explicitly extracting the Clenshaw-Curtis nodes and weights, like this:

[s,w] = chebpts(Npts);
Iclenshawcurtis = w*f(s)

Iclenshawcurtis =
   0.336732834781728

Or we could try Gauss quadrature with the same number of points and weights.

[s,w] = legpts(Npts);
Igauss = w*f(s)

Igauss =
   0.336732834781725

Though this value of Npts is in the hundreds, Chebfun can handle values in the millions without difficulty. This is achieved by the algorithm of Glaser, Liu and Rokhlin [1] as implemented by Nick Hale. See Examples quad/GaussQuad.m.

Let's take a look at the accuracy as a function of Npts. Gauss quadrature converges geometrically, since f is analytic:

figure, tic, err = [];
NN = 10:10:500;
for Npts = NN
  [s,w] = legpts(Npts);
  Igauss = w*f(s);
  err = [err abs(Igauss-Ichebfun)];
end
semilogy(NN,err,'.-',LW,1,MS,16), grid on
ylim([1e-18 1])
xlabel('Npts',FS,12), ylabel('Error',FS,12)
title('Gauss quadrature convergence',FS,16), toc

Elapsed time is 0.566076 seconds.

Let's add another curve to the plot for Clenshaw-Curtis:

hold on, tic, err = [];
for Npts = NN
  [s,w] = chebpts(Npts);
  Iclenshawcurtis = w*f(s);
  err = [err abs(Iclenshawcurtis-Ichebfun)];
end
semilogy(NN,err,'.-r',LW,1,MS,16)
title('Gauss and Clenshaw-Curtis',FS,16)
legend('Gauss','Clenshaw-Curtis','location','southwest'), toc

Elapsed time is 0.506196 seconds.

Clenshaw-Curtis quadrature also converges geometrically for analytic functions. In some circumstances Gauss converges up to twice as fast as C-C, with respect to Npts, but as this example suggests, the two formulas are often closer than that. The computer time is often faster with C-C. For details of the cmoparison, see [2] and Chapter 19 of [3].

References:

[1] A. Glaser, X. Liu and V. Rokhlin, A fast algorithm for the calculation of the roots of special functions, SIAM J. Sci. Comp. 29 (2007), 1420-1438.

[2] L. N. Trefethen, Is Gauss quadrature better than Clenshaw-Curtis?, SIAM Review 50 (2008), 67-87.

[3] L. N. Trefethen, Approximation Theory and Approximation Practice, draft textbook available at http://www.maths.ox.ac.uk/chebfun/ATAP.