function output = hermite(x,y,s,xx) % yi = hermite(x,y,s,xi) % For given vectors x and y find the piecewise cubic Hermite interpolation % and evaluate it at the points xi, yielding the vector yi. % % pp = hermite(x,y,s) returns the pp-form of the interpolant instead. % This can be used to evaluate the interpolant later at xi % using yi = ppval(pp,xi). % % EXAMPLE: for sin(x) and the nodes 0, pi/2, pi evaluate the interpolating function % at 1000 equidistant points in [0,pi] and plot it. % % x = [0,pi/2,pi]; y = sin(x); s = cos(x); % xi = linspace(0,pi,1000); yi = hermite(x,y,s,xi); % plot(xi,sin(xi),xi,yi); legend('sin(x)','Hermite interpolation') output=[]; n=length(x); if n<2, error('There should be at least two data points.'), end if any(diff(x)<0), [x,ind]=sort(x); else, ind=1:n; end x=x(:); dx = diff(x); if all(dx)==0, error('The data abscissae should be distinct.'), end [yd,yn] = size(y); % if Y happens to be a column matrix, change it to % the expected row matrix. if yn==1, yn=yd; y=reshape(y,1,yn); s=reshape(s,1,yn); yd=1; end yi=y(:,ind).';s=s(:,ind).'; dd = ones(1,yd); dx = diff(x); divdif = diff(yi)./dx(:,dd); % convert to pp form c4=(s(1:n-1,:)+s(2:n,:)-2*divdif(1:n-1,:))./dx(:,dd); c3=(divdif(1:n-1,:)-s(1:n-1,:))./dx(:,dd) - c4; pp=mkpp(x.',... reshape([(c4./dx(:,dd)).' c3.' s(1:n-1,:).' yi(1:n-1,:).'],(n-1)*yd,4),yd); if nargin==3 output=pp; else output=ppval(pp,xx); end