matlab-三次样条插值

function [ yi ] = sanciyangtiao( x,y,xi,dx1,dx2 )
x=input(' please input x ='); % 输入插值节点
y=input(' please input y ='); % 输入插值节点
dx1=input(' please input dx1 ='); % 输入插值节点
dx2=input(' please input dx2 ='); % 输入插值节点
xi=input(' please input xi =');% 输入要求的x
if length(x)==length(y)
n=length(x);
else
error('the data of x and y must be equal')
end
H=diff(x);
D=zeros(n);
D(1)=6*((y(2)-y(1))/H(1)-dx1)/H(1);
D(n)=6*(dx2-(y(n)-y(n-1))/H(n-1))/H(n-1);
for i=2:n-1
D(i,1)=6*((y(i+1)-y(i))/H(i)-(y(i)-y(i-1))/H(i-1))/(H(i)+H(i-1));
end
q=zeros(n);
for i=1:n
q(i,i)=2;
end
q(2,1)=1;
q(n,n-1)=1;
for i=1:n-2
q(i+1,i)=H(i)/(H(i+1)+H(i));
end
for i=2:n-1
q(i,i+1)=H(i)/(H(i)+H(i-1));
end
m=zeros(n,1);
m=q\D;
for i=1:n
if xi>x(i)
S=m(i,1)*(x(i+1)-xi)^3/(6*H(i))+m(i+1,1)*(xi-x(i))^3/(6*H(i))+...
(y(i)-(m(i,1)*H(i)^2)/6)*(x(i+1)-xi)/H(i)+(y(i+1)-(m(i+1,1)*H(i)^2)/6)*(xi-x(i))/H(i);
end
end
yi=S;
end