三次样条插值的算法推导与matlab实现

最近有一个数值分析上机作业，脑子一热便用matlab做了实现（仅考虑了边界条件为一阶导数），下面先对求解过程进行推导。

目录

算法推导

matlab实现

---

算法推导

个人感觉数值分析教材已经比较详细了，把我用的教材贴个图：（喻文健教授编写的，第三版）

matlab实现

matlab写的my_spline函数代码如下：

function [Y,p]=my_spline(x,y,xx,flag,condition)
%x，y为插值节点坐标
%xx为需要输出S(x)值得x坐标
%flag为边界条件选取，仅做了1，边界为一阶导数情况
%condition为边界条件

%Y为XX对应S(x)值
%p为个段多项式得系数

n=length(x);
u=ones(1,n-1);
q=ones(1,n-1);
A=zeros(n);
d=zeros(1,n);
h=zeros(1,n-1);

for i=1:1:n-1
    h(i)=x(i+1)-x(i);
end
for i=1:1:n-2
    u(i)=h(i)/(h(i)+h(i+1));
    q(i+1)=1-u(i);
end

%自然边界
if(flag==1)
    for i=1:1:n
        for j=1:1:n
            if(i==j)
                A(i,j)=2;
            elseif(i-j==1)
                A(i,j)=u(j);
            elseif(i-j==-1)
                A(i,j)=q(i);
            end
        end
    end
    d(1)=6/h(1)*((y(2)-y(1))/h(1)-condition(1));
    d(n)=6/(h(length(h)))*(condition(2)-(y(n)-y(n-1))/(h(length(h))));
    for i=2:1:n-1
        d(i)=6*(y(i-1)/(h(i-1)*(h(i-1)+h(i)))+y(i+1)/(h(i)*(h(i-1)+h(i)))-y(i)/(h(i-1)*h(i)));
    end
end

%求解AM=d
M=A\d';
%计算并将S(x)系数放入矩阵
p=zeros(4,n-1)';
for i=1:1:n-1
    a=(y(i+1)-y(i))/h(i)-(M(i+1)-M(i))*h(i)/6;
    b=(y(i)*x(i+1)-y(i+1)*x(i))/h(i)+(M(i+1)*x(i)-M(i)*x(i+1))*h(i)/6;
    p(i,1)=(M(i+1)-M(i))/(6*h(i));
    p(i,2)=(M(i)*x(i+1)-M(i+1)*x(i))/(2*h(i));
    p(i,3)=(M(i+1)*x(i)^2-M(i)*x(i+1)^2)/(2*h(i))+a;
    p(i,4)=(M(i)*x(i+1)^3-M(i+1)*x(i)^3)/(6*h(i))+b;
end
Y=zeros(1,length(xx));
for i=1:1:length(xx)
    for j=1:1:n-1
        if(xx(i)>=x(j)&&xx(i)<=x(j+1))
            Y(i)=p(j,1)*xx(i)^3+p(j,2)*xx(i)^2+p(j,3)*xx(i)+p(j,4);
        end
    end
    
end

end

 调用实例代码：

clc
clear
clf
x=[0.52 3.1 8.0 17.95 28.65 39.62]；
y=[5.288 9.4 13.84 20.2 24.9 28.44];
condition=[1.86548,-0.046115];
xx=linspace(0.52,520,2000);
[y1,p]=my_spline(x,y,xx,1,condition);
plot(x,y,'x',xx,y1)

实测可用