三次样条插值.txt

来自「数值计算方法的一些题」· 文本 代码 · 共 36 行

function S=myspline(x,y,dy1,dy2,xi)
[m,n]=size(x);
for i=1:n-1
    h(i)=x(i+1)-x(i);
end
for i=1:n-1
    b(i)=(y(i+1)-y(i))/h(i);
end
d(1)=b(1)-dy1;
d(n)=dy2-b(n-1);
for i=2:n-1
    d(i)=b(i)-b(i-1);
end
for i=2:n-1
    B(i,i)=2*(h(i-1)+h(i));
end
 B(1,1)=2*h(1);
 B(n,n)=2*h(n-1);
for i=2:n
    B(i,i-1)=h(i-1);
end
for i=1:n-1
    B(i,i+1)=h(i);
end
B(:,n+1)=6*d';
m=mypursue(B);
for i=1:n-1
    A(1,i)=y(i);
    A(2,i)=b(i)-(m(i+1)+2*m(i))*h(i)/6;
     A(3,i)=m(i)/2;
     A(4,i)=(m(i+1)-m(i))/(6*h(i));
end
A=A';
k=find((xi-x)<0);
k=k(1)-1;
S=A(k,1)+A(k,2)*(xi-x(k))+A(k,3)*(xi-x(k))^2+A(k,4)*(xi-x(k))^3;