[2018年最新整理]三次样条插值法《数值分析》上机实验作业.docxVIP

昆明理工大学研究生《数值分析》上机实验作业姓名： 学号： 专业： 课题名称课题七 三次样条插值法1、问题提出设已知数据如下：0.20.40.60.81.009177710063860920.3843735求f(x)的三次样条插值函数S(x)。2、要求（1）满足自然边界条件；（2）满足第一类边界条件，。（3）打印输出用追赶法解出的弯曲向量(,,,)和 （i=0,1,2,3,4,5,6,7,8）的值。并画出的图形。班级、姓名、学号三、目的和意义由于航空、造船等工程设计的需要而发展起来所谓样条插值方法，既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性，而且具有较好的稳定性。今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越广泛地应用。其中，尤以三次样条插值函数应用最为广泛，如在高速飞机的机翼形体和船体放样等方面的应用，同时在计算机作图方面更是大有作为。它能够解决一些既有二阶光滑度，又有二阶连续导数的方程，具有良好的收敛性和稳定性。通过本次实验进一步了解三次样条插值函数，并通过求解三弯矩方程组得出曲线函数组；通过MATLAB编程实现求三次样条插值函数的算法，分别考虑不同的边界条件，同时用追赶法解出弯曲向量和 （i=0,1,2,3,4,5,6,7,8）的值。四、计算公式首先我们利用的二阶导数值表达,因为在区间上是不高于三次的多项式，其二阶导数必是线性函数，所以可表示为： 对积分两次并利用，可定出积分常数，于是得三次样条表达式。 这里是未知的。为了确定，对求导得 由此可得类似的可求出在区间上的表达式，进而得 利用可得 ，（三弯矩方程）其中 ，， .其中有（）个未知数，而方程只有（n-1）个，当满足第一种边界条件时，可得另两个方程，如果令，将上述方程综合后的一下矩阵形式： 可以证明此方程组满足追赶法的条件，我们用追赶法可得的值，将其带入公式即得。对第二种边界条件，直接的端点方程并且令，则又得三弯矩方程同理即可求得解。五、结构程序设计1.满足自然边界条件时 自定义函数：followup.m %追赶法求m %A为线性方程组的系数矩阵 %b为常数向量 function 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(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.9798652 0.9177710 0.8080348 0.6386092 0.3843735]; s02=0; s10=0; x0=a;n=length(x);for i=1:n if (x(i)=x0)(x(i+1)=x0) index=i; break; end end A=diag(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(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