三次样条插值法《数值分析》上机实验作业

1、昆明理工大学研究生数值分析上机实验作业姓名： 学号： 专业： 一、 课题名称课题七 三次样条插值法1、问题提出设已知数据如下：0.20.40.60.81.00.0.0.0.0.求f(x)的三次样条插值函数S(x)。2、要求（1）满足自然边界条件；（2）满足第一类边界条件，。（3）打印输出用追赶法解出的弯曲向量(,)和 （i=0,1,2,3,4,5,6,7,8）的值。并画出的图形。二、 班级、姓名、学号三、目的和意义由于航空、造船等工程设计的需要而发展起来所谓样条插值方法，既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性，而且具有较好的稳定性。今天，样条插值方法已成为数值逼近的一个

2、极其重要的分支，在许多领域里得到越来越广泛地应用。其中，尤以三次样条插值函数应用最为广泛，如在高速飞机的机翼形体和船体放样等方面的应用，同时在计算机作图方面更是大有作为。它能够解决一些既有二阶光滑度，又有二阶连续导数的方程，具有良好的收敛性和稳定性。1. 通过本次实验进一步了解三次样条插值函数，并通过求解三弯矩方程组得出曲线函数组；2. 通过MATLAB编程实现求三次样条插值函数的算法，分别考虑不同的边界条件，同时用追赶法解出弯曲向量和 （i=0,1,2,3,4,5,6,7,8）的值。四、计算公式首先我们利用的二阶导数值表达,因为在区间上是不高于三次的多项式，其二阶导数必是线性函数，所以可表示

3、为： 对积分两次并利用，可定出积分常数，于是得三次样条表达式。 这里是未知的。为了确定，对求导得 由此可得类似的可求出在区间上的表达式，进而得 利用可得 ，（三弯矩方程）其中 ， .其中有（）个未知数，而方程只有（n-1）个，当满足第一种边界条件时，可得另两个方程，如果令，将上述方程综合后的一下矩阵形式： 可以证明此方程组满足追赶法的条件，我们用追赶法可得的值，将其带入公式即得。对第二种边界条件，直接的端点方程并且令，则又得三弯矩方程同理即可求得解。五、结构程序设计1.满足自然边界条件时 自定义函数：followup.m %追赶法求m %A为线性方程组的系数矩阵 %b为常数向量 functio

