Code covered by the BSD License

# Chebfun V4

### Chebfun Team (view profile)

30 Apr 2009 (Updated )

Numerical computation with functions instead of numbers.

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

legpts(n,int,meth)
function [x w v] = legpts(n,int,meth)
%LEGPTS  Legendre points and Gauss Quadrature Weights.
%  LEGPTS(N) returns N Legendre points X in (-1,1).
%
%  [X,W] = LEGPTS(N) returns also a row vector W of weights for Gauss quadrature.
%
%  LEGPTS(N,D) scales the nodes and weights for the domain D. D can be
%  either a domain object or a vector with two components. If the interval
%  is infinite, the map is chosen to be the default 'unbounded map' with
%  mappref('parinf') = [1 0] and mappref('adaptinf') = 0.
%
%  [X,W,V] = LEGPTS(N) returns additionally a column vector V of weights in
%  the barycentric formula corresponding to the points X.
%
%  [X,W] = LEGPTS(N,METHOD) allows the user to select which method to use.
%    METHOD = 'FASTSMALL' uses the recurrence relation for the Legendre
%       polynomials and their derivatives to perform Newton iteration
%       on the WKB approximation to the roots.
%    METHOD = 'FAST' uses the Glaser-Liu-Rokhlin fast algorithm, which
%       is much faster for large N.
%    METHOD = 'GW' will use the traditional Golub-Welsch eigenvalue method,
%       which is maintained mostly for historical reasons.
%       By default LEGPTS uses 'FASTSMALL' when N<=256 and FAST when N>256
%

%  Copyright 2011 by The University of Oxford and The Chebfun Developers.
%  See http://www.maths.ox.ac.uk/chebfun/ for Chebfun information.

%  'GW' by Nick Trefethen, March 2009 - algorithm adapted from [1].
%  'FAST' by Nick Hale, April 2009 - algorithm adapted from [2].
%
%  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.

% Defaults
interval = [-1,1];
method = 'default';

if n < 0
error('CHEBFUN:legpts:n','First input should be a positive number.');
end

% Return empty vector if n == 0
if n == 0
x = []; w = []; v = []; return
end

% Check the inputs
if nargin > 1
if nargin == 3
interval = int; method = meth;
elseif nargin == 2
if ischar(int), method = int; else interval = int; end
end
if ~any(strcmpi(method,{'default','GW','fast','fastsmall'}))
error('CHEBFUN:legpts:inputs',['Unrecognised input string.', method]);
end
if isa(interval,'domain')
interval = interval.endsandbreaks;
end
if numel(interval) > 2,
warning('CHEBFUN:legpts:domain',...
'Piecewise intervals not supported and will be ignored.');
interval = interval([1 end]);
end
end

% Decide to use GW or FAST
if n == 1
% Trivial case when N = 1
x = 0; w = 2; v = 1;
elseif strcmpi(method,'GW')
% GW, see [1]
beta = .5./sqrt(1-(2*(1:n-1)).^(-2)); % 3-term recurrence coeffs
T = diag(beta,1) + diag(beta,-1);     % Jacobi matrix
[V,D] = eig(T);                       % Eigenvalue decomposition
x = diag(D); [x,i] = sort(x);         % Legendre points
w = 2*V(1,i).^2;                      % Quadrature weights
v = sqrt(1-x.^2).*abs(V(1,i))';       % Barycentric weights
v = v./max(v);

