Code covered by the BSD License

# Interpolation Utilities

### Joe Henning (view profile)

22 May 2012 (Updated )

A variety of interpolation utilities

floaterhormann(x, y, xi, c)
```function [yi, u, ypi, yppi] = floaterhormann(x, y, xi, c)

% FLOATERHORMANN Rational interpolation using the Floater-Hormann Method
%    FLOATERHORMANN(X,Y,XI,C) interpolates to find YI, the values of the
%    underlying function Y at the points in the array XI, using
%    Floater-Hormann rational interpolation.  X and Y must be
%    vectors of length N.
%
%    C specifies the order of the interpolating polynomial
%    (0 <= C <= N).  The default is C = 0 (Berrut rational interpolation).
%
%    [YI,U] = FLOATERHORMANN() also returns the barycentric interpolation
%    weights U.
%
%    [YI,U,YPI,YPPI] = FLOATERHORMANN() also returns the interpolated
%    derivative YPI and the interpolated second derivative YPPI.

% Joe Henning - Fall 2011

% Barycentric Rational Interpolation with no Poles and High Rates of Approximation
% Michael S. Floater and Kai Hormann
% Ifl Technical Report Series, Ifl-06-06

if (nargin < 4)
c = 0;
end

n = length(x);

if (n-1 < c)
fprintf('??? Bad c input to floaterhormann ==> c <= length(x)-1\n');
yi = [];
u = [];
ypi = [];
yppi = [];
return
end

% calculate Floater-Hormann coefficients
u = [];
for k = 1:n
u(k,1) = 0;
for m = k-c:k
prod = 1;
if (m < 1 || m > n-c)
continue
end
for j = m:m+c
if (j ~= k)
prod = prod/(x(k)-x(j));
end
end
u(k,1) = u(k,1) + (-1)^m*prod;
end
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 floaterhormann ==> 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
```