%% LEGPTS_DEMO - Computing Gauss quadrature nodes and weights with the chebfun system.
%
% This demo is designed to be used with the chebfun system
% http://www.comlab.ox.ac.uk/chebfun
% and compiled with 'publish'.
%
% e.g. web(publish('legpts_demo'))
%
% -----------------------------------------------
%%
% The chebfun system can efficiently compute Legendre points (the roots of
% Legendre polynomials), and the corresponding weights for Gauss-Legendre
% quadrature in two ways.
%%
% The first, typically used when N (the number of points to be computed) is
% small uses a standard tridiagonal eigenvalue problem of Golub and Welsch [1].
N = 32;
x = legpts(N);
%%
% When N is large (say greater than 128), the fast algorithm [2] developed by
% Glaser et al is used. Although a little less accurate than 'GW', this
% algorithm has a complexity cost of only O(N), rather than O(N^3) for the
% eigenvalue problem. (Type "help legpts" for more information).
N = 256;
x = legpts(N);
%%
% To see the difference in speeds, we could type
N = 128;
tic
[x w] = legpts(N,'GW');
toc
tic
[x w] = legpts(N,'Fast');
toc
%%
% Here the double output "[x w] = ..." will return also the weights needed
% for Gauss quadrature.
%
% The difference in times becomes even more dramtic when N is made larger.
N = 512;
tic
[x w] = legpts(N,'GW');
toc
tic
[x w] = legpts(N,'Fast');
toc
%%
% The fast algorithm is still pretty quick, even when N is very large.
N = 100000;
tic
[x w] = legpts(N,'Fast');
toc
%%
% We could now compute the Gauss quadrature approximation to the integral
% of the following recursively defined function
% Eqn (1.0)
f = @(x) sin(pi*x);
s = f;
for j = 1:10;
f = @(x) (3/4)*(1 - 2*f(x).^4);
s = @(x) s(x) + f(x);
end
f = @(x) s(x);
%%
% over [-1,1] by
IGQ = w*f(x)
toc
%%
% Of course, we could use the chebfun system itself to compute the same
% integral, using Clenshaw-Curtis rather than Gauss quadrature [3].
tic
g = chebfun(f)
ICC = sum(g)
toc
%%
% Chebfun allows us to plot the function easily
plot(g)
title('The recursive function in Eqn (1.0)');
%%
% We can check the accuracy of these computations by comparing their values
abs(IGQ-ICC)
%%
% or by plotting convergence as N is increased
NN = 32:32:2*length(g)/3;
errGQ = []; errCC = [];
for N = NN
[x w] = legpts(N);
errGQ = [errGQ ; abs(w*f(x)-ICC)];
errCC = [errCC ; abs(sum(chebfun(f,N))-ICC)];
end
semilogy(NN, errGQ, 'b', NN, errCC, 'r');
legend('Gauss Quadrature','Clenshaw-Curtis via Chebfun')
title('Convergence of quadratures for recursive function Eqn (1.0)');
%%
% See http://www.comlab.ox.ac.uk/chebfun for chebfun information.
%
% Copyright 2009 by The Chebfun Team.
%
% References:
% [1] G. H. Golub and J. A. Welsch, "Calculation of Gauss quadrature
% rules", Math. Comp. 23:221-230, 1969,
% [2] A. Glaser, X. Liu and V. Rokhlin, "A fast algorithm for the
% calculation of the roots of special functions", SIAM Journal
% on Scientific Computing, 29(4):1420-1438, 2007.
% [3] L. N. Trefethen, "Is Gauss quadrature better than Clenshaw-Curtis?",
% SIAM Review, 50(1):67–87, 2008.
