Code covered by the BSD License

Interpolation Utilities

Joe Henning (view profile)

22 May 2012 (Updated )

A variety of interpolation utilities

shermite(x, y, yp, ypp, yppp, xi)
```function [yi, ypi, yppi] = shermite(x, y, yp, ypp, yppp, xi)

% SHERMITE 1-D piecewise septic Hermite spline
%    SHERMITE(X,Y,YP,YPP,YPPP,XI) interpolates to find YI, the values of
%    the underlying function Y at the points in the array XI, using
%    piecewise septic Hermite splines.  X and Y must be vectors
%    of length N.
%
%    First, second, and third differences are calculated in a finite
%    difference fashion when YP, YPP, and YPPP are not specified.
%
%    [YI,YPI,YPPI] = SHERMITE() also returns the interpolated
%    septic derivative and quintic second derivative of the underlying
%    function Y at points XI.

% Joe Henning - Fall 2012

if (isempty(yp) && isempty(ypp) & isempty(yppp))
yp = lfindiff (x, y);
ypp = lfindiff (x, yp);
yppp = lfindiff (x, ypp);
elseif (isempty(yp) && isempty(ypp))
yp = lfindiff (x, y);
ypp = lfindiff (x, yp);
elseif (isempty(yp) && isempty(yppp))
yp = lfindiff (x, y);
yppp = lfindiff (x, ypp);
elseif (isempty(yp))
yp = lfindiff (x, y);
elseif (isempty(ypp) && isempty(yppp))
ypp = lfindiff (x, yp);
yppp = lfindiff (x, ypp);
elseif (isempty(ypp))
ypp = lfindiff (x, yp);
elseif (isempty(yppp))
yppp = lfindiff (x, ypp);
else
m = -1;
end

n = length(x);

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 shermite ==> x values must be distinct\n');
yi(i) = NaN;
ypi(i) = NaN;
yppi(i) = NaN;
continue;
end

a = yp(klo);
b = yp(khi);
c = ypp(klo);
d = ypp(khi);
e = yppp(klo);
f = yppp(khi);

% Evaluate septic Hermite polynomial
t = (xi(i) - x(klo))/h;
t2 = t*t;
t3 = t2*t;
t4 = t3*t;
t5 = t4*t;
t6 = t5*t;
t7 = t6*t;
h2 = h*h;
h3 = h2*h;
z0 = 20*t7 - 70*t6 + 84*t5 - 35*t4 + 1;
z1 = 10*t7 - 36*t6 + 45*t5 - 20*t4 + t;
z2 = 2*t7 - 7.5*t6 + 10*t5 - 5*t4 + 1/2*t2;
z3 = 1/6*t7 - 2/3*t6 + t5 - 2/3*t4 + 1/6*t3;
z4 = 1 - z0;
z5 = 10*t7 - 34*t6 + 39*t5 - 15*t4;
z6 = -2*t7 + 6.5*t6 - 7*t5 + 2.5*t4;
z7 = 1/6*t7 - 1/2*t6 + 1/2*t5 - 1/6*t4;

yi(i) = z0*y(klo) + z1*h*a + z2*h2*c + z3*h3*e + z4*y(khi) + z5*h*b + z6*h2*d + z7*h3*f;

% Differentiate to find the sixth-order interpolant
z0 = 140*t6 - 420*t5 + 420*t4 - 140*t3;
z1 = 70*t6 - 216*t5 + 225*t4 - 80*t3 + 1;
z2 = 14*t6 - 45*t5 + 50*t4 - 20*t3 + t;
z3 = 7/6*t6 - 4*t5 + 5*t4 - 8/3*t3 + 1/2*t2;
z4 = -z0;
z5 = 70*t6 - 204*t5 + 195*t4 - 60*t3;
z6 = -14*t6 + 39*t5 - 35*t4 + 10*t3;
z7 = 7/6*t6 - 3*t5 + 5/2*t4 - 2/3*t3;

ypi(i) = z0*y(klo)/h + z1*a + z2*h*c + z3*h2*e + z4*y(khi)/h + z5*b + z6*h*d + z7*h2*f;

% Differentiate to find the fifth-order interpolant
z0 = 840*t5 - 2100*t4 + 1680*t3 - 420*t2;
z1 = 420*t5 - 1080*t4 + 900*t3 - 240*t2;
z2 = 84*t5 - 225*t4 + 200*t3 - 60*t2 + 1;
z3 = 7*t5 - 20*t4 + 20*t3 - 8*t2 + t;
z4 = -z0;
z5 = 420*t5 - 1020*t4 + 780*t3 - 180*t2;
z6 = -84*t5 + 195*t4 - 140*t3 + 30*t2;
z7 = 7*t5 - 15*t4 + 10*t3 - 2*t2;

yppi(i) = z0*y(klo)/h2 + z1*a/h + z2*c + z3*h*e + z4*y(khi)/h2 + z5*b/h + z6*d + z7*h*f;
end

% 3-pt finite difference
function [yp] = lfindiff(x, y)
% Finite Difference Formulae for Unequal Sub-Intervals Using Lagrange's Interpolation Formula
% Singh, Ashok K. and B. S. Bhadauria
% Int. Journal of Math. Analysis, Vol. 3, 2009, no. 17, 815-827

n = length(x);

if (n ~= length(y))
fprintf('Error, len(x) != len(y), returning\n');
yp = NaN;
return;
end

for k = 1:n
if (k == 1)
h1 = x(k+1) - x(k);
h2 = x(k+2) - x(k+1);
yp(k) = -(2*h1+h2)/(h1*(h1+h2))*y(k) + (h1+h2)/(h1*h2)*y(k+1) - h1/(h2*(h1+h2))*y(k+2);
elseif (k == n)
h1 = x(k-1) - x(k-2);
h2 = x(k) - x(k-1);
yp(k) = h2/(h1*(h1+h2))*y(k-2) - (h1+h2)/(h1*h2)*y(k-1) + (h1+2*h2)/(h2*(h1+h2))*y(k);
else
h1 = x(k) - x(k-1);
h2 = x(k+1) - x(k);
yp(k) = -h2/(h1*(h1+h2))*y(k-1) - (h1-h2)/(h1*h2)*y(k) + h1/(h2*(h1+h2))*y(k+1);
end
end
```