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

jacpts(n,alpha,beta,varargin)
function [x,w,v] = jacpts(n,alpha,beta,varargin)
%JACPTS  Gauss-Jacobi Abscissae and Quadrature Weights.
%  X = JACPTS(N,ALPHA,BETA) returns the N roots of the degree N Jacobi
%       polynomial with parameters ALPHA and BETA (which must both be
%       greater than or equal -1)
%
%  [X,W] = JACPTS(N,ALPHA,BETA) returns also a row vector W of weights for
%
%  [X,W,V] = JACPTS(N,ALPHA,BETA) returns additionally a column vector V of
%       weights in the barycentric formula corresponding to the points X.
%
%  [X,W,V] = JACPTS(N,ALPHA,BETA,METHOD) allows choice in which method is used.
%       METHOD = 'GW' will use the traditional Golub-Welsch eigenvalue method,
%       which is best suited for when N is small. METHOD = 'FAST' will use
%       the Glaser-Liu-Rokhlin fast algorithm, which is much faster for large N.
%       By default JACPTS will use 'GW' when N < 128.
%
%  [X,W,V] = JACPTS(N,ALPHA,BETA,[A,B]) scales the nodes and weights for the
%       finite interval [A,B].
%
%  The cases ALPHA = BETA = -.5 and ALPHA = BETA = .5 correspond to
%  Gauss-Chebyshev nodes and quadrature, and are treated specially
%  (as a closed form of the nodes and weights is available).
%  ALPHA = BETA = 0 calls LEGPTS, which is a more efficient code.
%

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

%  'FAST' by Nick Hale, April 2009 - algorithm adapted from [1].
%
%  References:
%   [1] 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:jacpts:n','First input should be a positive number.');
end

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

if alpha <= -1 || beta <= -1,
error('CHEBFUN:jacpts:SizeAB','alpha and beta must be greater than -1')
end
a = alpha; b = beta;

% check inputs
if nargin > 3
if isa(varargin{1},'double') && length(varargin{1}) == 2
interval = varargin{1};
elseif isa(varargin{1},'domain')
interval = varargin{1}.ends;
elseif isa(varargin{1},'char')
method = varargin{1};
end
if length(varargin) == 2,
if isa(varargin{2},'double') && length(varargin{2}) == 2
interval = varargin{2};
elseif isa(varargin{1},'domain')
interval = varargin{2}.ends;
elseif isa(varargin{2},'char')
method = varargin{2};
end
end
end

if nargout > 1 && any(isinf(interval)) % infinite intervals not yet supported
% (How do we scale the weights?)
error('CHEBFUN:jacpts:infinterval', ...
'jacpts does not yet support infinite intervals');
end

% % Special cases
% Legendre: alpha = beta = 0
if ~(a || b) % The case alpha = beta = 0 is treated by legpts
[x w v] = legpts(n,varargin{:});
return
end
% Gauss-Chebyshev: alpha = beta = -.5
if a == -.5 && b == -.5
[x ignored v] = chebpts(n,interval,1);
w = repmat(pi/n,1,n);
return
end
% Gauss-Chebyshev2:  alpha = beta = .5
if a == .5 && b == .5
x = chebpts(n+2,2);     x = x(2:n+1);
w = pi/(n+1)*(1-x.^2);  w = w';
if nargout == 3, v = (1-x.^2);  v(2:2:end) = -v(2:2:end); end
[x w] = rescale(x,w,interval,alpha,beta);
return
end

% Fix n == 1 case for GW
if n == 1, method = 'fast'; end

% decide to use GW or FAST
if (n < 128 || strcmpi(method,'GW')) && ~strcmpi(method,'fast')
ab = a + b;
ii = (2:n-1)';
abi = 2*ii + ab;
aa = [(b - a)/(2 + ab) ; (b^2 - a^2)./((abi - 2).*abi) ; (b^2 - a^2)./((2*n - 2+ab).*(2*n+ab))];
bb = [2*sqrt( (1 + a)*(1 + b)/(ab + 3))/(ab + 2) ;
2*sqrt(ii.*(ii + a).*(ii + b).*(ii + ab)./(abi.^2 - 1))./abi];
TT = diag(bb,1) + diag(aa) + diag(bb,-1); % Jacobi matrix
[V x] = eig( TT );                        % Eigenvalue decomposition
x = diag(x);                              % Jacobi points
w = V(1,:).^2*( 2^(ab + 1) * gamma(a + 1) * gamma(b + 1) / gamma(2 + ab) ); % Quadrature weights
v = sqrt(1-x.^2).*abs(V(1,:))';           % Barycentric weights
v = v./max(v); v(2:2:n) = -v(2:2:end);
else   % Fast, see [2]
[x ders] = alg0_Jac(n,a,b);                % Nodes and P_n'(x)
w = 1./((1-x.^2).*ders.^2)';               % Quadrature weights
if a && b                                  % Get the right constant
C = 2^(a+b+1)*gamma(2+a)*gamma(2+b)/(gamma(2+a+b)*(a+1)*(b+1));
w = C*w/sum(w);
else
w = 2^(a+b+1)*w;
end
v = 1./ders; v = v./max(abs(v));           % Barycentric weights
if ~mod(n,2), v = -v;   end
end

