《三次样条插值法《数值分析》上机实验作业》由会员分享，可在线阅读，更多相关《三次样条插值法《数值分析》上机实验作业（13页珍藏版）》请在金锄头文库上搜索。
1、昆明理工大学研究生数值分析上机实验作业姓名：学号：专业：、课题名称课题七三次样条插值法1、问题提出设已知数据如下：0.20.40.60.81.00.97986520.91777100.80803480.63860920.3843735求f(x)的三次样条插值函数S(x)。2、要求(1) 满足自然边界条件s(0.2)s(1.0)0;(2) 满足第一类边界条件s(0.2)0.20271,s(1.0)=1.55741。(3)打印输出用追赶法解出的弯曲向量(M,M,M,M,M)和01234S(0.2,O.li)(i=0,1,2,3,4,5,6,7,8的值。并画出y=S(x)的图形。二、班级、姓名、学号
2、三、目的和意义由于航空、造船等工程设计的需要而发展起来所谓样条插值方法既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性，而且具有较好的稳定性。今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越广泛地应用。其中，尤以三次样条插值函数应用最为广泛，如在高速飞机的机翼形体和船体放样等方面的应用，同时在计算机作图方面更是大有作为。它能够解决一些既有二阶光滑度，又有二阶连续导数的方程，具有良好的收敛性和稳定性。1.通过本次实验进一步了解三次样条插值函数，并通过求解三弯矩方程组得出曲线函数组；2.通过MATLAB编程实现求三次样条插值函数的算法，分别考虑不同的边界条
3、件，同时用追赶法解出弯曲向量和S(0.2O.li)(i=0,l,2,3,4,5,6,7,8)的值。四、计算公式首先我们利用S(x)的二阶导数值S，(x)二M(j二0,1,2,n)表jj达S(x),因为在区间xj，Xj1上S(x)二S(x)是不高于三次的多项式,j其二阶导数S(x)必是线性函数，所以可表示为：”x一xx一xS(x)Mj+1+Mj,xgx,xjhj+1hjj+1jj对S(x)积分两次并利用S(x)y,S(x)y，可定出积分j+1j常数，于是得三次样条表达式。(X一X)3S(x)二Mj1Mj6hj+1jMx一x+(y一/1h2)jj+16jhj,j0,厶2,n一1(x一x)36hj(
4、yj2)X一Xj1hj-j+1这里M(j0,1,n)是未知的。为了确定M(j0,1,，n)，对S(x)jj求导得(x一x)2(x一x)2y一yhS，(x)Mj+1+Mj+j+1jj(MM)j2hj+12hh6j+1jjjj由此可得hhS(x,0)二j_MM6j+13jhj-iM3j类似的可求出S(x)在区间xj_i，xj上的表达式，进而得y-yhS(x0)=j+亠1Mjh6j-1j-1利用S(x-0)二S(x+0)可得jj“M+2M+久M二djj-1jjj+1jj=1,n-1,(三弯矩方程)hh+hj-1jj-1jj，1丄八1，-1,6nhnn-1nn-1其中“二j，久jh+hjj-1jfx,
5、x-fx,xd二6j4i-Ljjh+hj-1j其中有(n+1)个未知数M,M，M，而方程只有(nT)个,01n当满足第一种边界条件时，可得另两个方程2M,M=A(flx,x广)，01h0100M,2M=(广-fx,x)n-1nhnn-1nn-1如果令九=1,d=(ftx,x-广)卩=1,g=(广-ftx,x)，将00h0100上述方程综合后的一下矩阵形式：2X-M_d-000卩2XMd11112Mdn-1n-1n-1n-1卩2M_dnnn可以证明此方程组满足追赶法的条件，我们用追赶法可得M的(n+1)值，将其带入公式即得s(x)。对第二种边界条件，直接的端点方程M=f,M=f并且令00nn九二
6、,二0,d二2fd二2f，则又得三弯矩方程同理即可求得解。0n00nn五、结构程序设计1. 满足自然边界条件时自定义函数：fOllowup.m%追赶法求m%A为线性方程组的系数矩阵%b为常数向量functionm=f0llowup(A,b)n=rank(A);fOri=l:nifA(i,i)=0disp(error:对角元素中有数据为O);return;endendd=ones(n,l);a=ones(n-l,l);c=ones(n-l);fori=l:n-la(i,l)=A(i+l,i);c(i,l)=A(i,i+l);d(i,l)=A(i,i);endd(n,l)=A(n,n);fori=
7、2:nd(i,l)=d(i,l)-(a(i-l,l)/d(i-l,l)*c(i-l,l);b(i,l)=b(i,l)-(a(i-l,l)/d(i-l,l)*b(i-l,l);endm(n,l)=b(n,l)/d(n,l);fori=(n-l):-l:lm(i,l)=(b(i,l)-c(i,l)*m(i+l,l)/d(i,l);end自定义函数：thrsample2.m%a为要求的插值点%f为区间内的插值函数%fO为输入点处的插值%m为追赶法解出的弯矩向量functionthrsample2(a)x=O.2:O.2:l.O;y=O.9798652O.9l777lOO.8O8O348O.6386O
8、92O.3843735;s02=0;slO=O;xO=a;n=length(x);fori=1:nif(x(i)=x0)index=i;break;endendA=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);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)*(
9、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)-(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)A2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index)A2/
10、h/h/h+m(index)*(t-x(index)*(x(index+1)-t)A2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index)A2/h/hf0=subs(f,t,x0)运了-C;UsusXDe-s1&21O3001rs-ampll-e-Zi.mfolksfkirsarripJem运行结果(方程S、弯矩M和插值函数f的值)为：50 =0.979865200051 =0.952572670352 =0.917771000053 =0.869504939054 =0.8080348000d2iM4Q.674755 =0.733408548556 =0.63
11、86092000-It57 =0.518730429658 =0.3843735000f=125*(8825841099186633*t)/4503599627370496-8825841099186633/45035996273704960)*(t-2/5)人2-125*(8266546267222895*t)/4503599627370496-8266546267222895/9007199254740992)*(t-1/5)A2-25*(7396583675403915*t)/18014398509481984-1479316735080783/9007199254740992)*(t-
12、1/5)人2-25*(t-2/5)A2*(2290591133669*t)/8796093022208-2290591133669/43980465111040)2. 满足第一类边界条件时自定义函数：thrsamplel.m%a为要求的插值点%f为区间内的插值函数%f0为输入点处的插值%m为追赶法解出的弯矩向量functionthrsamplel(a)x=0.2:0.2:l.0;y=0.97986520.9l777l00.80803480.63860920.3843735;s02=0.2027l;sl0=l.5574l;x0=a;n=length(x);fori=l:nif(x(i)=x0)i
13、ndex=i;break;endendA=diag(2*ones(l,n);u=zeros(n-2,l);lamda=zeros(n-l,l);c=zeros(n,l);fori=2:n-lu(i-l)=(x(i)-x(i-l)/(x(i+l)-x(i-l);lamda(i)=(x(i+l)-x(i)/(x(i+l)-x(i-l);c(i)=3*lamda(i)*(y(i)-y(i-l)/(x(i)-x(i-l)+3*u(i-l)*(y(i+l)-y(i)/(x(i+l)-x(i);A(i,i+l)=u(i-l);A(i,i-l)=lamda(i);endc(l)=2*s02;c(n)=2*sl0;m=followup(A,c)h=x(index+l)-x(index);symst;f=y(index)*(2*(t-x(index)+h)*(t-x(index+l)A2/h/h/h+y(index+l)*(2*(x(index+l)-t)+h)*(t-x(index)A2/h/h/h+m(index)*(t-x(index)*(x(index+l)-t)A2/h/h-m(index+l)*(x(index+l)-t)*(t-x(index)A2/h/hf0=subs(f,t,x0)运行结果(方程S、弯矩M和插值函数f的值)为:50 =0.979865200051 =0.9
