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

《三次样条插值法《数值分析》上机实验作业.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

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为常数向量functionm

4、=followup(A,b)n=rank(A);fori=1:nifA(i,i)==0disp('error:对角元素中有数据为0');return;endendd=ones(n,1);a=ones(n-1,1);c=ones(n-1);fori=1:n-1a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1)=A(n,n);fori=2:nd(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(i-1,1);endm(n,1)

5、=b(n,1)/d(n,1);fori=(n-1):-1:1m(i,1)=(b(i,1)-c(i,1)*m(i+1,1))/d(i,1);end自定义函数：thrsample2.m%a为要求的插值点%f为区间内的插值函数%f0为输入点处的插值%m为追赶法解出的弯矩向量functionthrsample2(a)x=[0.2:0.2:1.0];y=[0.0.0.0.0.];s02=0;s10=0;x0=a;n=length(x);fori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endendA=diag(2*ones(1,n));A(1,2)

6、=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);endc(1)=3*(y(2)-y(1))/(x(2)-x(1))

7、-(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);symst;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