function
v = pchiptx(x,y,u)
%PCHIPTX Textbook piecewise cubic Hermite interpolation.
% v = pchiptx(x,y,u) finds the shape-preserving piecewise
% P(x), with P(x(j)) = y(j), and returns v(k) = P(u(k)).
% See PCHIP, SPLINETX.
% First
derivatives
h = diff(x);
delta = diff(y)./h;
d = pchipslopes(h,delta);
%
Piecewise polynomial coefficients
n = length(x);
c = (3*delta - 2*d(1:n-1) -
d(2:n))./h;
b = (d(1:n-1) - 2*delta +
d(2:n))./h.^2;
% Find
subinterval indices k so that x(k) <= u
< x(k+1)
k = ones(size(u));
for j = 2:n-1
k(x(j) <= u) = j;
end
%
Evaluate interpolant
s = u - x(k);
v = y(k) + s.*(d(k) + s.*(c(k) + s.*b(k)));
function d
= pchipslopes(h,delta)
% PCHIPSLOPES Slopes for shape-preserving Hermite cubic
% pchipslopes(h,delta) computes d(k) = P'(x(k)).
% Slopes at interior points
% delta = diff(y)./diff(x).
% d(k) = 0 if delta(k-1) and delta(k) have opposites signs or
either is zero.
% d(k) = weighted harmonic mean of delta(k-1) and delta(k) if they
have the same sign.
n = length(h)+1;
d = zeros(size(h));
k =
find(sign(delta(1:n-2)).*sign(delta(2:n-1))>0)+1;
w1 = 2*h(k)+h(k-1);
w2 = h(k)+2*h(k-1);
d(k) = (w1+w2)./(w1./delta(k-1) + w2./delta(k));
% Slopes at
endpoints
d(1) = pchipend(h(1),h(2),delta(1),delta(2));
d(n) = pchipend(h(n-1),h(n-2),delta(n-1),delta(n-2));
function d = pchipend(h1,h2,del1,del2)
%
Noncentered, shape-preserving, three-point formula.
d = ((2*h1+h2)*del1 - h1*del2)/(h1+h2);
if sign(d) ~= sign(del1)
d = 0;
elseif
(sign(del1)~=sign(del2))&(abs(d)>abs(3*del1))
d = 3*del1;
end