spline3.m

一个时滞系统的工具箱

function yi=spline3(x,y,xi);
%SPLINE3 Interpolation and extrapolation by degenerate cubic splines.
%	 Given data vectors x and y and a new absissa vector xi,
%	 the function yi=spline3(x,y,xi) uses interpolation and,
%	 if needed, extrapolation by the degenerate splines 
%	 to find a vector yi corresponding to xi.
 
% Interpolation and extrapolation for the simplest cases
n=length(x);
if n==0  return, end;
if n==1 
   yi=ones(length(xi),1)*y;
   return
end;
if n==2
   L1=(xi(:)-x(2))/(x(1)-x(2));
   L2=(xi(:)-x(1))/(x(2)-x(1));
   yi=[L1, L2]*y;
   return
end;
if n==3
   L1=(xi(:)-x(2)).*(xi(:)-x(3))/((x(1)-x(2))*(x(1)-x(3)));
   L2=(xi(:)-x(1)).*(xi(:)-x(3))/((x(2)-x(1))*(x(2)-x(3)));
   L3=(xi(:)-x(1)).*(xi(:)-x(2))/((x(3)-x(1))*(x(3)-x(2)));
   yi=[L1 L2 L3]*y;
   return
end;

m = length(xi);
if  m == 0,  return, end;
yi = zeros(m,size(y,2));

% Partition of array x
temp = rem(n,3);
if temp == 0,       k=[1, 2, 3:3:n];
elseif temp == 1,   k=1:3:n;
elseif temp == 2,   k=[1, 2:3:n];
end;

% Building the spline for every point of array xi
for l=1:m
  xx = xi(l);
  j=find(x(k) > xx);
  if isempty(j) 
    % Extrapolation to the right side from array x
    i=n;
  else
    % Interpolation and extrapolation to the left side from array x
    i=k(j(1));
  end;
  if i >= 4  
    L1=(xx-x(i-2))*(xx-x(i-1))*(xx-x(i))/((x(i-3)-x(i-2))*(x(i-3)-x(i-1))*(x(i-3)-x(i)));
    L2=(xx-x(i-3))*(xx-x(i-1))*(xx-x(i))/((x(i-2)-x(i-3))*(x(i-2)-x(i-1))*(x(i-2)-x(i)));
    L3=(xx-x(i-3))*(xx-x(i-2))*(xx-x(i))/((x(i-1)-x(i-3))*(x(i-1)-x(i-2))*(x(i-1)-x(i)));
    L4=(xx-x(i-3))*(xx-x(i-2))*(xx-x(i-1))/((x(i)-x(i-3))*(x(i)-x(i-2))*(x(i)-x(i-1)));
    yi(l,:) = [L1 L2 L3 L4]*[y(i-3,:); y(i-2,:); y(i-1,:); y(i,:)];
  else
    L1=(xx-x(2))*(xx-x(3))*(xx-x(4))/((x(1)-x(2))*(x(1)-x(3))*(x(1)-x(4)));
    L2=(xx-x(1))*(xx-x(3))*(xx-x(4))/((x(2)-x(1))*(x(2)-x(3))*(x(2)-x(4)));
    L3=(xx-x(1))*(xx-x(2))*(xx-x(4))/((x(3)-x(1))*(x(3)-x(2))*(x(3)-x(4)));
    L4=(xx-x(1))*(xx-x(2))*(xx-x(3))/((x(4)-x(1))*(x(4)-x(2))*(x(4)-x(3)));
    yi(l,:) = [L1 L2 L3 L4]* y(1:4,:);
  end;
end;