4、n m=followup(A,b) n=rank(A); for i=1:n if A(i,i)=0 disp(error:对角元素中有数据为0); return; end end d=ones(n,1); a=ones(n-1,1); c=ones(n-1); for i=1:n-1 a(i,1)=A(i+1,i); c(i,1)=A(i,i+1); d(i,1)=A(i,i); end d(n,1)=A(n,n); for i=2:n d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1)*c(i-1,1); b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1)*b(

5、i-1,1); end m(n,1)=b(n,1)/d(n,1); for i=(n-1):-1:1 m(i,1)=(b(i,1)-c(i,1)*m(i+1,1)/d(i,1); end 自定义函数：thrsample2.m %a为要求的插值点 %f为区间内的插值函数 %f0为输入点处的插值 %m为追赶法解出的弯矩向量 function thrsample2(a) x=0.2:0.2:1.0; y=0. 0. 0. 0. 0.; s02=0; s10=0; x0=a;n=length(x);for i=1:n if (x(i)=x0) index=i; break; end end A=dia

6、g(2*ones(1,n); 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)-x(i); A(i,i+1)=u(i-1); A(i,i-1)=lamda(i); end c(1)=3*(y(

7、2)-y(1)/(x(2)-x(1)-(x(2)-x(1)*s02/2;c(n)=3*(y(n)-y(n-1)/(x(n)-x(n-1)-(x(n)-x(n-1)*s10/2; m=followup(A,c) h=x(index+1)-x(index); syms t; 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)

8、-t)*(t-x(index)2/h/h f0=subs(f,t,x0)运行结果（方程S、弯矩M和插值函数f的值）为：S0 = 0. S1 = 0. S2 = 0. S3 = 0.S4 = 0. S5 = 0. S6 = 0. S7 = 0.S8 = 0.f =125*(86633*t)/70496 - 86633/)*(t - 2/5)2 - 125*(22895*t)/70496 - 22895/40992)*(t - 1/5)2 - 25*(03915*t)/ - 80783/40992)*(t - 1/5)2 - 25*(t - 2/5)2*(69*t)/08 - 69/040)2.满

9、足第一类边界条件时 自定义函数：thrsample1.m %a为要求的插值点 %f为区间内的插值函数 %f0为输入点处的插值 %m为追赶法解出的弯矩向量 function thrsample1(a) x=0.2:0.2:1.0; y=0. 0. 0. 0. 0.; s02=0.20271; s10=1.55741; x0=a; n=length(x); for i=1:n if (x(i)=x0) index=i; break; end end A=diag(2*ones(1,n); u=zeros(n-2,1); lamda=zeros(n-1,1); c=zeros(n,1); for i

10、=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); end c(1)=2*s02; c(n)=2*s10;m=followup(A,c) h=x(index+1)-x(index); syms t; f=y(index)*(2*(t-x(index)+h)*(t

11、-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 f0=subs(f,t,x0)运行结果（方程S、弯矩M和插值函数f的值）为：S0 = 0. S1 = 0. S2 = 0. S3 = 0.S4 = 0. S5 = 0. S6 = 0. S7 = 0.S8 = 0.f =125*(86633*t)/70496 - 86633/)*(t - 2

12、/5)2 - 125*(22895*t)/70496 - 22895/40992)*(t - 1/5)2 + 25*(20271*t)/ - 20271/)*(t - 2/5)2 - 25*(50801*t)/85248 - 50801/13120)*(t - 1/5)23.画出y=S(x)的图形（1）满足自然边界条件时程序为：x=0.2:0.2:1; y=0. 0. 0. 0. 0.; xi=0.2:0.01:1.0; yi=interp1(x,y,xi,variational); plot(x,y,o,xi,yi,k-);输出图形为：图1 自然边界条件下的y=S(x)的图形（2）满足第一类

13、边界条件时程序为：x=0.2 0.4 0.6 0.8 1.0 ; y=0. 0. 0. 0. 0.; pp=csape(x,y,complete,0.20271,1.55741);%complete代表一次插值 xi=0.2:0.01:1.0;yi=ppval(pp,xi); pp.coefs plot(x,y,o,xi,yi);输出图形为：图2 第一类边界条件下的y=S(x)的图形六、结果讨论和分析在插值法中，拉格朗日插值法的公式结构整齐紧凑，在理论分析中十分方便，然而在计算中，当插值点增加或减少一个时，所对应的基本多项式就需要全部重新计算，于是整个公式都会变化，非常繁琐。此外，当插值点比较

14、多的时候，拉格朗日插值多项式的次数可能会很高，因此具有数值不稳定的特点，也就是说尽管在已知的几个点取到给定的数值，但在附近却会和“实际上”的值之间有很大的偏差，这类现象也被称为龙格现象。分段线性插值，逼近程度虽然好，但光滑性差，即不能保证节点处插值函数的导数连续，从而不能满足某些工程技术上的要求。分段三次Hermite插值，逼近程度好，光滑性也有所提高，但也增加了更多的条件，不太实用。而三次样条插值多项式，结合了二者的优点，即逼近程度好，光滑性强，不需要增加太多的条件，很实用。实验过程中遇到了许多问题，首先就是MATLAB基础不太扎实，编程过程比较繁琐，程序比较复杂，在查阅大量资料并询问同学后，最终得以解决。其次，在编写满足第一类边界条件下的y=S(x)图形时，没有注意到interp1语句无法使用边界条件,致使编写出的程序出现错误，最