Code covered by the BSD License

# Interpolation Utilities

### Joe Henning (view profile)

22 May 2012 (Updated )

A variety of interpolation utilities

schwerner(x, y, xi, p, q)
```function [yi, u, ypi, yppi] = schwerner(x, y, xi, p, q)

% SCHWERNER Rational interpolation using the Schneider-Werner Method
%    SCHWERNER(X,Y,XI,P,Q) interpolates to find YI, the values of the
%    underlying function Y at the points in the array XI, using
%    the Schneider Werner method and the second form (or true form)
%    of the barycentric interpolation formula.  X and Y must be
%    vectors of length N.
%
%    P and Q specify the degrees of the numerator and denominator
%    of the rational function.  P <= Q and P+Q=N-1.  If not specified,
%    values for P and Q are calculated.
%
%    [YI,U] = SCHWERNER() also returns the barycentric interpolation
%    weights U.
%
%    [YI,U,YPI,YPPI] = SCHWERNER() also returns the interpolated
%    derivative YPI and the interpolated second derivative YPPI.

% Joe Henning - Fall 2011

% Recent developments in barycentric rational interpolation
% Jean-Paul Berrut, Richard Baltensperger and Hans D. Mittelmann
% Trends and Applications in Constructive Approximation
% International Series of Numerical Mathematics Vol. 1??
% 2005 Birkhauser Verlag Basel (ISBN 3-7643-7124-2)

x = x(:).';
y = y(:).';
n = length(x);

if (nargin < 4)
p = floor((n-1)/2.0);
q = ceil((n-1)/2.0);
elseif (nargin < 5)
q = (n-1) - p;
end

if (p > q)
fprintf('??? Bad p, q input to schwerner ==> p <= q\n');
yi = [];
u = [];
ypi = [];
yppi = [];
return
end

if ((p+q) ~= (n-1))
fprintf('??? Bad p, q input to schwerner ==> p+q=n-1\n');
yi = [];
u = [];
ypi = [];
yppi = [];
return
end

if (p < 0)
fprintf('??? Bad p input to schwerner ==> p>=0\n');
yi = [];
u = [];
ypi = [];
yppi = [];
return
end

% construct A
A = ones(1,n);
for i = 1:p-1
A = [A; x.^i];
end
A = [A; y];
for i = 1:q-1
A = [A; y.*x.^i];
end

% calculate the kernel of A
%u = null(A);
u = kernel(A);

% check the dimension of the kernel
if (size(u,2) > 1)
fprintf('The kernel of A has a dimension larger than 1; recalling with p=%d-1=%d\n',p,p-1);
[yi, u, ypi, yppi] = schwerner(x, y, xi, p-1);
return
end

% remove unattainable points
[N,M] = size(u);
tol = reshape(eps*max([ones(size(u(:))) abs(u(:))],[],2),N,M);
i = find(abs(u) < tol);
if ~isempty(i)
for k = 1:length(i)
fprintf('(%f,%f) is an unattainable point; it will be eliminated\n',x(i(k)),y(i(k)));
end
k = find(abs(u) >= tol);
[yi, u, ypi, yppi] = schwerner(x(k), y(k), xi);
return
end

for i = 1:length(xi)
% Find the right place in the table by means of a bisection.
klo = 1;
khi = n;
while (khi-klo > 1)
k = fix((khi+klo)/2.0);
if (x(k) > xi(i))
khi = k;
else
klo = k;
end
end

h = x(khi) - x(klo);
if (h == 0.0)
fprintf('??? Bad x input to schwerner ==> x values must be distinct\n');
yi(i) = NaN;
ypi(i) = NaN;
yppi(i) = NaN;
continue;
end

isiny = 0;
for k = 1:n
if (xi(i) == x(k))
yi(i) = y(k);
num = 0;
for j = 1:n
if (j ~= k)
num = num + u(j)*(y(k)/(x(k)-x(j)) + y(j)/(x(j)-x(k)));
end
end
ypi(i) = -num/u(k);
num = 0;
for j = 1:n
if (j ~= k)
num = num + u(j)*(ypi(i)-(y(k)/(x(k)-x(j)) + y(j)/(x(j)-x(k))))/(x(k)-x(j));
end
end
yppi(i) = -2*num/u(k);
isiny = 1;
break
end
end

if (isiny)
continue
end

% Evaluate polynomial
num = 0;
den = 0;
for k = 1:n
num = num + y(k)*u(k)/(xi(i)-x(k));
den = den + u(k)/(xi(i)-x(k));
end
yi(i) = num/den;

num = 0;
den = 0;
for k = 1:n
term = yi(i)/(xi(i)-x(k)) + y(k)/(x(k)-xi(i));
num = num + term*u(k)/(xi(i)-x(k));
den = den + u(k)/(xi(i)-x(k));
end
ypi(i) = num/den;

num = 0;
den = 0;
for k = 1:n
term = (ypi(i)-(yi(i)/(xi(i)-x(k)) + y(k)/(x(k)-xi(i))))/(xi(i)-x(k));
num = num + term*u(k)/(xi(i)-x(k));
den = den + u(k)/(xi(i)-x(k));
end
yppi(i) = 2*num/den;
end

function Y = kernel(X)
% Calculate the kernel or null space of matrix X
[m,n] = size(X)
[U,S,V] = svd(X);
s = [];
for i = 1:m
s = [s; S(i,i)];
end
%tol = eps*max([ones(size(s)) abs(s)],[],2);
tol = eps*max([max(m,n) max(abs(s))]);
r = sum(s > tol);
Y = V(:,r+1:n);
```