數(shù)值分析實驗報告(共17頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上實驗報告實驗項目名稱 插值法 實驗室數(shù)學實驗室 所屬課程名稱 數(shù)值逼近 實 驗 類 型 算法設(shè)計 實 驗 日 期 班 級 學 號 姓 名 成 績 實驗概述：【實驗?zāi)康募耙蟆勘敬螌嶒灥哪康氖鞘炀殧?shù)值分析第二章“插值法”的相關(guān)內(nèi)容，掌握三種插值方法：牛頓多項式插值，三次樣條插值，拉格朗日插值，并比較三種插值方法的優(yōu)劣。本次試驗要求編寫牛頓多項式插值，三次樣條插值，拉格朗日插值的程序編碼，并在MATLAB軟件中去實現(xiàn)?！緦嶒炘怼繑?shù)值分析第二章“插值法”的相關(guān)內(nèi)容，包括：牛頓多項式插值，三次樣條插值，拉格朗日插值的相應(yīng)算法和相關(guān)性質(zhì)?！緦嶒灜h(huán)境】（使用的軟硬件）軟件：MA

2、TLAB 2012a硬件：電腦型號：聯(lián)想 Lenovo 昭陽E46A筆記本電腦操作系統(tǒng)：Windows 8 專業(yè)版 處理器：Intel（R）Core（TM）i3 CPU M 350 2.27GHz 2.27GHz實驗內(nèi)容：【實驗方案設(shè)計】第一步，將書上關(guān)于三種插值方法的內(nèi)容轉(zhuǎn)化成程序語言，用MATLAB實現(xiàn)；第二步，分別用牛頓多項式插值，三次樣條插值，拉格朗日插值求解不同的問題?！緦嶒炦^程】（實驗步驟、記錄、數(shù)據(jù)、分析）實驗的主要步驟是：首先分析問題，根據(jù)分析設(shè)計MATLAB程序，利用程序算出問題答案，分析所得答案結(jié)果，再得出最后結(jié)論。實驗一：已知函數(shù)在下列各點的值為xi0.20.40.6.0

3、.81.0f（xi）0.980.920.810.640.38試用4次牛頓插值多項式P4（x）及三次樣條函數(shù)S（x）（自然邊界條件）對數(shù)據(jù)進行插值。用圖給出（xi，yi），xi=0.2+0.08i，i=0，1, 11, 10，P4（x）及S（x）。(1)首先我們先求牛頓插值多項式，此處要用4次牛頓插值多項式處理數(shù)據(jù)。已知n次牛頓插值多項式如下：Pn=f(x0)+fx0,x1(x-x0)+ fx0,x1,x2(x-x0) (x-x1)+···+ fx0,x1，···xn(x-x0) ···(x-xn-1)我們

4、要知道牛頓插值多項式的系數(shù)，即均差表中得部分均差。在MATLAB的Editor中輸入程序代碼，計算牛頓插值中多項式系數(shù)的程序如下：function varargout=newtonliu(varargin)clear,clcx=0.2 0.4 0.6 0.8 1.0;fx=0.98 0.92 0.81 0.64 0.38;newtonchzh(x,fx);function newtonchzh(x,fx)%由此函數(shù)可得差分表n=length(x);fprintf('*差分表*n');FF=ones(n,n);FF(:,1)=fx'for i=2:n for j=i:n

