完整版數(shù)值分析作業(yè)-三次樣條插值

1、數(shù)值計(jì)算方法作業(yè)實(shí)驗(yàn)名稱實(shí)驗(yàn)4.3三次樣條插值函數(shù) (P126)4.5三次樣條插值函數(shù)的收斂性 (P127)實(shí)驗(yàn)時(shí)間姓名班級(jí)學(xué)號(hào)成績(jī)實(shí)驗(yàn)4.3三次樣條差值函投r實(shí)驗(yàn)?zāi)康模赫莆杖螛訔l插值函數(shù)的三彎矩方法實(shí)驗(yàn)函數(shù)：1 x工f (x)e 2dt2二二x0.00.10.20.30.4F(x)10.50000.53980.579310.61790.7554求f(0.13)和f(0.36)的近似值實(shí)驗(yàn)內(nèi)容：(1)編程實(shí)現(xiàn)求三次樣條插值函數(shù)的算法,分別考慮不同的邊界條件；(2)計(jì)算各插值節(jié)點(diǎn)的彎矩值；(3)在同一坐標(biāo)系中繪制函數(shù)f(x),插值多項(xiàng)式,三次樣條插值多項(xiàng)式的曲線 比較插值結(jié)果.實(shí)驗(yàn)4.5三次樣

2、條差值函數(shù)的收斂性實(shí)驗(yàn)?zāi)康模憾囗?xiàng)式插值不一定是收斂的,即插值的節(jié)點(diǎn)多,效果不一定好.對(duì)三次樣條插值 函數(shù)如何呢理論上證實(shí)三次樣條插值函數(shù)的收斂性是比較困難的, 通過(guò)本實(shí)驗(yàn) 可以證實(shí)這一理論結(jié)果.實(shí)驗(yàn)內(nèi)容：根據(jù)一定的規(guī)那么分別選擇等距或非等距的插值節(jié)點(diǎn),并不斷增加插值節(jié)點(diǎn)的個(gè)數(shù).實(shí)驗(yàn)要求：(1)隨著節(jié)點(diǎn)個(gè)數(shù)的增加,比較被逼近函數(shù)和三樣條插值函數(shù)的誤差變化情 況,分析所得結(jié)果并與拉格朗日插值多項(xiàng)式比較；(2)三次樣條插值函數(shù)的思想最早產(chǎn)生于工業(yè)部門(mén). 作為工業(yè)應(yīng)用的例子,考Xk0i23456789ioyk0.00.79i.532.i92.7i3.033.272.893.063.i93.29yk0.

