三次样条插值函数matlab代码实现

xi=[27.7 28 29 30]
yi=[4.1 4.3 4.1 3]
si=[3 -4]
diffq(xi,yi)
s=splin3(xi,yi,si)
simplify(s)
%均差表计算
function p=diffq(xi,yi)
n=length(xi);
f=ones(n,n);
%对差商表第一列赋值
for kk=1:n
    f(kk)=yi(kk);
end

for i1=2:n
    for j1=i1:n
        f(j1,i1)=(f(j1,i1-1)-f(j1-1,i1-1))/(xi(j1)-xi(j1-(i1-1)));
    end
end
disp('均差表如下：');
disp(f);
p=f;
end
%三次样条
function s=splin3(xi,yi,si)
p=diffq(xi,yi)
%xi是插值节点横坐标
%yi是插值节点纵坐标
%si是边界条件
n=length(xi);
a=zeros(n,n);
hi=zeros(n-1,1);
%第一种边界条件
u=zeros(n-1,1);
r=zeros(n-1,1);
d=zeros(n,1);
for i=1:n-1
    hi(i)=xi(i+1)-xi(i);
end
u(n-1)=1;
for ii=1:n-2
    u(ii)=(hi(ii)/(hi(ii)+hi(ii+1)));
end
r(1)=1;
for iii=2:n-1
    r(iii)=(hi(iii)/(hi(iii-1)+hi(iii)));
end
d(1)=(6*(p(2,2)-si(1)))/hi(1);
d(n)=(6*(si(2)-(p(n,2))))/hi(n-1);
for i2=2:n-1
    d(i2)=6*p(i2+1,3);
end
a(n,n)=2;
for k=1:n-1
    for j=1:n-1
        if(k==j)
            a(k,j)=2;
            a(k,j+1)=r(k);
            a(k+1,j)=u(j);
        end
    end
end
m=a\d;
syms x
c=zeros(n,1);
b=sym(c);
for i=1:n-1
    b(i,1)=(m(i)*(xi(i+1)-x)^3)/(6*hi(i))+(m(i+1)*(x-xi(i))^3)/(6*hi(i))+((yi(i)-(m(i)*(hi(i)^2)/6))*(xi(i+1)-x)/hi(i))+((yi(i+1)-(m(i+1)*(hi(i)^2)/6))*(x-xi(i))/hi(i));
end 
s=b;  
end