样条插值实验报告

四、三次样条插值1.样条函数插值的原理给定区间[,]ab上划分011:nnaxxxxb，若分段函数()Sx满足：1.()Sx在各个子区间1[,]iixx，0,1,,1in上均为x的三次多项式；2.()Sx在整个区间[,]ab上有直至二阶的连续导数。则称()Sx为[,]ab上依次划分的三次样条函数，简称样条函数。具体地有分段表达式：3200000132111112322222233211111,[,],[,](),[,](1),[,]nnnnnnaxbxcxdxxxaxbxcxdxxxSxaxbxcxdxxxaxbxcxdxxx共有4n个参数,,,,0,1,,iiiiabcdin，它们在内节点处满足00''00''''00()(),()(),1,2,,1.(2)()(),iiiiiiSxSxSxSxinSxSx满足样条函数定义的函数集合称为分划上的三次样条函数空间，记为(3,)S，可以证明(3,)S为线性空间。若()(3,)SxS，且进一步满足插值条件()(),0,1,,(3)iiiSxyfxin其中iy为节点ix处的给定函数值（若被插函数()fx已知，则用()ifx代替之），则称()Sx为以011,,,,nnxxxx为节点的三次样条函数。其中式（3）插值节点提供了1n个约束条件，加上式（2）的33n个，合起来共有42n个，欲求4n个待定参数的唯一解，尚缺两个条件。这两个条件一般由样条函数的边界条件提供。常用三类边界条件，他们分别与三次样条函数，构成不同边界条件的样条函数插值问题。2.三类样条函数插值问题2.1第二类边界条件给定边界条件两端的一阶导数值：''''000(),()nnnSxymSxym这相当于样条两短处的方向给定（压铁在两端点的压力方向确定），对应的插值问题如下：对于分划011:nnaxxxxb，给定节点对应的函数值012,,,,nyyyy，以及两端点处的一阶导数值'00ym，'nnym，求三次样条函数()Sx，使''00(),0,1,,(),()iinnSxyinSxmSxm2.2第一类边界条件给定边界两端的二阶导数值：''''''''000(),()nnnSxyMSxyM这相当于在样条两端处外加一个力矩，使梁两端点处有相应的曲率。对应的插值问题如下：对于划分011:nnaxxxxb，给定节点对应的函数值012,,,,nyyyy，以及两端点处的二阶导数值'00ym，'nnym，求三次样条函数()Sx，使''''00(),0,1,2,,(),()iinnSxyinSxMSxM特别地，若''''000,0nnyMyM，这相当于样条边界上不加力矩，样条在边界处是自由的，这样的样条称为自由样条，边界条件称为“自由边界条件”。2.3第三类边界条件被插函数()fx是以0nxx为周期的周期函数时，则要求()Sx也是周期函数，此时边界条件应满足：000''000''''000()(),()(),()(),nnnSxSxSxSxSxSx而且还要加上0nyy。这样得出的()Sx称为周期样条插值函数。3.三次样条插值问题的实际求解本部分分别编写了三次样条插值在第一类边界条件（二阶导数），第二类边界条件（一阶导数），第三类边界条件（周期函数）下求解的matlab程序，并给出了实际例子进行求解。3.1第一类边界条件的插值求解程序function[f,f0]=spline1(x,y,y_1,y_N,x0)symst;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y维数不相等');return;end%维数检查fori=1:nif(x(i)=x0)&&(x(i+1)=x0)index=i;break;endend%找到待求解x0值所在的区间A=diag(2*ones(1,n));%求解m的系数矩阵u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);fori=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+...3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);%形成系数矩阵及向量Cendc(1)=2*y_1;c(n)=2*y_N;m=followup(A,c);%用追赶法求解三对角方程组h=x(index+1)-x(index);%求x0所在区间长度f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+...y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+...m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-...m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h;%求x0所在区间的插值函数f1=simplify(f);f=collect(f1);%化简f并合并同类项。f0=subs(f,'t',x0);%x0处的插值end其中[f,f0]=spline1(x,y,y_1,y_N,x0)中的x,y分别是相应的插值向量，y_1和y_N为对应的边界条件值，x0为待求插值点。当x,y的维度n不是很大的时候，我们可以选择n-1个分别位于不同且无重合区域的区间段里的点Xi，分别插值，就可以得出整个区间里的三次样条插值函数。实例求解已知函数，()yfx，在若干点处的函数值如下表所示，求满足二阶导数别边界条件的样条插值函数，其中''''(0)0,(4)24SS。表4-1X和Y对应的数值X01234Y-8-701956解：分别取待插值节点为0.5,1.5,2.5,3.5进行四次运算，程序调用格式如下：以X=0.5处为例：调用格式为：[f,f0]=spline1(x,y,0,24,0.5)，运行结果截图如图4-1图4-1第一种边界条件区间[0,1]的插值结果即结果为：在区间[0,1]上的插值函数为231103()8(01)77tSxtx在x=0.5处的插值为-7.9286。同理在x=1.5,x=2.5,x=3.5处进行运算得相应区间的插值函数为图4-2第一种边界条件下区间[1,2]的插值结果图4-3第一种边界条件下区间[2,3]的插值结果图4-4第一种边界条件的区间[3,4]的插值结果综合以上四种计算结果，得出整个区间[0,4]上关于x的插值函数为：，232323231038,[0,1]772333644,[1,2]7777()40219468380,[2,3]7777116118537443822,[3,4]7777tttttttSttttttttt，3.2第二类边界条件插值的求解程序function[f,f0]=spline2(x,y,y2_1,y2_N,x0)symst;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y的维度不相等');return;endfori=1:nif(x(i)=x0)&&(x(i+1)=x0)index=i;break;endend%确定x0所在的区间A=diag(2*ones(1,n));%求解m的系数矩阵A(1,2)=1;A(n,n-1)=1;u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);fori=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+...3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);%形成系数矩阵及向量Cendc(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*y2_1/2;c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*y2_N/2;m=followup(A,c);%追赶法求解系数方程组h=x(index+1)-x(index);%求x0所在区间的长度f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+...y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+...m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-...m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h;%f1=simplify(f);f=collect(f1);%化简f以及进行合并同类项f0=subs(f,'t',x0);%求x0处的插值其中[f,f0]=spline2(x,y,y_1,y_N,x0)中的x,y分别是相应的插值向量，y_2和y2_N为对应的边界条件值，x0为待求插值点。当x,y的维度n不是很大的时候，我们可以选择n-1个分别位于不同且无重合区域的区间段里的点Xi，分别插值，就可以得出整个区间里的三次样条插值函数。实例求解构造适合下列数据表的三次插值样条x-1013y-1135'y6＼＼1解：其插值思路及程序调用同第一类插值，将区间分为三段，分别在x=-0.5，x=0.5，x=2处进行插值，得插值结果如下图4-5第二类边界条件在区间[-1,0]的插值结果图4-6第二类边界条件在区间[0,1]的插值结果图4-7第二类边界条件在区间[1,3]的插值结果综合图4-5，图4-6，图4-7的结果可知，基于表格数组的第二类边界插值结果为32323285161751,[1,0]6923691116175()1,[0,1]69236913723111,[1,3]276921292tttttttStttttt3.3第三类边界条件插值的求解程序function[f,f0]=spline3(x,y,x0)symst;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不想当！');return;end%维数检查fori=1:nif(x(i)=x0)&&(x(i+1)=x0)index=i;break;endend%找到x0所在的区间A=diag(2*ones(1,n-1));%形成系数求解m的系数矩阵h0=x(2)-x(1);h1=x(3)-x(2);he=x(n)-x(n-1);A(1,2)=h0/(h0+h1);A(1,n-1)=1-A(1,2);A(n-1,1)=he/(h0+he);A(n-1,n-2)=1-A(n-1,1);c=zeros(n-1,1);c(1)=3*A(1,n-1)*(y(2)-y(1))/h0+3*A(1,2)*(y(3)-y(2))/h1;fori=2:n-2u=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda=(x(i+1)-x(