% Enforce symmetry
ii = 1:floor(n/2);  x = x(ii);  w = w(ii);
vmid = v(floor(n/2)+1); v = v(ii);
if mod(n,2)
x = [x ; 0 ; -x(end:-1:1)];  w = [w  2-sum(2*w) w(end:-1:1)];
v = [v ; vmid ; v(end:-1:1)];
else
x = [x ; -x(end:-1:1)];      w = [w w(end:-1:1)];
v = [v ; v(end:-1:1)];
end
v(2:2:n) = -v(2:2:end);
elseif (n < 256 && ~strcmpi(method,'fast')) || strcmpi(method,'fastsmall')
% Fastsmall
[x ders] = fastsmall(n);              % Nodes and P_n'(x)
w = 2./((1-x.^2).*ders.^2)';          % Quadrature weights
v = 1./ders; v = v./max(abs(v));      % Barycentric weights
if ~mod(n,2), ii = (floor(n/2)+1):n; v(ii) = -v(ii);   end
else
% Fast, see [2]
[x ders] = alg0_Leg(n);               % Nodes and P_n'(x)
w = 2./((1-x.^2).*ders.^2)';          % Quadrature weights
v = 1./ders; v = v./max(abs(v));      % Barycentric weights
if ~mod(n,2), ii = (floor(n/2)+1):n;  v(ii) = -v(ii);   end
end

% Normalise so that sum(w) = 2
w = (2/sum(w))*w;

