drawspline.m

来自「Incorporating Prior Knowledge in Cubic S」· M 代码 · 共 53 行

function [yv,dyv] = drawspline(sp,xv,s)
% DRAWSPLINE Cubic spline interpolation function
%
% [yv,dyv] = drawspline(sp,xv)  or
% [yv,dyv] = drawspline(sp,xv,s)
%
%   sp: spline
%   xv: vector of points where funtcion will interpolate yv and dyv
%   s:  character string for PLOT function
%          (if you don't use it, function won't draw)
%   yv: spline values
%   dyv: spline derivatives
%
% example:
%       hold on
%       drawspline(sp,[0:0.1:10],'b-');

yv = zeros(size(xv));
dyv = yv;
n = length(sp.x);
m = length(xv);
i = 1;
j = 1;
while (j<=m),
  while (i<n-1) & (xv(j)>=sp.x(i+1)),
    i = i+1;
  end
  if (xv(j)<sp.x(i)) | (xv(j)>sp.x(i+1)),
    warning('Out of intervals (xv)');
    yv(j)=0;
    i = 1;
  else
    k1 = sp.x(i); k2 = sp.x(i+1);
    xx = xv(j); h = k2-k1;
    a = (k2-xx)^2*(xx-k1)/h^2;  
    b = -(xx-k1)^2*(k2-xx)/h^2;
    c = (k2-xx)^2*(2*(xx-k1)+h)/h^3;
    d = (xx-k1)^2*(2*(k2-xx)+h)/h^3;
    da = (-2*(k2-xx)*(xx-k1)+(k2-xx)^2)/h^2;
    db = -(2*(xx-k1)*(k2-xx)-(xx-k1)^2)/h^2;
    dc = (-2*(k2-xx)*(2*(xx-k1)+h)+(2*(k2-xx)^2))/h^3;
    dd = (2*(xx-k1)*(2*(k2-xx)+h)-(2*(xx-k1)^2))/h^3;
    yv(j) = c*sp.y(i) + a*sp.dy(i) + d*sp.y(i+1) + b*sp.dy(i+1);
    dyv(j) = dc*sp.y(i) + da*sp.dy(i) + dd*sp.y(i+1) + db*sp.dy(i+1);
  end
  j = j+1;
end

% Drawing
if (nargin==3),
  plot(xv,yv,s);
end