样条插值实验报告

1、四、三次样条插值1.样条函数插值的原理给定区间上划分，若分段函数满足：1. 在各个子区间，上均为的三次多项式；2. 在整个区间上有直至二阶的连续导数。则称为上依次划分的三次样条函数，简称样条函数。具体地有分段表达式：共有个参数，它们在内节点处满足满足样条函数定义的函数集合称为分划上的三次样条函数空间，记为，可以证明为线性空间。若，且进一步满足插值条件其中为节点处的给定函数值（若被插函数已知，则用代替之），则称为以为节点的三次样条函数。其中式（3）插值节点提供了个约束条件，加上式（2）的个，合起来共有个，欲求个待定参数的唯一解，尚缺两个条件。这两个条件一般由样条函数的边界条件提供。常用三类边界条

2、件，他们分别与三次样条函数，构成不同边界条件的样条函数插值问题。2.三类样条函数插值问题2.1 第二类边界条件给定边界条件两端的一阶导数值：这相当于样条两短处的方向给定（压铁在两端点的压力方向确定），对应的插值问题如下：对于分划，给定节点对应的函数值，以及两端点处的一阶导数值，求三次样条函数，使2.2 第一类边界条件给定边界两端的二阶导数值：这相当于在样条两端处外加一个力矩，使梁两端点处有相应的曲率。对应的插值问题如下：对于划分，给定节点对应的函数值，以及两端点处的二阶导数值，求三次样条函数，使特别地，若，这相当于样条边界上不加力矩，样条在边界处是自由的，这样的样条称为自由样条，边界条件称为“

3、自由边界条件”。2.3 第三类边界条件被插函数是以为周期的周期函数时，则要求也是周期函数，此时边界条件应满足：而且还要加上。这样得出的称为周期样条插值函数。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

4、 disp(x和y维数不相等); return;end %维数检查 for i=1:n if(x(i)=x0) index = i; break; end end %找到待求解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)

5、/(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(in

6、dex)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个分别位于不同且无重

7、合区域的区间段里的点Xi，分别插值，就可以得出整个区间里的三次样条插值函数。 实例求解已知函数，在若干点处的函数值如下表所示，求满足二阶导数别边界条件的样条插值函数，其中。表4-1 X和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上的插值函数为在x=0.5处的插值为-7.9286。同理在x=1.5,x=2.5,x=3.5处进行运算得相

8、应区间的插值函数为图4-2 第一种边界条件下区间1,2的插值结果图4-3 第一种边界条件下区间2,3的插值结果图4-4 第一种边界条件的区间3,4的插值结果综合以上四种计算结果，得出整个区间0,4上关于x的插值函数为：，3.2 第二类边界条件插值的求解程序function f,f0 = spline2 (x,y,y2_1, y2_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:n if(x(i)=x0) index = i;

9、 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);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)

10、-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

11、)*(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不是很大的时候，我

12、们可以选择n-1个分别位于不同且无重合区域的区间段里的点Xi，分别插值，就可以得出整个区间里的三次样条插值函数。实例求解构造适合下列数据表的三次插值样条-1013-113561解：其插值思路及程序调用同第一类插值，将区间分为三段，分别在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的结果可知，基于表格数组的第二类边界插值结果为3.3 第三类边界条件插值的求解程序function f,f0 = spline3 (x,y

13、,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:n if(x(i)=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)

14、 = 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;for i=2:n-2 u = (x(i)-x(i-1)/(x(i+1)-x(i-1); lamda = (x(i+1)-x(i)/(x(i+1)-x(i-1); c(i) = 3*lamda*(y(i)-y(i-1)/(x(i)-x(i-1)+ . 3*u*(y(i+1)-y(i)/(x(i+1)-x(i); A(i, i+1) = u; A(i, i-1) =

15、lamda; %形成系数矩阵及向量Cendc(n-1) = 3*( he*(y(2)-y(1)/h0+h0*( y(n)-y(n-1)/he)/(h0+he);m = zeros(n,1);m(2:n) = qrxq(A,c); %用QR分解法求解方程m(1) = m(n);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);f0 = subs(f,t,x0); %x0处的插值其中f,f0 =spline2 (x,y,x0)中的x,y分别是相应的插值向量， x0为待求插值点。当x,y的维度n不是很大的时候，我们可以选择n-1个分别位于不同且无重合区域的区间段里的点Xi，分别插值，就可以得出整个区间里的三次样条插值函数。实例求解函数f(x)的数值表如下：x2