cubic_spline.m

三次样条插值的MATLAB程序（三弯矩方程）

function yi=cubic_spline(x,y,ydot,xi)
%三次样条插值公式（三弯矩方程，第一类边界条件），其中，
%x为向量,全部的插值节点;
%y为向量,插值节点处的函数值;
%ydot为向量，端点处的导数值,
%如果此处值缺省,则用均差代替导数,
%xi为标量，自变量x;
%yi为xi处的函数估计值。
n=length(x); ny=length(y);
%输入的插值点与它的函数值应有相同的个数。
if n~=ny
    error('The lengths of X and Y must be equal!');
    return;
end
%如果没有给出y的导数值，则用均差代替导数
if isempty(ydot)==1
    ydot=[(y(2)-y(1))/(x(2)-x(1)) (y(n)-y(n-1))/(x(n)-x(n-1))];
end
h=zeros(1,n);lambda=ones(1,n);mu=ones(1,n);
M=zeros(n,1);d=zeros(n,1);
for k=2:n
    h(k)=x(k)-x(k-1);
    %插值点必须互异
    if abs(h(k))<eps
        error('the DATA is error!');
        return;
    end
end
for k=2:n-1
    lambda(k)=h(k+1)/(h(k)+h(k+1)); mu(k)=1-lambda(k);
    d(k)=6/(h(k)+h(k+1))...
        *((y(k+1)-y(k))/h(k+1)-(y(k)-y(k-1))/h(k));
end
d(1)=6/h(2)*((y(2)-y(1))/h(2)-ydot(1));
d(n)=6/h(n)*(ydot(2)-(y(n)-y(n-1))/h(n));
A=diag(2*ones(1,n));
for i=1:n-1
    A(i,i+1)=lambda(i);A(i+1,i)=mu(i+1);
end
M=A\d;
for k=2:n
    if x(k-1)<=xi & xi<=x(k)
        yi=M(k-1)/6/h(k)*((x(k)-xi)^3)...
            +M(k)/6/h(k)*((xi-x(k-1))^3)...
            +1/h(k)*(y(k)-M(k)*h(k)^2/6)*(xi-x(k-1))...
            +1/h(k)*(y(k-1)-M(k-1)*h(k)^2/6)*(x(k)-xi);
        return;
    end
end