[x w] = rescale(x,w,interval,alpha,beta);

function [x w] = rescale(x,w,interval,alpha,beta)
% rescale to arbitrary interval
if all(interval == [-1 1])
% Nothing to do
return
end

if ~any(isinf(interval))
% finite interval
c1 = .5*sum(interval);
c2 = .5*diff(interval);
w = c2^(alpha+beta+1)*w;
x = c1+c2*x;
else
% infinite interval (not yet supported)
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

function [roots ders] = alg0_Jac(n,a,b)
if abs(a)<=.5 && abs(b)<=.5 % use asymptotic formula
r = ceil(n/2); % Choose a root near the middle
C = (2*r+a-.5)*pi/(2*n+a+b+1);
T = C + 1/(2*n+a+b+1)^2*((.25-a^2)*cot(.5*C)-(.25-b^2)*tan(.5*C));
x1 = cos(T);

% Make accurate using Newton
[u up] = eval_Jac(x1,n,a,b);
while abs(u) > eps
x1 = x1 - u/up;
[u up] = eval_Jac(x1,n,a,b);
end

if n == 1, roots = x1; ders = up; return, end

[rootsl dersl] = alg1_Jac_as(n,x1,up,a,b,1); % Get roots to the left
if a ~= b
[rootsr dersr] = alg1_Jac_as(n,x1,up,a,b,0); % To the right
else
rootsr = -rootsl; % Use symmetry.
dersr = dersl;
if rootsl(1) > 0, rootsl(1) = []; dersl(1) = []; end
end

roots = [rootsl(end:-1:2) ; rootsr];
ders = [dersl(end:-1:2) ; dersr];
else

%    % Attempt 1: Starts at the left and works through
%     [x1 up] = alg3_Jac(n,1-eps,a,b) ;
%     d1 = up;
%     if n == 1, roots = x1; ders = up; return, end
%     if n > 1 % Get the 2nd root
%         % Initial guess
%         x2 = rk2_Jac(pi/2,-pi/2,x1,n,a,b);
%         % Make accurate using Newton
%         [u up] = eval_Jac(x2,n,a,b);
%         while abs(u) > eps
%             x2 = x2 - u/up;
%             [u up] = eval_Jac(x2,n,a,b);
%         end
%
%         if a ~= b
%             [roots ders] = alg1_Jac(n,x2,up,a,b,1);
%             roots = [x1 ; roots];
%             ders = [d1 ; ders];
%         else
%             [roots ders] = alg1_Jac(n,x2,up,a,b,0);
%             roots = [x1 ; roots ; -roots((end-mod(n,2):-1:1)) ; x1];
%             ders = [d1 ; ders ; ders((end-mod(n,2):-1:1)) ; d1];
%         end
%     end

%    % Attempt 2: Starts at the middle (more accurate)
[x1 up] = alg3_Jac(n,0,a,b);

if n == 1, roots = x1; ders = up; return, end

[rootsl dersl] = alg1_Jac(n,x1,up,a,b,1); % Get roots to the left
if a ~= b
[rootsr dersr] = alg1_Jac(n,x1,up,a,b,0); % To the right
else
rootsr = -rootsl(1+mod(n,2):end); % Use symmetry.
dersr = dersl(1+mod(n,2):end);
if rootsl(1) > 0, rootsl(1) = []; dersl(1) = []; end
end

roots = [rootsl(end:-1:2) ; rootsr];
ders = [dersl(end:-1:2) ; dersr];
end

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

function [roots ders] = alg1_Jac(n,x1,up,a,b,flag)
ab = a + b;
% if flag, N = n-2; else N = ceil(n/2)-2; end
if flag, sgn = -1; else sgn = 1; end
N = n-1;
roots = [x1 ; zeros(N,1)]; ders = [up ; zeros(N,1)];
m = 30; % number of terms in Taylor expansion
hh1 = ones(m+1,1); u = zeros(1,m+1); up = zeros(1,m+1);
for j = 1:N
x = roots(j);
h = rk2_Jac(sgn*pi/2,-sgn*pi/2,x,n,a,b) - x;

if abs(x+h) > 1, roots(j+1:end) = []; ders(j+1:end) = []; return, end

% scaling
M = 1/h;

