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

barymat(M,N,map,w)
function B = barymat(M,N,map,w) %BARYMAT(Y,X,W)
% BARYMAT  Barycentric Interpolation Matrix
%  BARYMAT(Y,X,W), where Y is a vector of length(M) and X and W are vectors
%  of length(N), returns the M*N matrix which interpolates data from the grid
%  X to the grid Y using barycentric weights W. If W is not supplied it is
%  assumed to be the Chebyshev weights W(j) = (-1)^j, W([1 N]) = 0.5*W([1 N]).
%
%  The points Y are assumed to lie within the interval [X(1) X(end)].
%
%  BARYMAT(M,N) returns the M*N matrix which interpolates a set of values
%  on an N-point Chebyshev grid to values on an M-point Chebyshev grid.
%
%  BARYMAT(M,N,MAP) is the same as the above, but for the mapped Chebyshev
%  grid defined by the map structure or function handle MAP.
%
%  BARYMAT(M,N,W) or BARYMAT(M,N,MAP,W) is similar, but uses the barycentric
%  weights given in W, which should be a vector of length N.

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

if length(M) == 1 && length(M) == length(N)
% Size of matrix is given.
if M == N    % Nothing to do here.
B = eye(M);
return
else
if nargin < 4, w = []; end
if nargin < 3,
map = [];          % No map given. Default to Chebyshev.
elseif isstruct(map)
map = map.for;     % Extract the forward map from map structure.
elseif ~isempty(map) && isnumeric(map)
w = map; map = []; % map is actually the barycentric weights.
end
end

x = chebpts(N);            % The original grid
y = chebpts(M);            % The up/down-sampled grid
if ~isempty(map), x = map(x); y = map(y); end % Map the points

elseif isempty(M) || isempty(N)
% Not much to do here.
B = []; return,
else
% Grid points are given.

x = N; N = length(x);
y = M; M = length(y);
if nargin == 3, w = map; else w = []; end
if M == N && all(x==y)     % Nothing to do here
B = eye(N);
return
end

end

if isempty(w) % Default to the Chebyshev barycentric weights
w = ones(N,1); w(2:2:end) = -1; w([1 N]) = 0.5*w([1 N]);
end

% Construct the matrix
if M >= 500 && N >= 1000         % <-- Experimentally determined.
% Testing shows BSXFUN is faster in this regime
B = bsxfun(@minus,y,x');     % Repmat(Y-X')
B = bsxfun(@rdivide,w',B);  % w(k)/(y(j)-x(k))
c = 1./sum(B,2);             % Normalisation ('denom' in bary-speak).
B = bsxfun(@times,B,c);
else
% Else use for loops
B = bsxfun(@minus,y,x');     % Repmat(Y-X')
for k = 1:N
B(:,k) = w(k)./B(:,k);   % w(k)/(y(j)-x(k))
end
c = 1./sum(B,2);             % Normalisation ('denom' in bary-speak).
for j = 1:M
B(j,:) = B(j,:)*c(j);
end
end

% Where points conincide there will be division by zeros (as with bary.m).
% Replace these entries with the identity.
B(isnan(B)) = 1;