數(shù)值分析試驗報告1

1、百度文庫-讓每個人平等地提升自我實驗報告實驗項目名稱插值法實驗室數(shù)學實驗室所屬課程名稱數(shù)值逼近實驗類型算法設計實驗日期班級學號姓名成績11百度文庫-讓每個人平等地提升自我實驗概述：【實驗目的及要求】本次實驗的目的是熟練數(shù)值分析第二章“插值法”的相關內容，掌握三種插值方法：牛頓多本次試驗要求編寫牛頓多項式插值，三次樣條插值，拉格朗日插值的程序編碼，并在MATLABm【實驗原理】數(shù)值分析第二章“插值法”的相關內容，包括：牛頓多項式插值，三次樣條插值，拉格朗日插值【實驗環(huán)境】(使用的軟硬件)軟件：MATLAB2012a硬件：電腦型號：聯(lián)想Lenovo昭陽E46A筆記本電腦、操作系統(tǒng)：Windows8

2、專業(yè)版處理器：Intel(R)Core(TMi3CPUM350實驗內容：【實驗方案設計】第一步，將書上關于三種插值方法的內容轉化成程序語言，用MATLA取現(xiàn)；第二步，分別用【實驗過程】(實驗步驟、記錄、數(shù)據(jù)、分析)實驗的主要步驟是：首先分析問題，根據(jù)分析設計MATLABi序，利用程序算出問題答案，分析所實驗一：已知函數(shù)在下列各點的值為Xi/.、f(Xi)試用4次牛頓插值多項式P4(x)及三次樣條函數(shù)S(x)(自然邊界條件)對數(shù)據(jù)進行插值。用圖給出(Xi,(1)首先我們先求牛頓插值多項式，此處要用4次牛頓插值多項式處理數(shù)據(jù)。已知n次牛頓插值多項式如下：Pn=f(x0)+fX0,x1(x-x0)+

3、fx0,x1,X2(x-x0)(x-x1)+zfx0,x1,Xn(x-X0)(x-xn-1)我們要知道牛頓插值多項式的系數(shù)，即均差表中得部分均差。在MATLAB勺Editor中輸入程序代碼，計算牛頓插值中多項式系數(shù)的程序如下：functionvarargout=newtonliu(varargin)/clear,clcx=；22百度文庫-讓每個人平等地提升自我fx=;newtonchzh(x,fx);functionnewtonchzh(x,fx)他此函數(shù)可得差分表n=length(x);fprintf(*差分表*n);FF=ones(n,n);FF(:,1)=僅;fori=2:nforj=i

4、:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-i+1);endendfori=1:nfprintf(%4.2f,x(i);forj=1:ifprintf(%10.5f,FF(i,j);endfprintf(n);end由MATLA即算得:xif(xi)一階差商二階差商三階差商四階差商所以有四次插值牛頓多項式為:P4(x)=)(2)接下來我們求三次樣條插值函數(shù)。用三次樣條插值函數(shù)由上題分析知，要求各點的M值:33百度文庫-讓每個人平等地提升自我0.5000.500MoMi-3.75000.500020.5000.50020.500-4.5000-6.750

5、0三次樣條插值函數(shù)計算的程序如下:functiontgsanci(n,s,t)%nM4代表元素數(shù)，s,t代表端點的一階導。x=；y=;n=5forj=1:1:n-1h(j)=x(j+1)-x(j);endforj=2:1:n-1r(j)=h(j)/(h(j)+h(j-1);endforj=1:1:n-1u(j)=1-r(j);endforj=1:1:n-1f(j)=(y(j+1)-y(j)/h(j);endforj=2:1:n-1d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);endd(1)=0d(n)=0a=zeros(n,n);forj=1:1:na(j,j)=2;end

6、r(1)=0;u(n)=0;forj=1:1:n-1a(j+1,j)=u(j+1);44百度文庫-讓每個人平等地提升自我endzb=inv(a)/m=b*d/p=zeros(n-1,4);%p矩陣為S(x)函數(shù)的系數(shù)矩陣forj=1:1:n-1p(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)A2/6)/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)、2/6)/h(j);end得到m=(00)即M0=O;M1=;M2=;M3=;M4=0則根據(jù)三次樣條函數(shù)定義，可得：0(0.4x)31.3393x0

