selfspline.m

第一、二种边界条件的三次样条插值程序

function [s,m]=selfspline(x0,y0,df,x,conds)
%函数selfspline是三次样条函数按第一和第二种边界条件求x所对应的函数值s.
%其中x是自变量的值构成的向量,输出的m为弯矩向量.
%x0为插值节点构成的向量,y0为对应于插值节点的函数值构成的向量.
%conds为边界条件,其值取1时, 为第一种边界条件,此时df为边界点的一阶
%导数值构成的向量;其值取2时, 为第二种边界条件,df为边界点的二阶导数值构成的向量.
n=length(x0);
h=diff(x0);
b=ones(1,n)*2;
for j=2:n-1
    a(j)=h(j-1)/(h(j-1)+h(j));
    c(j)=h(j)/(h(j-1)+h(j));
    d(j)=6*(h(j-1)*(y0(j+1)-y0(j))-h(j)*(y0(j)-y0(j-1)))/(h(j)*h(j-1)*(h(j)+h(j-1)));
end
a(1:n-2)=a(2:n-1);
switch conds
case 1
       a(n-1)=1;
       c(1)=1;
       d(1)=6*((y0(2)-y0(1))/(x0(2)-x0(1))-df(1))/(x0(2)-x0(1));     
       d(n)=6*(df(2)-(y0(n)-y0(n-1))/(x0(n)-x0(n-1)))/(x0(n)-x0(n-1)); 
case 2
       a(n-1)=0;
       c(1)=0;
       d(1)=2*df(1);
       d(n)=2*df(2);
otherwise
       error('conds的值输入错误，conds的值只能取1或者2')
end
for i=2:n
    r=a(i-1)/b(i-1);
    b(i)=b(i)-r*c(i-1);
    d(i)=d(i)-r*d(i-1);
end
m(n)=d(n)/b(n);
for i=n-1:-1:1
    m(i)=(d(i)-c(i)*m(i+1))/b(i);
end
for j=1:length(x)
    for i=1:n-1
        if x(j)>=x0(i)&x(j)<=x0(i+1)
           s(j)=m(i)*(x0(i+1)-x(j))^3/(6*h(i))+m(i+1)*(x(j)-x0(i))^3/(6*h(i))+...
                (y0(i)-m(i)*h(i)^2/6)*(x0(i+1)-x(j))/h(i)+...
                (y0(i+1)-m(i+1)*h(i)^2/6)*(x(j)-x0(i))/h(i);        
        end
    end
end