3、80.2慮如下例子：某汽車(chē)制造商根據(jù)三次樣條插值函數(shù)設(shè)計(jì)車(chē)門(mén)曲線, 其中一 段數(shù)據(jù)如下：算法描述:拉格朗日插值：錯(cuò)誤！未找到引用源.n其中錯(cuò)誤！未找到引用源.是拉格朗日基函數(shù),其表達(dá)式為：ii(x)=nj=0HiX - Xj(xi - Xj )牛頓插值:Nn(X)f(Xo)fXo,Xi(X -X0)fXo,Xi,X2(X - xO)(x - Xi) .fXo,Xi,X2,.Xn(X -Xo)(X -Xi).(X -Xnj)f(Xi)-f(Xj) f Xi, Xj :Xi -XjfXj,Xk - fx ,Xj f Xi , Xj ,Xk:Xk -Xi其中fXo,Xi.Xn =(fXi,X2,.X

4、j 一 fXo,X,.Xn.)/(Xn - %)三樣條插值：所謂三次樣條插值多項(xiàng)式Sn(x)是一種分段函數(shù),它在節(jié)點(diǎn)Xi(a<X0<Xi <Xn<b)分成的每個(gè)小區(qū)間xi-i,xi上是三次多項(xiàng)式,其在此區(qū)間上的表達(dá)式如下：S( x) =M/、3(Xi -X)6hi(x.yi - yi4 hi(MMt)MliI-IX6hihi6.津hiM'(hiMihi 66yi-一),x xl hi式中 Mi= S (Xi).因此,只要確定了 Mi的值,就確定了整個(gè)表達(dá)式,Mi的計(jì)算方法如下:dihihi 4 1hi 1hi - hi 16(hi - hi 1yi 1 一小 y

5、i -Yiihi 1hi)= 6fXij,Xi,Xii那么Mi滿足如下n-1個(gè)方程:iM i2M i iM i 1 = di,i =1,2,.n -1常用的邊界條件有如下幾類(lèi):(D給定區(qū)間兩端點(diǎn)的斜率 m0,mn,即S r(x0) = y0 = m0, S r( xn) = yn = mn(2)給止區(qū)間兩端點(diǎn)的一階導(dǎo)數(shù) M0,Mn,即S"(x0) = y0 = M0,S*(xn) = y； = M n(3)假設(shè)y=f(x)是以b-a為周期的周期函數(shù),那么要求三次樣條插值函數(shù)S (x)也為周期函數(shù),對(duì) S(X)加上周期條件 S(p)(Xo +0)=S(p)(Xn-0),p = 0,1,

6、2對(duì)于第一類(lèi)邊界條件有6 y1 - y0 2M0 +M1 =( hhi6M nJ 2Mn = (mn - hn- mo)yn - ynl)hn對(duì)于第二類(lèi)邊界條件有2M0 + 九 0M1 =d0AMna +2Mn =dnd其中d6 0h16,(fX°,Xj -m0) 2(1 - 0川 0hn(m0 - fXn,Xn)2(1 -Un)Mn那么解就可以為2.九0.因.2.九1TmMM-d d dn 1.2. n 4.n.2Mndn./n對(duì)于第三類(lèi)邊界條件,y° = yn' M 0 = M n, S(X0 +0) =S(Xn 0),由此推得2M0 +%Mi +,Mn -1=

7、d.,其中h-, % =-h,d0 =-6一(fX0,Xi - f Xn -1, Xn ),那么解就可以為: hi - hnhi hnhi - hn.,*0八加.出.2.兒. Mi M 2dodid2.dn/:dn程序代碼:i拉格朗日插信函數(shù)Lang.mfunction f=lang(X,Y,xi)%X 為數(shù)據(jù)的橫坐標(biāo)%Y為數(shù)據(jù)的縱坐標(biāo)%xi插值點(diǎn)處的橫坐標(biāo)%f求得的拉格朗日插值多項(xiàng)式的值 n=length(X);f=0;for i=i:nl=i;for j=i:i-il=l.*(xi-X(j)/(X(i)-X(j);end ;for j=i+i:nl=l.*(xi-X(j)/(X(i)-X(

8、j);end ; %拉格朗日基函數(shù) f=f+l*Y(i);endfprintf( '%dn',f)return2牛頓插信函數(shù)newton.mfunction f=newton(X,Y,xi)%X 為數(shù)據(jù)的橫坐標(biāo)%Y為數(shù)據(jù)的縱坐標(biāo)%xi插值點(diǎn)處的橫坐標(biāo)%f求得的拉格朗日插值多項(xiàng)式的值 n=length(X);newt=X',Y'%計(jì)算差商表for j=2:nfor i=n:-1:1 if i>=jY(i)=(Y(i)-Y(i-1)/(X(i)-X(i-j+1);else Y(i)=0; endend newt=newt,Y'end% 計(jì)算牛頓插值f=

9、newt(1,2);for i=2:nz=1;for k=1:i-1 z=(xi-X(k)*z;endf=f+newt(i-1,i)*z;endfprintf( '%dn',f)return3三次樣條插值第一類(lèi)邊界條件Threch.mfunction S=Threch1(X,Y,dy0,dyn,xi)%X為數(shù)據(jù)的橫坐標(biāo)%Y為數(shù)據(jù)的縱坐標(biāo)%xi插值點(diǎn)處的橫坐標(biāo)%領(lǐng) 得的三次樣條插值函數(shù)的值%dy0左端點(diǎn)處的一階導(dǎo)數(shù)% dyn右端點(diǎn)處的一階導(dǎo)數(shù) n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(

10、1,n-2);for i=1:n%求函數(shù)的一階差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i)/h(i);endfor i=2:n%求函數(shù)的二階差商f2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1)/h(n-1);%?賦初值A(chǔ)=zeros(n+1,n+1); B=zeros(1,n-1); C=zeros(1,n-1); for i=1:n-1 B(i)=h(i)/(h(i)+h(i+1); C(i)=1-B(i); end

11、A(1,2)=1; A(n+1,n)=1; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=Ad; syms x; for i=1:n Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x-X(i) +M(i)/2*(x-X(i)A2+(M(i+1)-M(i)/(6*h(i)*(x-X(i)A3); digits(4); Sx(i)=vpa(Sx(i);%三樣條插值函數(shù)表達(dá)式end for i=1:n disp( 'S(x)='

12、;); fprintf( '%s (%d,%d)n',char(Sx(i),X(i),X(i+1);end for i=1:n if xi>=X(i)&&xi<=X(i+1)S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(xi-X(i)+M(i)/2*(xi-X(i)A2+(M (i+1)-M(i)/(6*h(i)*(xi-X(i)A3; end end disp( 'xi S' ); fprintf( '%d,%dn',xi,S);return4三次樣條插值第二類(lèi)邊界條件Threch2.m

13、functionSx=Threch2(X,Y,d2y0,d2yn,xi)X為數(shù)據(jù)的橫坐標(biāo)%Y為數(shù)據(jù)的縱坐標(biāo)%xi插值點(diǎn)處的橫坐標(biāo)%S得的三次樣條插值函數(shù)的值%d2y0左端點(diǎn)處的二階導(dǎo)數(shù)% d2yn右端點(diǎn)處的二階導(dǎo)數(shù)n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求一階差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i)/h(i);endfor i=2:n%求二階差商f2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*

14、f2(i);endd(1)=2*d2y0;d(n+1)=2*d2yn;%® 初值A(chǔ)=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1);C(i)=1-B(i);endA(1,2)=0;A(n+1,n)=0;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=Ad;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x

15、-X(i) +M(i)/2*(x-X(i)A2+(M(i+1)-M(i)/(6*h(i)*(x-X(i)A3);digits(4);Sx(i)=vpa(Sx(i);endfor i=1:ndisp( 'S(x)=');fprintf( '%s (%d,%d)n',char(Sx(i),X(i),X(i+1);endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(xi-X(i)+M(i)/2*(xi-X(i)2 +(M(i+1)-M(i)/

16、(6*h(i)*(xi-X(i)3;endenddisp( 'xi S' );fprintf( '%d,%dn',xi,S);return5插值節(jié)點(diǎn)處的插值結(jié)果main3.mclearcicX=0,0,0,1,0,2,0,3,0,4;Y=0.5000,0.5398,0.5793,0.6179,0.7554;xi=0.13;%xi=0.36;disp( 'xi=0.13');%disp('xi=0.36');disp(拉格朗日插值結(jié)果)； lang(X,Y,xi);disp('牛頓插彳I結(jié)果)；newton(X,Y,xi);

17、disp('三次樣條第一類(lèi)邊界條件插值結(jié)果 )；Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分別為兩端點(diǎn)處的一階導(dǎo)數(shù)disp('三次樣條第二類(lèi)邊界條件插值結(jié)果 )；Threch2(X,Y,0,-0.136,xi);%0,-0.136分別為兩端點(diǎn)處的二階導(dǎo)數(shù)6將多種插值函數(shù)即原函數(shù)圖像畫(huà)在同一張圖上main2.mclearcicX=0.0,0.1,0.2,0.3,0.4;Y=0.5000,0.5398,0.5793,0.6179,0.7554;a=iinspace(0,0.4,21);NUM=21;L=zeros(1,NUM);N=zeros(1,NUM

18、);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUMxi=a(i);L(i)=lang(X,Y,xi);%拉格朗日插值N(i)=newton(X,Y,xi);% 牛頓插值B(i)=normcdf(xi,0,1);%原函數(shù)S(i)=Threch1(X,Y,0.4,0.36,xi);%三次樣條函數(shù)第一類(lèi)邊界條件endplot(a,B,'-r');holdon;plot(a,L,'b');holdon;plot(a,N,'r');holdon;plot(a,S, 'r+'); hold on;legen

19、d'原函數(shù),拉格朗日插值,牛頓插值,三次樣條插值',2;hold off7增加插值節(jié)點(diǎn)觀察誤差變化main4.mclear;clc;N=5;%4.5第一問(wèn)Ini=zeros(1,1001);a=linspace(-1,1,1001);Ini=1./(1+25*a.A2);for i=1:3%節(jié)點(diǎn)數(shù)量變化次數(shù)N=2*N;t=linspace(-1,1,N+1);%插值節(jié)點(diǎn)ft=1./(1+25*t.A2);%f 值節(jié)點(diǎn)函數(shù)值val=linspace(-1,1,101);for j=1:101L(j)=lang(t,ft,val(j);S(j)=Threch1(t,ft,0.074

20、,-0.074,val(j);end%三樣條第一類(lèi)邊界條件插值plot(a,Ini, hold on'k' %原函數(shù)圖象plot(val,L,'r')%立格朗日插值函數(shù)圖像hold onplot(val,S,'b' ) %E次樣條插值函數(shù)圖像str=sprintf( '插值節(jié)點(diǎn)為d時(shí)的插值效果',N);title(str);legend('原函數(shù),拉格朗日插值,三次樣條插值)；%顯示圖例 hold offfigureend8車(chē)門(mén)曲線main5.mclearclcX=0,1,2,3,4,5,6,7,8,9,10;Y=0.0,

21、0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29;dy0=0.8;dyn=0.2;n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:nh(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i)/h(i);endfor i=2:nf2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1

22、)/h(n-1); A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1);C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=Ad;x=zeros(1,n);S=zeros(1,n);for i=1:nx(i)=X(i)+0.5;S(i)=丫(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x(i)-X(i)+M(

23、i)/2*(x(i)-X(i)A2+(M(i+1)-M(i)/(6*h(i)*(x(i)-X(i)3;endplot(X,Y,'k' ); hold on ;plot(x,S,'o');title('三次樣條插值效果圖')；legend('插值節(jié)點(diǎn),三次樣條插值)；hold off實(shí)驗(yàn)結(jié)果：4.31計(jì)算插值節(jié)點(diǎn)處的函數(shù)值xi=0.13 時(shí)1 3拉格蜩日炳值一黑 5.301坤幀插值精累三樣綏那一類(lèi)邊界奉件插值經(jīng)果 S?注) =.5000+.400a*K-.1, OOOOOOe-OO1>S=.5 1 4T- 1 O, 9 曰*W 口了

24、曰年一< 1 , OOOOOOe-OO i ODOOOOe-OO 1 5S=.1 093-1-39. 68*k 326. 392-1-6. 口4口*冗 <2. OOOOOOe-OO Ij. 3. OCQDOOe-00 1 ) S <xi =2. 4 7 6-4?. 93*h 3*豈之,4.6*k:"2-17. til*K <3. OOUDOOe-OO1, Q, ODOOOUe-UO 1 Jxi S1 . 30aOOOe-ODlF 5. S31403e-0Ol 三樣奈第二類(lèi)邊界舉件插值結(jié)果 S <ic> =.5Q00-I-. 3805*k4-1. 7

25、64*3： - 3(0, 1. OOOOOOs-OOl JS <京) =.5 10 8-9, O7O + k*3+3. 2472-t- 5S73e-l*x< 1, OOOOOOe-OO U 2, OOOOOOe-OO 1 >2 =.I 6fi3-»-r33, 93*x " 3-22. 55率¥2+5- 215 增 x<2-門(mén)口口口口口已一0口1, 3. 口口口門(mén) 口口門(mén)一口口1)S=工.SO7 2S.&Q*ie * 3-1-32. 14" N L 1 . 1 0*督(3, 口 口口口欄一0口 1, 4- OOOOOOeOD 1 >S1. SOOOOOe-OOl, O5. 5302 1 1 e-OOlj. >>Xi=0.36 時(shí)拉格朗日插值絹黑#頓捕值結(jié)果 ?. 1 693S3e-OO1 三才羊傍弟一類(lèi)邊界多仔括值絡(luò)果.50-«-* dODO*K. 3E9B*3E " Z-I-3. 09S