三次样条插值法数值分析上机试验作业

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°, M1, M2 , M3, M4)和S(0.2+0.1i) (i=0,1,2,3,4,5,6,7,8)的值

2、。并画出 y = S(x)的图形。班级、姓名、学号三、目的和意义由于航空、造船等工程设计的需要而发展起来所谓样条插值方法， 既保留了分段低次插值多项式的各种优点， 又提高了插值函数的光滑 性，而且具有较好的稳定性。今天，样条插值方法已成为数值逼近的 一个极其重要的分支，在许多领域里得到越来越广泛地应用。其中， 尤以三次样条插值函数应用最为广泛，如在高速飞机的机翼形体和船 体放样等方面的应用，同时在计算机作图方面更是大有作为。 它能够 解决一些既有二阶光滑度，又有二阶连续导数的方程，具有良好的收敛性和稳定性。1,通过本次实验进一步了解三次样条插值函数，并通过求解三弯矩方程组得出曲线函数组；2.通

3、过MATLA褊程实现求三次样条插值函数的算法，分别考虑不同的边界条件，同时用追赶法解出弯曲向量和S(0.2 + 0.1i)(i=0,1,2,3,4,5,6,7,8)的值。四、计算公式首先我们利用S(x)的二阶导数值STXj) = M(j = 0,1,2,，n)表 达S(x)，因为在区间xj, Xj+1上s(x) = Sj(x)是不高于三次的多项式， 其二阶导数S(x)必是线性函数，所以可表示为：x； , xx - x；S"(x) = Mjj + M j+L , x w x j,为中对S(x)积分两次并j hjj hj利用S (xj ) = yj, S (xj书)=y书，可定出积分常数

4、，于是得三次样条表达式M j 2(yj%6)x±_L_ hj,、3,、3(xj 1 - x ) 一 (x - xj )S (x) = M j M j 1 6hj6%gh6(yj 12 X - xj j ) L hjj = 0,1,2,这里M(j=0,1,n)是未知的为了确定M(j=0,1，n),对S(x)求导得(xj 1 - x)2S (x) - -M -j 2hj(x - xj )2 j - y2hjhjhj- U (M j 1 - M j )由此可得S (Xj 0) = yj Lyjjhhjhj类似的可求出Stx）在区间Xj-i,Xj上的表达式，进而得hjjyj -Yj JhjR

5、,S (Xj - 0) -'M j Jhj.16禾I用S (Xj 0)SXj +0）可得gMj+2Mj + Mj 书=dj j = 1,，n 1,（三弯矩方程）其中j?hj j hjf Xj，xj i - f Xj，xjdj = 6 二 6f Xj/，Xj, Xj i L j = 1, n - 1,hj _ihj其中有（n+1）个未知数M0,M1,.Mn,而方程只有（n-1 ）个,当满足第一种边界条件时，可得另两个方程2M0 +M1 =" Lx 】- f0), h06Mn12Mn fn - f Xn，4 1 hn 1如果令-=1,d0 =g(f k0,X1 L f0)Pn =

6、1,gn =-(f； f【Xn,Xn ),将 h0hn.1上述方程综合后的一下矩阵形式：2 b11 M0 1 j d0 I H 2 入| M1 I &I XI ： =l：l-2"1|Mn_1|dn_11n 2 Ji Mn J dn J可以证明此方程组满足追赶法的条件，我们用追赶法可得（n+1）值，将其带入公式即得s（X）对第二种边界条件，直接的端点方程M o = f0：M n = f；并且令,0 =，=0,或=2f°；dn =2f；;则又得三弯矩方程同理即可求得解。五、结构程序设计1 .满足自然边界条件时自定义函数：followup.m%追赶法求m%A为线性方程组的

7、系数矩阵%b为常数向量function m=followup(A,b)n=rank(A);fori=1:nif A(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

8、,1)*b(i-1,1);endm(n,1)=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为追赶法解出的弯矩向量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);fori=1:nif (x(i)<=x0)&

