三次樣條插值法數(shù)值分析上機試驗作業(yè)

1、昆明理工大學(xué)研究生數(shù)值分析上機實驗作業(yè)姓名:學(xué)號:專業(yè):課題名稱課題七三次樣條插值法1、問題提出設(shè)已知數(shù)據(jù)如下:0.20.40.60.81.00.97986520.91777100.80803480.63860920.3843735求f(x)的三次樣條插值函數(shù)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)的圖形。班級、姓名、學(xué)號三、目的和意義由于航空、造船等工程設(shè)計的需要而發(fā)展起來所謂樣條插值方法， 既保留了分段低次插值多項式的各種優(yōu)點， 又提高了插值函數(shù)的光滑 性，而且具有較好的穩(wěn)定性。今天，樣條插值方法已成為數(shù)值逼近的 一個極其重要的分支，在許多領(lǐng)域里得到越來越廣泛地應(yīng)用。其中， 尤以三次樣條插值函數(shù)應(yīng)用最為廣泛，如在高速飛機的機翼形體和船 體放樣等方面的應(yīng)用，同時在計算機作圖方面更是大有作為。 它能夠 解決一些既有二階光滑度，又有二階連續(xù)導(dǎo)數(shù)的方程，具有良好的收斂性和穩(wěn)定性。1,通過本次實驗進一步了解三次樣條插值函數(shù)，并通過求解三彎矩方程組得出曲線函數(shù)組；2.通

3、過MATLA褊程實現(xiàn)求三次樣條插值函數(shù)的算法，分別考慮不同的邊界條件，同時用追趕法解出彎曲向量和S(0.2 + 0.1i)(i=0,1,2,3,4,5,6,7,8)的值。四、計算公式首先我們利用S(x)的二階導(dǎo)數(shù)值STXj) = M(j = 0,1,2,，n)表 達S(x)，因為在區(qū)間xj, Xj+1上s(x) = Sj(x)是不高于三次的多項式， 其二階導(dǎo)數(shù)S(x)必是線性函數(shù)，所以可表示為：x； , xx - x；S"(x) = Mjj + M j+L , x w x j,為中對S(x)積分兩次并j hjj hj利用S (xj ) = yj, S (xj書)=y書，可定出積分常數(shù)

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)求導(dǎo)得(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）在區(qū)間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）個未知數(shù)M0,M1,.Mn,而方程只有（n-1 ）個,當(dāng)滿足第一種邊界條件時，可得另兩個方程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；;則又得三彎矩方程同理即可求得解。五、結(jié)構(gòu)程序設(shè)計1 .滿足自然邊界條件時自定義函數(shù)：followup.m%追趕法求m%A為線性方程組的

7、系數(shù)矩陣%b為常數(shù)向量function m=followup(A,b)n=rank(A);fori=1:nif A(i,i)=0disp('error:對角元素中有數(shù)據(jù)為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自定義函數(shù)：thrsample2.m%a為要求的插值點%f為區(qū)間內(nèi)的插值函數(shù)%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)運行結(jié)果(方程S0 = 0.9798652000S1 = 0.9525726703S2 = 0.9177710000S3 = 0.8695049390S4 = 0.8080348000S5 = 0.7334085485S6 = 0.6386092000S7 = 0.5187304296S8 =0.3843735000S、彎矩M和插值函數(shù)f的值)為：rsnsuliciho62"103GO 1 1呂r5r口 mpleZ.m fol

12、lowup, mTh rsa mpleZ. m 產(chǎn) +市定義的釉或塞里，>> 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 .滿足第一類邊界條件時自定義函數(shù)：thrsample1.m%a為要求的插值點%f為區(qū)間內(nèi)的插值函數(shù)%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)運行結(jié)果(方程S0 = 0.9798652000S1 = 0.9685577748S2 = 0.9177710000S3 = 0.8590474261S4 = 0.8080348000S5 = 0.7592534958S6 = 0.6386092000S7 = 0.4258081535S8 =0.3843735000S、彎矩M和插值函數(shù)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)的圖形六、結(jié)果討論和分析在插值法中，拉格朗日插值法的公式結(jié)構(gòu)整齊緊湊，在理論分析中十分方便，然而在計算中，當(dāng)插值點增加或減少一個時，所對應(yīng)的 基本多項