cfspline.m

Interpolation routines in matlab

%
% Compute the free cubic spline.
%
% function [a, b, c, d] = cfspline( x, f )
%
% input:
%        x:     vector containing the 
%               interpolation points
%
%        f:     vector containing the 
%               function values to be interpolated 
%
%
% output:
%
%        a:     vector containing the 
%               coefficients a(i)
%
%        b:     vector containing the 
%               coefficients b(i)
%
%        c:     vector containing the 
%               coefficients c(i)
%
%        d:     vector containing the 
%               coefficients d(i)
%
%        iflag  error flag
%               iflag = 0  spline is computed
%               iflag = 1  number of interpolation points and function
%                          values are not the same

%

function [a, b, c, d, iflag] = cfspline( x, f );

iflag = 0;

n = size(x(:),1);
if ( n ~= size(f(:),1) )
   iflag = 1;
   return
end

% compute interval length
for i = 1:n-1
    h(i) = x(i+1) - x(i);
end

% Compute the  (n-2) by (n-2)  tridiagonal matrix  T
for i = 1:n-2
    d(i) = 2*( h(i)+h(i+1) );
    e(i) = h(i+1);
end

% right hand side of the system (stored in a)
t1   =  3*(f(2) - f(1)) / h(1);
for i = 1:n-2
    t2   =  3*(f(i+2) - f(i+1)) / h(i+1);
    a(i) = t2 - t1;
    t1   = t2;
end

%  solve  T c = a
[d, e, ierr] = tridiagLDL( d, e );    % factor T
[c, ierr] = tridiagLDLsl( d, e, a );  % solve  T c = a

%  shift vector c by one
for i = n-2:-1:1
    c(i+1) = c(i);
end

% boundary conditions for free cubic spline
c(1) = 0;
c(n) = 0;


%  compute a(i), d(i)
for i = 1:n-1
    a(i) = f(i);
    d(i) = (c(i+1) - c(i)) / (3*h(i));
end
a(n) = f(n);

%  compute b(i)
for i = 1:n-1
    b(i) = (a(i+1) - a(i))/h(i) - h(i)*(2*c(i) + c(i+1))/3;
end