9、;&(x(i+1)>=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)*(y(i+1)-y(i)/(x(i+1)-x(i);A(i,i+1)=u(i-1

10、);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);syms t;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(ind ex)A2/h/h/h+m(index)*(t-x(index)*(x(index+1

11、)-t)A2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index) A2/h/hf0=subs(f,'t',x0)运行结果(方程S0 = 0.9798652000S1 = 0.9525726703S2 = 0.9177710000S3 = 0.8695049390S4 = 0.8080348000S5 = 0.7334085485S6 = 0.6386092000S7 = 0.5187304296S8 =0.3843735000S、弯矩M和插值函数f的值)为：rsnsuliciho62"103GO 1 1吕r5r口 mpleZ.m fol

12、lowup, mTh rsa mpleZ. m 产 +市定义的釉或塞里，>> tlur=3JHple2CO. 2)-O- 2604-o,-1_ OTB1-1. 3677f =125*(8825841099186633*t)/4503599627370496 - 8825841099186633/45035996273704960)*(t - 2/5)A2 - 125*(8266546267222895*t)/4503599627370496 - 8266546267222895/9007199254740992)*(t - 1/5)人2 - 25*(7396583675403915

13、*t)/18014398509481984 - 1479316735080783/9007199254740992)*(t - 1/5)A2 - 25*(t - 2/5)A2*(2290591133669*t)/8796093022208 - 2290591133669/43980465111040)2 .满足第一类边界条件时自定义函数：thrsample1.m%a为要求的插值点%f为区间内的插值函数%f0为输入点处的插值%m为追赶法解出的弯矩向量function thrsample1(a) x=0.2:0.2:1.0;y=0.9798652 0.9177710 0.8080348 0.638

14、6092 0.3843735;s02=0.20271;s10=1.55741;x0=a;n=length(x);fori=1:nif (x(i)<=x0)&&(x(i+1)>=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);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(

15、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)=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-x(index+1)A2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index)A2/h/h/h+m(index)*(t-x(index)*(x(index+1)-t)A2/h/

16、h-m(index+1)*(x(index+1)-t)*(t-x(index)A2/h/hf0=subs(f,'t',x0)运行结果(方程S0 = 0.9798652000S1 = 0.9685577748S2 = 0.9177710000S3 = 0.8590474261S4 = 0.8080348000S5 = 0.7592534958S6 = 0.6386092000S7 = 0.4258081535S8 =0.3843735000S、弯矩M和插值函数f的值)为：0- OU0162103001amplei. mfollowupn thrwmpi»1jn x 1

17、+*专行1口» thrsaaplel (O.Q. 2027-0. 5669-0,4327-la 8S991. 5674f =125*(8825841099186633笛/4503599627370496 - 8825841099186633/45035996273704960)*(t - 2/5)A2 - 125*(8266546267222895笛/4503599627370496 - 8266546267222895/9007199254740992)*(t - 1/5)八2 +25*(20271*t)/100000 - 20271/500000)*(t - 2/5)八2 - 2

18、5*(1321529499150801*t)/2251799813685248 -1321529499150801/5629499534213120)*(t - 1/5)八23.画出y=S(x)的图形(1)满足自然边界条件时程序为：x=0.2:0.2:1;y=0.9798652 0.9177710 0.8080348 0.6386093 0.3843735;xi=0.2:0.01:1.0;yi=interp1(x,y,xi,'variational');plot(x,y,'o',xi,yi,'k-');输出图形为：图1自然边界条件下的y=S(x)

19、的图形(2)满足第一类边界条件时 程序为：x=0.2 0.4 0.6 0.8 1.0 ;y=0.9798652 0.9177710 0.8080348 0.6386093 0.3843735;pp=csape(x,y,'complete',0.20271,1.55741);%complete代表一次插值xi=0.2:0.01:1.0;yi=ppval(pp,xi);pp.coefsplot(x,y,'o',xi,yi);输出图形为：Figure 1 X文件旧 病电旧 者看俚 ffiACD Tfll0 嗯面 窗口故)海勖(HJ 013 冷4， 巧耍*一0口国口- r I QlB rOlT 卜0.B r也5 -44也 3«1110.20.30.40.5 Q.S 0.70.8091图2第一类边界条件下的 y=S(x)的图形六、结果讨论和分析在插值法中，拉格朗日插值法的公式结构整齐紧凑，在理论分析中十分方便，然而在计算中，当插值点增加或减少一个时，所对应的 基本多项