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

lagpts(n,int,meth)
function [x w v] = lagpts(n,int,meth)
%LAGPTS  Laguerre points and Gauss-Laguerre Quadrature Weights.
%  LAGPTS(N) returns N Laguerre points X in (0,inf).
%
%  [X,W] = LAGPTS(N) returns also a row vector W of weights for Gauss-Laguerre
%  quadrature. [X,W,V] = LAGPTS(N) returns in addition a column vector V
%  of the barycentric weights corresponding to X.
%
%  LAGPTS(N,D) scales the nodes and weights for the semi-infinite domain D.
%  D can be either a domain object or a vector with two components.
%
%  [X,W] = LAGPTS(N,METHOD) allows the user to select which method to use.
%  METHOD = 'GW' will use the traditional Golub-Welsch eigenvalue method,
%  which is best for when N is small. METHOD = 'FAST' will use the
%  Glaser-Liu-Rokhlin fast algorithm, which is much faster for large N.
%  By default LEGPTS uses 'GW' when N < 128.
%

%  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, March 2010 - 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.

method = 'default';
interval = [0 inf];

if n < 0
error('CHEBFUN:lagpts: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 ~(strcmpi(method,'default') || strcmpi(method,'GW') || strcmpi(method,'fast'))
error('CHEBFUN:lagpts:inputs',['Unrecognised input string.', method]);
end
if isa(interval,'domain')
interval = interval.endsandbreaks;
end
if numel(interval) > 2,
warning('CHEBFUN:lagpts:domain',...
'Piecewise intervals not supported and will be ignored.');
interval = interval([1 end]);
end
end

if ~sum(isinf(interval)) == 1
error('CHEBFUN:lagpts:inf','lagpts only supports semi-infinite domains.');
end

% decide to use GW or FAST
if (n < 128 || strcmpi(method,'GW')) && ~strcmpi(method,'fast') % GW, see [1]
alpha = 2*(1:n)-1;  beta = 1:n-1;     % 3-term recurrence coeffs
T = diag(beta,1) + diag(alpha) + diag(beta,-1);  % Jacobi matrix
[V,D] = eig(T);                       % eigenvalue decomposition
[x,indx] = sort(diag(D));             % Laguerre points
w = V(1,indx).^2;                     % Quadrature weights
v = sqrt(x).*abs(V(1,indx)).';        % Barycentric weights
v = v./max(v); v(2:2:n) = -v(2:2:n);
else                                                            % Fast, see [2]
[x ders] = alg0_Lag(n);               % Nodes and L_n'(x)
w = exp(-x)./(x.*ders.^2); w = w';    % Quadrature weights
v = exp(-x/2)./ders;                  % Barycentric weights
v = -v./max(abs(v));
end
w = (1/sum(w))*w;                        % Normalise so that sum(w) = 1

% Nonstandard interval
if ~all(interval == [0 inf])
a = interval(1); b = interval(2);
if isinf(b)
x = x + a;
w = w*exp(-a);
else
x = - x + b;
w = w*exp(b);
end
end

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

function [x ders] = alg0_Lag(n)
ders = zeros(n,1);
xs = 1/(2*n+1);
n1 = 20;
n1 = min(n1, n);
x = zeros(n,1);
for k = 1:n1
[xs ders(k)] = alg3_Lag(n,xs);
x(k) = xs;
xs = 1.1*xs;
end
[x ders] = alg1_Lag(x,ders,n,n1);

% --------------------------- UNSCALED VERSION ------------------------------

function [roots ders] = alg1_Lag2(roots,ders,n,n1)
m = 30;
u = zeros(1,m+1); up = zeros(1,m+1);
for j = n1:n-1
x = roots(j);
h = rk2_Lag(pi/2,-pi/2,x,n) - x;

r = x.*(n+.5 -.25*x);  p = x.^2;
u(1:2) = [0; ders(j)];
u(3) = (-x*u(2)/2-r*u(1))/p;
u(4) = (-x*u(3)-(1+r)*u(2)/6-(n+.5-.5*x)*u(1))/p;
up(1:3) = [u(2) ; 2*u(3) ; 3*u(4)];

for k = 2:m-2
%         u(k+3)=(-x*(2*k+1)*u(k+2)-(k*k+r)*u(k+1)-k*(n+.5-.5*x)*u(k)+.25*k*(k-1)*u(k-1))/p;
u(k+3)=(-x*(2*k+1)*(k+1)*u(k+2)-(k*k+r)*u(k+1)-(n+.5-.5*x)*u(k)+.25*u(k-1))/(p*(k+2)*(k+1));
up(k+2) = (k+2)*u(k+3);
end

hh = [1;cumprod(h+zeros(m,1))];

for l = 1:5
%         hh = [1;cumprod(h*ones(m,1))];
hh = [1;cumprod(h+zeros(m,1))];
h = h - (u*hh)./(up*hh);
end

roots(j+1) = roots(j) + h;
ders(j+1) = up*[1;cumprod(h*ones(m,1))];

end

% --------------------------- SCALED VERSION ------------------------------

function [roots ders] = alg1_Lag(roots,ders,n,n1)
m = 30;
% Storage
hh1 = ones(m+1,1); zz = zeros(m,1); u = zeros(1,m+1); up = zeros(1,m+1);
x = roots(n1);
for j = n1:n-1
h = rk2_Lag(pi/2,-pi/2,x,n) - x;

M = 1/h; M2 = M^2; M3 = M^3; M4 = M^4;

r = x*(n + .5 - .25*x);  p = x^2;
u(1:2) = [0; ders(j)/M];
u(3) = -.5*u(2)/(M*x)-(n + .5 - .25*x)*u(1)/(x*M^2);
u(4) = -u(3)/(M*x)+(-(1+r)*u(2)/6/M^2-(n+.5-.5*x)*u(1)/M^3)/p;
up(1:3) = [u(2) ; 2*u(3)*M ; 3*u(4)*M];

for k = 2:m-2
%         u(k+3) = (-x*(2*k+1)*u(k+2)-(k*k+r)*u(k+1)-k*(n+.5-.5*x)*u(k)+.25*k*(k-1)*u(k-1))/p;
u(k+3) = (-x*(2*k+1)*(k+1)*u(k+2)/M-(k*k+r)*u(k+1)/M2-(n+.5-.5*x)*u(k)/M3+.25*u(k-1)/M4)/(p*(k+2)*(k+1));
up(k+2) = (k+2)*u(k+3)*M;
end
up(m+1) = 0;

% Flip for more accuracy in inner product calculation.
u = u(m+1:-1:1);  up = up(m+1:-1:1);

% Newton iteration
hh = hh1; hh(end) = M;    step = inf;  l = 0;

if M == 1
Mhzz = (M*h)+zz;
hh = [M;cumprod(Mhzz)];
hh = hh(end:-1:1);
end

while (abs(step) > eps) && (l < 10)
l = l + 1;
step = (u*hh)/(up*hh);
h = h - step;
Mhzz = (M*h)+zz;
hh = [M;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) = up*hh;

end

% --------------------------- SCALED VERSION ------------------------------

function [roots ders] = alg4_Lag(roots,ders,n,n1)
m = 30;
% Storage
hh1 = ones(m+1,1); zz = zeros(m,1); u = zeros(1,m+1); up = zeros(1,m+1);
x = roots(n1);
for j = n1:n-1
h = rk2_Lag(pi/2,-pi/2,x,n) - x;

M = 1/h; M2 = M^2; M3 = M^3; M4 = M^4;

r = x*(n + .5 - .25*x);  p = x^2;
u(1:2) = [0; ders(j)/M];
u(3) = -.5*u(2)/(M*x)-(n + .5 - .25*x)*u(1)/(x*M^2);
u(4) = -u(3)/(M*x)+(-(1+r)*u(2)/6/M^2-(n+.5-.5*x)*u(1)/M^3)/p;
up(1:3) = [u(2) ; 2*u(3)*M ; 3*u(4)*M];

for k = 2:m-2
%         u(k+3) = (-x*(2*k+1)*u(k+2)-(k*k+r)*u(k+1)-k*(n+.5-.5*x)*u(k)+.25*k*(k-1)*u(k-1))/p;
u(k+3) = (-x*(2*k+1)*(k+1)*u(k+2)/M-(k*k+r)*u(k+1)/M2-(n+.5-.5*x)*u(k)/M3+.25*u(k-1)/M4)/(p*(k+2)*(k+1));
up(k+2) = (k+2)*u(k+3)*M;
end
up(m+1) = 0;

% Flip for more accuracy in inner product calculation.
u = u(m+1:-1:1);  up = up(m+1:-1:1);

% Newton iteration
hh = hh1; hh(end) = M;    step = inf;  l = 0;

if M == 1
Mhzz = (M*h)+zz;
hh = [M;cumprod(Mhzz)];
hh = hh(end:-1:1);
end

while (abs(step) > eps) && (l < 10)
l = l + 1;
step = (u*hh)/(up*hh);
h = h - step;
Mhzz = (M*h)+zz;
hh = [M;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) = up*hh;

end

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

function [x1 d1] = alg3_Lag(n,xs)
[u up] = eval_Lag(xs,n);
theta = atan(sqrt(xs/(n+.5-.25*xs))*up/u);
x1 = rk2_Lag(theta,-pi/2,xs,n);

% Newton iteration
step = inf;  l = 0;
while (abs(step) > eps || abs(u) > eps) && (l < 200)
l = l + 1;
[u up] = eval_Lag(x1,n);
step = u/up;
x1 = x1 - step;
end

[ignored d1] = eval_Lag(x1,n);

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

function [L Lp] = eval_Lag(x,n)
Lm2 = 0; Lm1 = exp(-x/2); Lpm2 = 0; Lpm1 = 0;
for k = 0:n-1
L = ((2*k+1-x).*Lm1-k*Lm2)/(k+1);
Lp = ((2*k+1-x).*Lpm1-Lm1-k*Lpm2)/(k+1);
Lm2 = Lm1; Lm1 = L; Lpm2 = Lpm1; Lpm1 = Lp;
end

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

function x = rk2_Lag(t,tn,x,n)
m = 10; h = (tn-t)/m;
for j = 1:m
f1 = (n+.5-.25*x);
k1 = -h/(sqrt(f1/x)+.25*(1/x-.25/f1)*sin(2*t));
t = t+h;  x = x+k1;   f1 = (n+.5-.25*x);
k2 = -h/(sqrt(f1/x)+.25*(1/x-.25/f1)*sin(2*t));
x = x+.5*(k2-k1);
end