hermiteint.m

Interpolation routines in matlab

%
% Compute Hermite interpolation polynomial.
%
% function [a] = HermiteInt( x, f )
%
% input:
%        x:     vector containing the 
%               interpolation points
%        f:     vector containing the values 
%               to be interpolated
%
% output:
%        a:     vector of coefficients of Newton's 
%               interpolating polynomial
%
%

function [a] = HermiteInt( x, f );


N = size(x(:),1);


if size(f(:),1) ~= N 
   ifail = 1;
   return
end

A = zeros(N,N);


%  Compute first column of divided difference table

z      = x(1);
A(1,1) = f(1); 
i1     = 1;             %  New interpolation point appears first at index i1 
for j = 2:N
    if ( z == x(j) )
       A(j,1) = f(i1);
    end
    if ( z ~= x(j) ) 
       z      = x(j);   % New interpolation point
       i1     = j;
       A(j,1) = f(i1);
     end
end  

%  Compute remaining columns of the divided differences table

fac = 1;
for j = 2:N
    fac = fac * (j-1);   % fac = (j-1)!
    i1  = j;             % first time equal nodes are detected
    for i = j:N
        if ( x(i) == x(i-j+1) )
             A(i,j) = f(i1) / fac ;
        end
        if ( x(i) ~= x(i-j+1) )
             i1 = i+1;
             A(i,j) = ( A(i,j-1) - A(i-1,j-1) ) / (x(i) - x(i-j+1)) ;
        end
    end
end

for j = 1:N
    a(j) = A(j,j);
end

A