Splineinterpl.m

function [M_Vdmon,M_cft,Splvalue,Svx]=Splineinterpl(vec_x,vec_y,y01,yn1,x)
% 三次样条函数的多项式计算方法,默认问题1
k=3;
n=length(vec_x);p=length(x);
A=zeros(n+2,k+1);H=zeros(n+2,n-2);Bp=zeros(p,n+k-1);

    j=0:k;
    A(1:n,:)=vec_x'.^j;
    A(n+1,:)=j.*vec_x(1).^(j-1);%一阶求导前半部，左边界点取值
    A(n+2,:)=j.*vec_x(n).^(j-1);%一阶求导前半部，右边界点取值

for i=3:n
    H(i,1:i-2)=(vec_x(i)-vec_x(2:i-1)).^k;
end

    H(n+1,:)=0;
    H(n+2,:)=k*(vec_x(n)-vec_x(2:n-1)).^(k-1);%一阶求导后半部，右边界点取值
    M_Vdmon=[A,H];
    yb=[vec_y,y01,yn1];
    M_cft=M_Vdmon\yb';%拟合系数M_cft=[a0,a1,...c1,c2...cn-1]

for s=1:p
     b=x(s)-vec_x(2:n);
     [b2,~]=findplus(b);
     B=x(s).^j;
     Bp(s,:)=[B,b2.^k];
end

Splvalue=Bp*M_cft;%三次样条逼近值

digits(5);

syms t

bt=t-vec_x;

Bt=t.^j;
Btm=Bt*M_cft(1:k+1);

Svx(1,2)=vpa(Btm);Svx(1,1)=1;st=0;
for r=2:n-1
    st=st+bt(r)^k*M_cft(k+r);
    Svx(r,2)=vpa(st+Btm);Svx(r,1)=r;
end