cubicspline.m

自行编写的三次样条插值函数源代码

function y=cubicspline(X,Y,m0,mn,x);
N=length(X);
M=zeros(1,N);
M(1)=m0;
M(N)=mn;
f=zeros(1,N-1);
for k=1:N-1
    f(k)=(Y(k+1)-Y(k))/(X(k+1)-X(k))
end;
h=zeros(1,N-1);
r=zeros(1,N-1);
u=zeros(1,N-1);
g=zeros(1,N-1);
for j=1:N-1
    h(j)=X(j+1)-X(j)
end;
for j=2:N-1;
    r(j)=h(j)/(h(j-1)+h(j))
    u(j)=h(j-1)/(h(j-1)+h(j))
    g(j)=3*(r(j)*f(j-1)+u(j)*f(j))
end;
A=zeros(N-2);
b=zeros(N-2,1);
for j=1:N-2
    A(j,j)=2;
    if j==1
        A(1,2)=u(2);
        b(1)=g(2)-r(2)*M(1);
    elseif j==N-2;
        A(N-2,N-3)=r(N-1);
        b(N-2)=g(N-1)-u(N-1)*M(N);
    else
        A(j,j+1)=u(j+1);
        A(j,j-1)=r(j+1);
        b(j)=g(j+1);
    end;
end;
    M(2:N-1)=inv(A)*b;

n=length(x);
y=zeros(1,n);
for j=1:N-1;
    h(j)=X(j+1)-X(j);
end
for j=1:n
    for k=1:N-1;
        if(x(j)>=X(k)&&x(j)<=X(k+1))
            t1=(x(j)-X(k+1))^2*(h(k)+2*(x(j)-X(k)))*Y(k)/h(k)^3;
            t2=(x(j)-X(k))^2*(h(k)+2*(X(k+1)-x(j)))*Y(k+1)/h(k)^3;
            t3=(x(j)-X(k+1))^2*(x(j)-X(k))*M(k)/h(k)^2;
            t4=(x(j)-X(k))^2*(x(j)-X(k+1))*M(k+1)/h(k)^2;
            y(j)=t1+t2+t3+t4;
            break;
        end
    end
end