5、FF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-i+1); endendfor i=1:n fprintf('%4.2f',x(i); for j=1:i fprintf('%10.5f',FF(i,j); end fprintf('n');end由MATLAB計算得：xi f(xi) 一階差商二階差商三階差商四階差商0.200.9800.400.920-0.300000.600.810-0.55000-0.625000.800.640-0.85000-0.75000-0.208331.000.380-1.300

6、00-1.12500-0.62500-0.52083所以有四次插值牛頓多項式為：P4（x）=0.98-0.3(x-0.2)-0.62500 (x-0.2)(x-0.4) -0.20833 (x-0.2)(x-0.4)(x-0.6)-0.52083 (x-0.2)(x-0.4)(x-0.6)(x-0.8）(2)接下來我們求三次樣條插值函數(shù)。用三次樣條插值函數(shù)由上題分析知,要求各點的M值：三次樣條插值函數(shù)計算的程序如下：function tgsanci(n,s,t) %n代表元素數(shù)，s,t代表端點的一階導。x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.

7、38; n=5for j=1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1 u(j)=1-r(j);endfor j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);end d(1)=0 d(n)=0 a=zeros(n,n);for j=1:1:n a(j,j)=2;end r(1)=0; u(n)=0;for j=1:1:n-1 a(j+1,j)=u(j+1); a

8、(j,j+1)=r(j);endb=inv(a)m=b*d'p=zeros(n-1,4); %p矩陣為S(x)函數(shù)的系數(shù)矩陣for j=1:1:n-1 p(j,1)=m(j)/(6*h(j); p(j,2)=m(j+1)/(6*h(j); p(j,3)=(y(j)-m(j)*(h(j)2/6)/h(j); p(j,4)=(y(j+1)-m(j+1)*(h(j)2/6)/h(j);end p得到m=(0 -1.6071 -1.0714 -3.1071 0)T 即M0=0 ;M1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0則根據(jù)三次樣條函數(shù)定義，可得：S

9、(x)= 接著，在Command Window里輸入畫圖的程序代碼，下面是畫牛頓插值以及三次樣條插值圖形的程序：x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38;plot(x,y)hold on for i=1:1:5y(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4) -0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endk=0 1 10 11x0=0.

10、2+0.08*kfor i=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4) -0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endplot( x0,y0,'o',x0,y0 )hold ony1=spline(x,y,x0)plot(x0,y1,'o')hold ons=csape(x,y,'variational') fnplt(s,&

11、#39;r')hold ongtext('三次樣條自然邊界')gtext('原圖像')gtext('4次牛頓插值')運行上述程序可知：給出的（xi，yi），xi=0.2+0.08i，i=0，1, 11, 10點，S（x）及P4（x）圖形如下所示：實驗二：在區(qū)間-1,1上分別取用兩組等距節(jié)點對龍格函數(shù)作多項式插值及三次樣條插值，對每個值，分別畫出插值函數(shù)即的圖形。我們先求多項式插值：在MATLAB的Editor中建立一個多項式的M-file,輸入如下的命令（如牛頓插值公式）：function C,D=newpoly(X,Y)n=length

12、(X);D=zeros(n,n) D(:,1)=Y' for j=2:n for k=j:n D(k,j)=(D(k,j-1)- D(k-1,j-1)/(X(k)-X(k-j+1); endendC=D(n,n);for k=(n-1):-1:1 C=conv(C,poly(X(k) m=length(C); C(m)= C(m)+D(k,k);end當n=10時，我們在Command Window中輸入以下的命令：clear,clf,hold on; X=-1:0.2:1; Y=1./(1+25*X.2); C,D=newpoly(X,Y); x=-1:0.01:1; y=polyv

13、al(C,x); plot(x,y,X,Y,'.'); grid on; xp=-1:0.2:1; z=1./(1+25*xp.2); plot(xp,z,'r')得到插值函數(shù)和f（x）圖形：當n=20時，我們在Command Window中輸入以下的命令：clear,clf,hold on; X=-1:0.1:1; Y=1./(1+25*X.2); C,D=newpoly(X,Y); x=-1:0.01:1; y=polyval(C,x); plot(x,y,X,Y,'.'); grid on; xp=-1:0.1:1; z=1./(1+25*

14、xp.2); plot(xp,z,'r')得到插值函數(shù)和f（x）圖形：下面再求三次樣條插值函數(shù)，在MATLAB的Editor中建立一個多項式的M-file,輸入下列程序代碼：function S=csfit(X,Y,dx0,dxn) N=length(X)-1;H=diff(X); D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N); C=H(2:N); U=6*diff(D);B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(-D(N); fo

15、r k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end M(N)=U(N-1)/B(N-1);for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1); S(k+1,2)=M(k+1)/2; S(k+1,3)=D(k+1)-H(k+1)*(2*

16、M(k+1)+M(k+2)/6; S(k+1,4)=Y(k+1);end當n=10時，我們在Command Window中輸入以下的命令：clear,clcX=-1:0.2:1;Y=1./(25*X.2+1); dx0= 0.14201;dxn= -0.14201;S=csfit(X,Y,dx0,dxn) x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1);x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2); x3=0:0.01:0.5; y3=polyval(S(3,:),x3-X(3);x4=0.5:0.01:1;y4=poly

17、val(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,'.')結(jié)果如圖：當n=20時，我們在Command Window中輸入以下的命令：clear,clcX=-1:0.1:1;Y=1./(25*X.2+1); dx0= 0.14201;dxn= -0.14201;S=csfit(X,Y,dx0,dxn) x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1);x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2); x3=0:0.01:0.5; y3=polyval(S

18、(3,:),x3-X(3);x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,'.')結(jié)果如圖：實驗三：下列數(shù)據(jù)點的插值x01491625364964y012345678可以得到平方根函數(shù)的近似，在區(qū)間0,64上作圖。（1）用這9各點作8次多項式插值L8(x).（2）用三次樣條（自然邊界條件）程序求S（x）。從結(jié)果看在0,64上，那個插值更精確；在區(qū)間0,1上，兩種哪個更精確？L8(x)可由公式Ln(x)=得出。 三次樣條可以利用自然邊界條件。寫成矩陣：其中j=，i=，dj=6f

