spline.m

% 三次样条插值
% 注释的地方需要修改
clear

xdata=[27.7,28,29,30];          % 插值节点，逗号分隔
ydata=[4.1,4.3,4.1,3.0];        % 插值节点对应的函数值，逗号分隔
xi=28.5;                        % 求值节点
conds='fir';                    % 边界条件类型，'fir'为一阶导数条件，
                                % 'sec'为二阶导数条件，'per'为周期条件
valconds=[3.0,-4.0];            % 边界条件值，为'per'时可不填

n=length(xdata);
hx=diff(xdata);
hy=diff(ydata);
mu=reshape(hx(1:n-2)./(hx(1:n-2)+hx(2:n-1)),[],1);
la=1-mu;
d=6*reshape(diff(hy./hx)./(hx(1:n-2)+hx(2:n-1)),[],1);
pp=zeros(n-1,6);
y=@(x)zeros(size(x));
if strcmpi(conds,'fir')
    M=chase([mu;1],2*ones(n,1),[1;la],...
        [6*(hy(1)/hx(1)-valconds(1))/hx(1);d;
        6*(valconds(2)-hy(n-1)/hx(n-1))/hx(n-1)]);
elseif strcmpi(conds,'sec')
    M1=chase(mu(2:n-2),2*ones(n-2,1),la(1:n-3),...
        [d(1)-mu(1)*valconds(1);d(2:n-3);d(n-2)-la(n-2)*valconds(2)]);
    M=[valconds(1);M1;valconds(2)];
elseif strcmpi(conds,'per')
    lan=hx(1)/(hx(1)+hx(end));
    mun=1-lan;
    A=2*eye(n-1)+diag(la,1)+diag([mu(2:n-2);mun],-1);
    A(1,n-1)=mu(1);
    A(n-1,1)=lan;
    b=[d;6*(hy(1)/hx(1)-hy(end)/hx(end))/(hx(1)+hx(end))];
    M1=A\b;
    M=[M1(n-1);M1];
else
    error('无效的边界条件.')
end
for k=1:n-1
    c=hy(k)/hx(k)-(1/3*M(k)+1/6*M(k+1))*hx(k);
    pp(k,1:4)=[0,0,0,ydata(k)]+...
        [0,0,c*poly(xdata(k))]+...
        [0,1/2*M(k)*poly(xdata([k,k]))]+...
        1/6/hx(k)*(M(k+1)-M(k))*poly(xdata([k,k,k]));
    pp(k,5:6)=xdata(k:k+1);
    y=@(x)y(x)+polyval(pp(k,1:4),x).*(x>=xdata(k) & x<xdata(k+1));
end
yi=feval(y,xi);
yi(xi==xdata(n))=ydata(n);

vpa(yi,10)                      % 结果有效数字位数

function varargout=chase(a,b,c,d)
n=length(d);
for k=2:n
    u=a(k-1)/b(k-1);
    b(k)=b(k)-u*c(k-1);
    d(k)=d(k)-u*d(k-1);
end
x=zeros(n,1);
x(n)=d(n)/b(n);
for k=n-1:-1:1
    x(k)=(d(k)-c(k)*x(k+1))/b(k);
end
[varargout{1:2}]=deal(x,...
    [diag(b)+diag(c,1),d]);
end