function interpMatrix = m20120703_03_polynomialInterpolationMatrix(xk, x, alpxk, alpx) % Returns the matrix for polynomial interpolation from the xk grid to the x % grid. % % This script is a modification of the polint.m routine from DMSuite, % by J.A.C. Weideman and S.C. Reddy, available here: % http://www.mathworks.com/matlabcentral/fileexchange/29 % or here: % http://dip.sun.ac.za/~weideman/research/differ.html % % The original polint.m routine performed the interpolation explicitly. % This modified function returns the interpolation matrix, for use in % implicit codes. % % Modification by Matt Landreman, MIT Plasma Science and Fusion Center, % July 3, 2012. % % Inputs % xk: Vector of x-coordinates of data (assumed distinct). % x: Vector of x-values where polynomial interpolant is to be evaluated. % alpxk: Vector of weight values sampled at the points xk. % alpx: Vector of weight values sampled at the points x. % % Output: % interpolation matrix % % The code implements the barycentric formula; see page 252 in % P. Henrici, Essentials of Numerical Analysis, Wiley, 1982. % (Note that if some fk > 1/eps, with eps the machine epsilon, % the value of eps in the code may have to be reduced.) % Note added May 2003: Except for certain nice node distributions % polynomial interpolation of high-degree is an ill-conditioned % problem. This code does not test for conditioning so use with % care. % J.A.C. Weideman, S.C. Reddy 1998. N = length(xk); M = length(x); if nargin == 2 alpx = 1; alpxk = 1; deWeightify = eye(N); elseif nargin > 2 deWeightify = diag(1./alpxk); else error('Unexpected number of inputs.') end x = x(:); % Make sure the data are column vectors xk = xk(:); alpxk = alpxk(:); alpx = alpx(:); L = logical(eye(N)); D = xk(:,ones(1,N))-xk(:,ones(1,N))'; % Compute the weights w(k) D(L) = ones(N,1); w = 1./prod(D)'; D = x(:,ones(1,N)) - xk(:,ones(1,M))'; % Compute quantities x-x(k) D = 1./(D+eps*(D==0)); % and their reciprocals. % Evaluate interpolant as matrix-vector products: interpMatrix = diag(alpx./(D*w))*D*diag(w)*deWeightify;