19、xj-1,xj,xj+1，n=0=0 d0=dn=0l0(x)=l1(x)= l2(x)= l3(x)= l4(x)= l5(x)= l6(x)= l7(x)= l8(x)= L8(x)= l1(x)+2 l2(x)+3 l3(x)+4 l4(x)+5 l5(x)+6 l6(x)+7 l7(x)+8 l8(x)求三次樣條插值函數(shù)由MATLAB計算：可得矩陣形式的線性方程組為：在MATLAB中的Editor中輸入程序代碼，以下是三次樣條函數(shù)的程序代碼：function tgsanci(n,s,t) %n代表元素數(shù)，s,t代表端點的一階導。y=0 1 2 3 4 5 6 7 8;x=0 1 4 9

20、 16 25 36 49 64; n=9for j=1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1 u(j)=1-r(j);endfor j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);end d(1)=0 d(n)=0 a=zeros(n,n);for j=1:1:n a(j,j)=2;end r(1)=0; u(n)=0;for j=1:1:n-1 a(j+1

21、,j)=u(j+1); a(j,j+1)=r(j);endb=inv(a)m=b*d't=ap=zeros(n-1,4); %p矩陣為S(x)函數(shù)的系數(shù)矩陣for j=1:1:n-1 p(j,1)=m(j)/(6*h(j); p(j,2)=m(j+1)/(6*h(j); p(j,3)=(y(j)-m(j)*(h(j)2/6)/h(j); p(j,4)=(y(j+1)-m(j+1)*(h(j)2/6)/h(j);end p解得:M0=0;M1=-0.5209;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=-0.0009;M

22、8=0，則三次樣條函數(shù)：S(x)= 下面進行畫圖，在Command Window中輸入畫圖的程序代碼：%畫圖形比較那個插值更精確的函數(shù)：x0=0 1 4 9 16 25 36 49 64;y0=0 1 2 3 4 5 6 7 8;x=0:64;y=lagr1(x0,y0,x);plot(x0,y0,'o')hold onplot(x,y,'r');hold on;pp=csape(x0,y0,'variational')fnplt(pp,'g');hold on;plot(x0,y0,':b');hold on%axis(0 2 0 1); %看0 1區(qū)間的圖形時加上這條指令gtext('三次樣條插值')gtext('原圖像')gtext('拉格朗日插值')%下面是求拉格朗