7、.2)4.9000.4x)4.6536(x0.2),x0.2,0,4S(x)=-1.33930.6x)30.8929(x-0.4)34.65360.6x)4.0857(x0.4),x0.4,0.6-0.89290.8x)3-2.5893x-0.6)34.0857(0.8x)3.3036x0.6),x0.6,0.81-2.58931.0x)3-0(x0.8)33.3036(1.0x)1.9(x0.8),x0.8.1.0接著，在CommandWindow!輸入畫圖的程序代碼,下面是畫牛頓插值以及三次樣條插值圖形的程序：x=;y=;plot(x,y)holdonfori=1:1:5、y(i)=*(x

8、(i)*(x(i)*(x(i)*(x(i)*(x(i)*(x(i)*(x(i)endk=011011x0=+*kfori=1:1:4y0(i)=*(x(i)*(x(i)*(x(i)*(x(i)*(x(i)*(x(i)*(x(i)endplot(x0,y0,o,x0,y0)55百度文庫-讓每個人平等地提升自我holdony1=spline(x,y,x0)人plot(x0,y1,o)/holdon/s=csape(x,y,variational)/fnplt(s,r)/holdon/gtext(三次樣條自然邊界,)gtext(原圖像)、gtext(,4次牛頓插值)運行上述程序可知：給出的(*,yi

9、),xi=+,i=0,1,11,10點，S(x)及P4(x)圖形如下所示:實驗二：1，在區(qū)間-1,1上分別取n10,20用兩組等距節(jié)點對龍格函數(shù)f(x)2作多項式插值及三次樣條插值,125x2我們先求多項式插值：在MATLAB勺Editor中建立一個多項式的Mfile,輸入如下的命令(如牛頓插值公式)：functionC,D=newpoly(X,Y)n=length(X);/D=zeros(n,n)/D(:,1)=Yforj=2:n66百度文庫-讓每個人平等地提升自我fork=j:nD(k,j)=(D(k,j-1)-D(k-1,j-1)/(X(k)-X(k-j+1);end/end/C=D(n

10、,n);fork=(n-1):-1:1/C=conv(C,poly(X(k)m=length(C);C(m)=C(m)+D(k,k);end當n=10時，我們在CommandWindow43輸入以下的命令:clear,clf,holdon;X=-1:1;Y=1./(1+25*X.A2);C,D=newpoly(X,Y);x=-1:1;y=polyval(C,x);plot(x,y,X,Y,.);gridon;xp=-1:1;z=1./(1+25*xp.A2);plot(xp,z,r)得到插值函數(shù)和f(x)圖形：77百度文庫-讓每個人平等地提升自我當n=20時，我們在CommandWindow4

11、3輸入以下的命令:clear,clf,holdon;X=-1:1;Y=1./(1+25*X.A2);C,D=newpoly(X,Y);x=-1:1;y=polyval(C,x);plot(x,y,X,Y,.);gridon;xp=-1:1;z=1./(1+25*xp.A2);plot(xp,z,r)得到插值函數(shù)和f(x)圖形：88百度文庫-讓每個人平等地提升自我i*1/1ki.ii|!i_,*LJU-50J.IiI!1i1d1*.*1i1電I11!111*覃1I|1-1JJH!i-111.*D-fcj1司84)6-0.4聞20020.406081下面再求三次樣條插值函數(shù)，在MATLAB勺Edi

12、tor中建立一個多項式的M-file,輸入下列程序代碼：functionS=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);fork=2:N-1temp=A(k-1)/B(k-1);、/B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);

13、/end/M(N)=U(N-1)/B(N-1);/fork=N-2:-1:199百度文庫-讓每個人平等地提升自我M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endzM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;/M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;/fork=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*M(k+1)+M(k+2)/6;S(k+1,4)=Y(k+1);end當n=10時，我們在CommandWindow

14、43輸入以下的命令:clear,clcX=-1:1;Y=1./(25*X.A2+1);dx0=;dxn=;S=csfit(X,Y,dx0,dxn)x1=-1:;y1=polyval(S(1,:),x1-X)；x2=:0;y2=polyval(S(2,:),x2-X(2);x3=0:;y3=polyval(S(3,:),x3-X(3);x4=:1;y4=polyval(S(4,:),x4-X(4);Plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,.)結果如圖:當n=20時，我們在CommandWindow43輸入以下的命令:1010百度文庫-讓每個人平等地提升自我clear,c

15、lcX=-1:1;Y=1./(25*X.A2+1);dx0=;dxn=;S=csfit(X,Y,dx0,dxn)/x1=-1:;y1=polyval(S(1,:),x1-X)；x2=:0;y2=polyval(S(2,:),x2-X(2);x3=0二；y3=polyval(S(3,:),x3-X(3);x4=:1;y4=polyval(S(4,:),x4-X(4);Plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,.)結果如圖：實驗三：x0149-162536y0123456下列數(shù)據(jù)點的插值47可以得到平方根函數(shù)的近似，在區(qū)間0,64上作圖。/(1)用這9各點作8次多項式插值L

16、8(x).(2)用三次樣條(自然邊界條件)程序求S(x)。從結果看在0,64上，那個插值更精確；在區(qū)間0,1nL8(x)可由公式Ln(x)=yklk(xj)得出。/k0三次樣條可以利用自然邊界條件。寫成矩陣：1111百度文庫-讓每個人平等地提升自我20000Mo12100Mi02220M200323M300042M4hghj其中j=_j1，,，=j,dj=6fxj-1,xj,xj+1,n=0=0d0=dn=0hkhjhj-1hj-xl0(x)=：(x-1)(x4)(x9)(x16)(x25)(x36)(x49)(x64)(01)(04)(09)(016)(025)(036)(049)(064)

17、11(x)=：(x-0)(x4)(x9)(x16)(x25)(x36)(x49)(x64)(10)(14)(19)(116)(125)(136)(149)(164)12(x)=:(x-0)(x1)(x9)(:x16)(x25)(x36)(x49)(x64)(40)(41)(49)(416)(425)(436)(449)(464)13(x)=：(x-0)(x1)(x4)(x16)(x25)(x36)(x49)(x64)(90)(91)(94)(916)(925)(936)(949)(964)14(x)=(x-0)(x1)(x4)(x19)(x25)(x36)(x49)(x64)(160)(161

18、)(164)(169)(1625)(1636)(1649)(1664)15(x)=(x-0)(x-1)(x4)(x9)(x16)(x36)(x49)(x64)(250)(251)(254)(259)(2516)(2536)(2549)(2564)16(x)=(x-0)(x1)(x4)(x9)(x16)(x25)(x49)(x64)(360)(361)(364)(369)(3616)(3625)(3649)(3664)17(x)=(x-0)(x-1)(x4)(x9)(x16)(x25)(x36)(x64)(490)(491)(494)(499)(49，16)(4925)(4936)(4964)1

19、8(x)=(x-0)(x-1)(x4)(x9)(x16)(x25)(x36)(x49)(640)(641)(644)(649)(6416)(6425)(6436)(6449)L8(x)=：1i(x)+212(x)+313(x)+414(x)+515(x)+616(x)+717(x)+818(x)求二次樣條插值函數(shù)由MAILA時算：可得矩陣形式的線性方程組為：2.000000000000Mo00.25002.00000.7500000000Ml1.000.37502.00000.625000000M20.1000.41672.00000.58330000M30.02860000.43752.00

20、000.5625000M40.011900000.45002.00000.550000M50.0061000000.45832.00000.54170M60.00350000000.46342.00000.5357M70.0022000000002.0000M80在MATLABHEditor中輸入程序代碼，以下是三次樣條函數(shù)的程序代碼functiontgsanci(n,s,t)%n代表兀素數(shù),s,t代表端點的一階導。1212百度文庫-讓每個人平等地提升自我y=012345678;x=01491625364964;n=9/forj=1:1:n-1/h(j)=x(j+1)-x(j);/end/fo

21、rj=2:1:n-1r(j)=h(j)/(h(j)+h(j-1);endforj=1:1:n-1u(j)=1-r(j);endforj=1:1:n-1f(j)=(y(j+1)-y(j)/h(j);endforj=2:1:n-1d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);endd(1)=0d(n)=0a=zeros(n,n);forj=1:1:na(j,j)=2;endr(1)=0;u(n)=0;forj=1:1:n-1a(j+1,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ù)的系

22、數(shù)矩陣1313百度文庫-讓每個人平等地提升自我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)A2/6)/h(j);P(j，4)=(y(j+1)-m(j+1)*(h(j)A2/6)/h(j);endP解得：M=0；Mi=；M2=；M3=；M4=；M5=；M6=；M7=；M8=0,則三次樣條函數(shù)：0(1x)30.0868(x0)30(1x)1.0868(x-0),x0,1-0.0289(4x)3-0.0031(x-1)30.5938(4x)0.6388(x1)0.0019(9x)30.0009(x-4)30.353

23、5(9x)0.6271(x4)x1,4x4,9S(x)=-0.0006(16x)30(x9)30.4590(160(25x)3-0.0001(x16)3-0.4436(250(36x)3-0(x25)30.4599(36x)0(48x)3-0(x36)30.4633(48x)0(64x)30(x48)30.4689(64x)x)-0.5708(x-9),x9,16x)0.5600(x-16),x16,250.5470(x-25),x25,360.5404(x-36),x36,480.5333(x-48),x48,64卜面進行畫圖，在CommandWindo則輸入畫圖的程序代碼:則圖形比較那個插值更精確的函數(shù):x0=01491625364964;y0=012345678;x=0:64;y=lagr1(x0,y0,x);p