% recurrence relation  (scaled)
u(1) = 0;   u(2) = ders(j)/M;  up(1) = u(2); up(m+1) = 0;
for k = 0:m-2
u(k+3) = ((x*(2*k+ab+2)+a-b)*u(k+2)/M + ...
(k*(k+ab+1)-n*(n+ab+1))*u(k+1)/M^2/(k+1))./((1-x.^2)*(k+2));
up(k+2) = (k+2)*u(k+3)*M;
end

% 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;
while (abs(step) > eps) && (l < 10)
l = l + 1;
step = (u*hh)/(up*hh);
h = h - step;
hh = [M;cumprod(M*h+zeros(m,1))]; % powers of h (This is the fastest way!)
hh = hh(end:-1:1);
end

if abs(h) < eps, roots(j+1:end) = []; ders(j+1:end) = []; return, end

% update
roots(j+1) = x + h;
ders(j+1) = up*hh;
end

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

function [roots ders] = alg1_Jac_as(n,x1,up,a,b,flag) % if |a|<=.5 && |b|<=.5 use asymptotic formula
ab = a + b;

% Approximate roots via asymptotic formula
nx1 = ceil(n/2);
if flag, r = (nx1+1):n; else r = (nx1-1):-1:1; end
C = (2*r+a-.5)*pi/(2*n+ab+1);
T = C + 1/(2*n+a+b+1)^2*((.25-a^2)*cot(.5*C)-(.25-b^2)*tan(.5*C));
roots = [x1 ; cos(T).']; ders = [up ; zeros(length(T),1)];

m = 30; % number of terms in Taylor expansion
hh1 = ones(m+1,1); u = zeros(1,m+1); up = zeros(1,m+1);
for j = 1:length(r)
x = roots(j); % previous root

% initial approx (via asymptotic foruma)
h = roots(j+1) - x;

% scaling
M = 1/h;

% recurrence relation for Jacobi polynomials (scaled)
cx = (1-x.^2);
u(1) = 0;   u(2) = ders(j)/M;  up(1) = u(2); up(m+1) = 0;
for k = 0:m-2
u(k+3) = ((x*(2*k+ab+2)+a-b)*u(k+2)/M + ...
(k*(k+ab+1)-n*(n+ab+1))*u(k+1)/M^2/(k+1))./(cx*(k+2));
up(k+2) = (k+2)*u(k+3)*M;
end

% 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;
while (abs(step) > eps) && (l < 10)
l = l + 1;
step = (u*hh)/(up*hh);
h = h - step;
hh = [M;cumprod(M*h+zeros(m,1))]; % powers of h (This is the fastest way!)
hh = hh(end:-1:1); % flip for more accuracy in inner product calculation
end

% update
roots(j+1) = x + h;
ders(j+1) = up*hh;
end

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

function [x1 d1] = alg3_Jac(n,xs,a,b)
[u up] = eval_Jac(xs,n,a,b);
theta = atan((1-xs.^2)*up/sqrt(n*(n+a+b+1)*(1-xs.^2))/u);
x1 = rk2_Jac(theta,-pi/2,xs,n,a,b);

for k = 1:10
[u up] = eval_Jac(x1,n,a,b);
x1 = x1 - u/up;
end

[u1 d1] = eval_Jac(x1,n,a,b);

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

function [P Pp] = eval_Jac(x,n,a,b)
ab = a + b;

P = .5*(a-b+(ab+2)*x);  Pm1 = 1;
Pp = .5*(ab+2);         Ppm1 = 0;

if n == 0;
P = Pm1; Pp = Ppm1;
end

for k = 1:n-1
A = 2*(k+1)*(k+ab+1)*(2*k+ab);
B = (2*k+ab+1)*(a^2-b^2);
C = prod(2*k+ab+(0:2)');
D = 2*(k+a)*(k+b)*(2*k+ab+2);

Pa1 = ( (B+C*x).*P - D*Pm1 ) / A;
Ppa1 = ( (B+C*x).*Pp + C*P - D*Ppm1 ) / A;

Pm1 = P; P = Pa1;  Ppm1 =  Pp; Pp = Ppa1;
end

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

function x = rk2_Jac(t,tn,x,n,a,b)
ab = a + b;
m = 10; h = (tn-t)/m;
for j = 1:m
f1 = (1-x.^2);
k1 = -4*h*f1./(4*sqrt(n*(n+ab+1)*f1) + (b-a-(ab+1)*x).*sin(2*t));
t = t+h;
f2 = (1-(x+k1).^2);
k2 = -4*h*f2./(4*sqrt(n*(n+ab+1)*f2) + (b-a-(ab+1)*(x+k1)).*sin(2*t));
x = x+.5*real(k1+k2);
end