样条插值实验报告..docx

四、三次样条插值1.样条函数插值的原理给定区间上划分，若分段函数满足：在各个子区间，上均为的三次多项式；在整个区间上有直至二阶的连续导数。则称为上依次划分的三次样条函数，简称样条函数。具体地有分段表达式：共有个参数，它们在内节点处满足满足样条函数定义的函数集合称为分划上的三次样条函数空间，记为，可以证明为线性空间。若，且进一步满足插值条件其中为节点处的给定函数值（若被插函数已知，则用代替之），则称为以为节点的三次样条函数。其中式（3）插值节点提供了个约束条件，加上式（2）的个，合起来共有个，欲求个待定参数的唯一解，尚缺两个条件。这两个条件一般由样条函数的边界条件提供。常用三类边界条件，他们分别与三次样条函数，构成不同边界条件的样条函数插值问题。2.三类样条函数插值问题第二类边界条件给定边界条件两端的一阶导数值：这相当于样条两短处的方向给定（压铁在两端点的压力方向确定），对应的插值问题如下：对于分划，给定节点对应的函数值，以及两端点处的一阶导数值，，求三次样条函数，使第一类边界条件给定边界两端的二阶导数值：这相当于在样条两端处外加一个力矩，使梁两端点处有相应的曲率。对应的插值问题如下：对于划分，给定节点对应的函数值，以及两端点处的二阶导数值，，求三次样条函数，使特别地，若，这相当于样条边界上不加力矩，样条在边界处是自由的，这样的样条称为自由样条，边界条件称为“自由边界条件”。第三类边界条件被插函数是以为周期的周期函数时，则要求也是周期函数，此时边界条件应满足：而且还要加上。这样得出的称为周期样条插值函数。3.三次样条插值问题的实际求解本部分分别编写了三次样条插值在第一类边界条件（二阶导数），第二类边界条件（一阶导数），第三类边界条件（周期函数）下求解的matlab程序，并给出了实际例子进行求解。3.1第一类边界条件的插值求解程序function [f,f0] =spline1 (x,y,y_1, y_N,x0)syms t;f = 0.0;f0 = 0.0;if(length(x) == length(y)) n = length(x);else disp(x和y维数不相等);return;end%维数检查for i=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);for i=2:n-1 u(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，分别插值，就可以得出整个区间里的三次样条插值函数。实例求解已知函数，，在若干点处的函数值如下表