Code covered by the BSD License

# Interpolation Utilities

### Joe Henning (view profile)

22 May 2012 (Updated )

A variety of interpolation utilities

newtint(x, y, xi, c)
```function [yi, p, pval] = newtint(x, y, xi, c)

% NEWTINT Interpolation of equally-spaced points.
%    NEWTINT(X,Y,XI,C) interpolates to find YI, the value
%    of the underlying function Y at the point XI, using
%    either Newton's forward interpolation formula or
%    Newton's backward interpolation formula.
%
%    C specifies the interpolation method:
%       0 : Newton's forward or backward (default)
%       1 : Newton's forward
%       2 : Newton's backward
%       3 : Stirling's method
%       4 : Everett's method
%
%    [YI,P,PVAL] = NEWTINT() also returns P, the coefficients
%    of the calculated interpolating polynomial, and PVAL,
%    which specifies the maximum degree of the interpolating
%    polynomial.

% Joe Henning - Fall 2011

if (nargin < 4)
c = 0;
end

pval = length(x)*(0+1) - 1;

n = length(x)-1;
D = zeros(n+1,n+1);
D(:,1) = y(:);

h = x(2)-x(1);

for i = 1:n
for j = 1:i
D(i+1,j+1) = D(i+1,j)-D(i,j);
end
end

if (c == 0)
if (xi < x(1))
c = 2;
elseif (xi > x(ceil((n+1)/2)) && xi <= x(n+1))
c = 2;
else
c = 1;
end
end

switch c
case 1
% forward difference interpolation
p = diag(D);
q = (xi-x(1))/h;
yi = p(1);
for i = 1:length(p)-1
term = 1;
for k = 0:(i-1)
term = term*(q-k);
end
term = term/factorial(i);
yi = yi + term*p(i+1);
end
case 2
% backward difference interpolation
p = D(n+1,:).';
q = (xi-x(n+1))/h;
yi = p(1);
for i = 1:length(p)-1
term = 1;
for k = 0:(i-1)
term = term*(q+k);
end
term = term/factorial(i);
yi = yi + term*p(i+1);
end
case 3
% Stirling's method
xk = 0;
ik = 0;
for i = 1:n+1
if (xi > x(i))
xk = x(i);
ik = i;
end
end
yi = 0;
q = (xi-xk)/h;
if (q > 0.5)
ik = ik + 1;
xk = x(ik);
q = (xi-xk)/h;
for i = 1:2:n+1
fi = floor(i/2);
term = 0;
%[term D(ik+i-1-floor(i/2),i)]s
for k = fi-1:fi
term = term + dnchoosek(q+k,i-1)/2.0;
end
p(i,1) = term;
if (ik+i-1-floor(i/2) <= n+1)
yi = yi + term*D(ik+i-1-floor(i/2),i);
end
end
for i = 2:2:n+1
fi = floor(i/2);
term = dnchoosek(q+fi-1,i-1);
p(i,1) = term;
%[term D(ik+i-1-floor(i/2),i) D(ik+i-floor(i/2),i)]
if (ik+i-floor(i/2) > n+1)
yi = yi + term*D(ik+i-1-floor(i/2),i)/2.0;
else
yi = yi + term*(D(ik+i-1-floor(i/2),i)+D(ik+i-floor(i/2),i))/2.0;
end
end
else
for i = 1:2:n+1
fi = floor(i/2);
term = 0;
for k = fi-1:fi
term = term + dnchoosek(q+k,i-1)/2.0;
end
p(i,1) = term;
%[term D(ik+i-1-floor(i/2),i)]
if (ik+i-1-floor(i/2) <= n+1)
yi = yi + term*D(ik+i-1-floor(i/2),i);
end
end
for i = 2:2:n+1
fi = floor(i/2);
term = dnchoosek(q+fi-1,i-1);
p(i,1) = term;
%[term D(ik+i-1-floor(i/2),i) D(ik+i-floor(i/2),i)]
if (ik+i-floor(i/2) > n+1)
yi = yi + term*D(ik+i-1-floor(i/2),i)/2.0;
else
yi = yi + term*(D(ik+i-1-floor(i/2),i)+D(ik+i-floor(i/2),i))/2.0;
end
end
end
case 4
% Everett's method
xk = 0;
ik = 0;
for i = 1:n+1
if (xi > x(i))
xk = x(i);
ik = i;
end
end
p = [];
q = (xi-xk)/h;
r = 1-q;
yi = 0;
for i = 1:2:n+1
fi = floor(i/2);
e = dnchoosek(r+fi,i);
f = dnchoosek(q+fi,i);
p = [p; e f];
%[e D(ik+i-1-floor(i/2),i) f D(ik+i-floor(i/2),i)]
if ((ik+i-floor(i/2)) > n+1)
yi = yi + e*D(ik+i-1-floor(i/2),i) + f*0;
else
yi = yi + e*D(ik+i-1-floor(i/2),i) + f*D(ik+i-floor(i/2),i);
end
end
end

function [d] = dnchoosek(n, k)

% DNCHOOSEK Binomial coefficient or all combinations.
%    DNCHOOSEK(N,K) returns the binomial coefficient
%    (N K) = N!/K!(N-K)!.

if (k == 0 || k == n)
d = 1;
return
end

if (k == 1 || k == n-1)
d = n;
return
end

d = gamma(n+1)/(gamma(k+1)*gamma(n-k+1));
```