三次样条插值.docx

三次样条插值
算法原理
由于在许多实际问题中，要求函数的二阶导数连续，人们便提出了三次样条插值函数，三次样条插值函数是由分段三次函数拼接而成的，在连接点处二阶导数连续。
设S(x)在节点处的二阶导数，其中为待定参数。由S（x）是分段三次多项式可知，是分段线性函数，在子区间上可以表示为

其中，对S（x）两端积分两次得

其中和为积分常数。由插值条件得

由此解得

代入得：

求导得：

令得在处的左导数

又令得在处右导数
，
从而有，由在节点处一阶导数的连续性知，即

两端同乘得，
记，，则关于的方程组写成。
三种边界条件的三弯矩方程：
（1）第一种边界条件：已知。取，这时方程组减少了两个未知量，变成只含n-1个未知量的n-1个方程的方程组，用矩阵表示为

可用追赶法求解出后，即得三条样条插值函数。
第二种边界条件，已知，记，则有，得
，即
，其中，得到方程组，用矩阵表示为
，该方程组的系数矩阵是严格三对角占优矩阵，可用追赶法求解。
（3）第三种边界条件：周期型边界条件。已知是以为周期的周期函数，则由周期性可知，，这时将点看成是内节点，则有，也即
，其中，方程组第1个方程为：，所以方程组为
，用矩阵表示为
，显然系数矩阵为严格对角占优矩阵，可用LU分解法求解。

程序框图

源程序
functionx=mchase(A,d)
%追赶法
n=length(d);
u=zeros(n,1);u(1)=A(1,1);
fork=2:n
l(k)=A(k,k-1)/u(k-1);
u(k)=A(k,k)-l(k)*A(k-1,k);
end
y=zeros(n,1);y(1)=d(1);
fori=2:n
y(i)=d(i)-l(i)*y(i-1);
end
x=zeros(n,1);
x(n)=y(n)/u(n);
fori=n-1:-1:1
x(i)=(y(i)-A(i,i+1)*x(i+1))/u(i);
end
x
end

functionT=mspline1(x0,y0,f21,f22,xx)
%三次样条插值函数第一种边界条件
%x0、y0分别为节点的横坐标和纵坐标；
%f21为左端点的二阶导数值,f22为右端点的二阶导数值；xx为由插值点组成的向量
n=length(x0)-1;%计算小区间数
fori=1:n
h(i)=x0(i+1)-x0(i);
end
fori=1:n-1
mu(i)=h(i)/(h(i)+h(i+1));
lamda(i)=1-mu(i);
d(i)=6*((y0(i+2)-y0(i+1))/h(i+1)-(y0(i+1)-y0(i))/h(i))/(h(i)+h(i+1));
end
A=zeros(n-1);
fori=1:n-2
A(i+1,i)=mu(i+1);%次下对角线
A(i,i+1)=lamda(i);%次上对角线
A(i,i)=2;%主对角线
end
A(n-1,n-1)=2;
dd=zeros(n-1,1);%右端列向量
fori=2:n-2
dd(i)=d(i);
end
dd(1)=d(1)-mu(1)*f21;dd(n-1)=d(n-1)-lamda(n-1)*f22;
M=mchase(A,dd);%追赶法求解M值
h
mu
lamda
A
dd
M=[f21,M',f22]'
t=sym('t');
a=zeros(n,1);b=zeros(n,1);c=zeros(n,1);e=zeros(n,1);
fori=1:n
a(i)=M(i)./(6*h(i));
b(i)=M(i+1)./(6*h(i));
W1(i)=b(i)-a(i);
W2(i)=3*(a(i).*x0(i+1)-b(i).*x0(i));
c(i)=y0(i)./h(i)-h(i).*M(i)/6;
e(i)=y0(i+1)./h(i)-h(i).*M(i+1)/6;
W3(i)=3*b(i).*x0(i).^2-3*a(i).*x0(i+1).^2+e(i)-c(i);
W4(i)=a(i).*x0(i+1).^3-b(i).*x0(i).^3+c(i).*x0(i+1)-e(i).*x0(i);
Si(t)=W1(i).*t^3+W2(i).*t^2+W3(i).*t+W4(i)%每个小区间的三次样条差值函数表达式
end
m=length(xx);T=zeros(m,1);
fork=1:m
forj=1:n
if((xx(k)>=x0(j))&(xx(k)<x0(j+1)))
T(k)=W1(j).*(xx(k).^3)+W2(j).*(xx(k).^2)+W3(j).*xx(k)+W4(j);
end
end
end
T
End

functionT=mspline2(x0,y