% Rescale to arbitrary interval
if ~all(interval == [-1 1])
if ~any(isinf(interval))
% Finite interval
dab = diff(interval);
x = (x+1)/2*dab + interval(1);
w = dab*w/2;
else
% Infinite interval
m = maps(fun,{'unbounded'},interval); % use default map
if nargout > 1
w = w.*m.der(x.');
end
x = m.for(x);
x([1 end]) = interval([1 end]);
end
end

% -------------------- Routines for FAST_SMALL algorithm ------------------

function [x PP] = fastsmall(n)

% Asymptotic formula (WKB) - only positive x.
if mod(n,2), s = 1; else s = 0; end
k = (n+s)/2:-1:1; theta = pi*(4*k-1)/(4*n+2);
x = (1-(n-1)/(8*n^3)-1/(384*n^4)*(39-28./sin(theta).^2)).*cos(theta);

% Initialise
Pm2 = 1; Pm1 = x;  PPm2 = 0; PPm1 = 1;
dx = inf; l = 0;

% Loop until convergence
while norm(dx,inf) > eps && l < 10
l = l + 1;
for k = 1:n-1,
P = ((2*k+1)*Pm1.*x-k*Pm2)/(k+1);           Pm2 = Pm1; Pm1 = P;
PP = ((2*k+1)*(Pm2+x.*PPm1)-k*PPm2)/(k+1);  PPm2 = PPm1; PPm1 = PP;
end
dx = -P./PP; x = x + dx;
Pm2 = 1; Pm1 = x; PPm2 = 0; PPm1 = 1;
end

% Once more for derivatives
for k = 1:n-1,
P = ((2*k+1)*Pm1.*x-k*Pm2)/(k+1);           Pm2 = Pm1; Pm1 = P;
PP = ((2*k+1)*(Pm2+x.*PPm1)-k*PPm2)/(k+1);  PPm2 = PPm1; PPm1 = PP;
end

% Reflect for negative values
x = [-x(end:-1:1+s) x].';
PP = [PP(end:-1:1+s) PP].';

% -------------------- Routines for FAST algorithm ------------------------

function [roots ders] = alg0_Leg(n) % Driver for 'Fast'.

% Compute coefficients of P_m(0), m = 0,..,N via recurrence relation.
Pm2 = 0; Pm1 = 1;
for k = 0:n-1, P = -k*Pm2/(k+1); Pm2 = Pm1; Pm1 = P; end

% Get the first roots and derivative values to initialise.
roots = zeros(n,1); ders = zeros(n,1);                      % Allocate storage
if mod(n,2)                                                 % n is odd
roots((n-1)/2) = 0;                                     % Zero is a root
ders((n+1)/2) = n*Pm2;                                  % P'(0)
else                                                        % n is even
[roots(n/2+1) ders(n/2+1)] = alg2_Leg(P,n);             % Find first root
end

[roots ders] = alg1_Leg(roots,ders);          % Other roots and derivatives

% -------------------------------------------------------------------------

function [roots ders] = alg1_Leg(roots,ders)  % Main algorithm for 'Fast'
n = length(roots);
if mod(n,2), N = (n-1)/2; s = 1; else N = n/2; s = 0; end

% Approximate roots via asymptotic formula.
k = (n-2+s)/2:-1:1; theta = pi*(4*k-1)/(4*n+2);
roots(((n+4-s)/2):end) = (1-(n-1)/(8*n^3)-1/(384*n^4)*(39-28./sin(theta).^2)).*cos(theta);
x = roots(N+1);

% Number of terms in Taylor expansion.
m = 30;

% Storage
hh1 = ones(m+1,1); zz = zeros(m,1); u = zeros(1,m+1); up = zeros(1,m+1);

% Loop over all the roots we want to find (using symmetry).
for j = N+1:n-1
% Distance to initial approx for next root (from asymptotic foruma).
h = roots(j+1) - x;

% Recurrence Taylor coefficients (scaled & incl factorial terms).
M = 1/h;                           % Scaling
c1 = 2*x/M; c2 = 1./(1-x^2);       % Some constants
% Note, terms are flipped for more accuracy in inner product calculation.
u([m+1 m]) = [0 ders(j)/M];  up(m+1) = u(m);
for k = 0:m-2
up(m-k) = (c1*(k+1)*u(m-k)+(k-n*(n+1)/(k+1))*u(m-k+1)/M^2)*c2;
u(m-(k+1)) = up(m-k)/(k+2);
end
up(1) = 0;

% Newton iteration
hh = hh1; step = inf;  l = 0;
while (abs(step) > eps) && (l < 10)
l = l + 1;
step = (u*hh)/(up*hh)/M;
h = h - step;
Mhzz = (M*h)+zz;
hh = [1;cumprod(Mhzz)];     % Powers of h (This is the fastest way!)
hh = hh(end:-1:1);          % Flip for more accuracy in inner product
end

% Update
x = x + h;
roots(j+1) = x;
ders(j+1) = M*(up*hh);

end

% Nodes are symmetric.
roots(1:N+s) = -roots(n:-1:N+1);
ders(1:N+s) = ders(n:-1:N+1);

% -------------------------------------------------------------------------

function [x1 d1] = alg2_Leg(Pn0,n) % Find the first root (note P_n'(0)==0)
% Approximate first root via asymptotic formula
k = ceil(n/2); theta = pi*(4*k-1)/(4*n+2);
x1 = (1-(n-1)/(8*n^3)-1/(384*n^4)*(39-28./sin(theta).^2)).*cos(theta);

m = 30; % Number of terms in Taylor expansion.

% Recurrence Taylor coefficients (scaled & incl factorial terms).
M = 1/x1; % Scaling
zz = zeros(m,1); u = [Pn0 zeros(1,m)]; up = zeros(1,m+1); % Allocate storage
for k = 0:2:m-2
up(k+2) = (k-n*(n+1)/(k+1))*u(k+1)/M^2;
u(k+3) = up(k+2)/(k+2);
end
% Flip for more accuracy in inner product calculation.
u = u(m+1:-1:1); up = up(m+1:-1:1);

% % Note, terms are flipped for more accuracy in inner product calculation.
% zz = zeros(m,1); u = [zeros(1,m) Pn0]; up = zeros(1,m+1); % Allocate storage
% for k = 0:2:m-2
%     up(m-k) = (k-n*(n+1)/(k+1))*u(m-k+1)/M^2;
%     u(m-(k+1)) = up(m-k)/(k+2);
% end

% Newton iteration
x1k = ones(m+1,1); step = inf; l = 0;
while (abs(step) > eps) && (l < 10)
l = l + 1;
step = (u*x1k)/(up*x1k)/M;
x1 = x1 - step;
x1k = [1;cumprod(M*x1+zz)]; % Powers of h (This is the fastest way!)
x1k = x1k(end:-1:1);        % Flip for more accuracy in inner product
end

% Get the derivative at this root, i.e. P'(x1).
d1 = M*(up*x1k);