數(shù)值分析實驗報告.docx

1、.實驗報告實驗項目名稱插值法實驗室數(shù)學實驗室所屬課程名稱數(shù)值逼近實驗類型算法設(shè)計實驗日期班級學號姓名成績優(yōu)質(zhì)范文.實驗概述：【實驗?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)境】（使用的軟硬件）軟件：MATLAB 2012a硬件：電腦型號：聯(lián)想Lenov

2、o 昭陽 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.81.0f （ xi ） 0.980.9

3、20.810.640.38試用 4 次牛頓插值多項式P4（ x）及三次樣條函數(shù)S（ x）（自然邊界條件）對數(shù)據(jù)進行插值。用圖給出 （ xi ，y）， x=0.2+0.08i， i=0， 1, 11, 10，P （ x）及 S（x）。ii4(1) 首先我們先求牛頓插值多項式，此處要用4 次牛頓插值多項式處理數(shù)據(jù)。已知 n 次牛頓插值多項式如下：P =f(x)+fx,x(x-x)+ fx0,x ,x (x-x0) (x-x)+ + fx,x， x (x-x) (x-xn-1)n001012101n0優(yōu)質(zhì)范文.我們要知道牛頓插值多項式的系數(shù)，即均差表中得部分均差。在 MATLAB的 Editor中輸

4、入程序代碼，計算牛頓插值中多項式系數(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);functionnewtonchzh(x,fx)%由此函數(shù)可得差分表n=length(x);fprintf(*差分表 *n);FF=ones(n,n);FF(:,1)=fx;for i=2:nfor j=i:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-i+1);endendfor i

5、=1:nfprintf(%4.2f,x(i);for j=1:ifprintf(%10.5f,FF(i,j);endfprintf(n);end由 MATLAB計算得：xif(x i )一階差商二階差商三階差商四階差商0.200.9800.400.920-0.300000.600.810-0.55000-0.625000.800.640-0.85000-0.75000-0.208331.000.380-1.30000-1.12500-0.62500-0.52083所以有四次插值牛頓多項式為：P4（ x） =0.98-0.3(x-0.2)-0.62500 (x-0.2)(x-0.4) -0.20

6、833 (x-0.2)(x-0.4)(x-0.6)-0.52083優(yōu)質(zhì)范文.(x-0.2)(x-0.4)(x-0.6)(x-0.8）(2) 接下來我們求三次樣條插值函數(shù)。用三次樣條插值函數(shù)由上題分析知, 要求各點的 M值：20000M 000.50020.50000M 1- 3.750000.50020.5000M 2- 4.5000000.50020.500M 3- 6.750000002M 40三次樣條插值函數(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

7、 0.64 0.38;n=5for j=1:1:n-1h(j)=x(j+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);endd(1)=0d(n)=0a=zeros(n,n);for j=1:1:na(j,j)=2;end優(yōu)質(zhì)范文.r(1)=0;u(n)=0;for j=1:1:n-1a(j+1,j)=u(j+1);a(

8、j,j+1)=r(j);endb=inv(a)m=b*dp=zeros(n-1,4); %p矩陣為 S(x) 函數(shù)的系數(shù)矩陣for j=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)2/6)/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)2/6)/h(j);endp得到 m=(0 -1.6071-1.0714 -3.10710)T即 M0=0 ;M 1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0則根據(jù)三次樣條函數(shù)定義，可得：（x3（）（），x0

9、.2,0.4）0 0.41.3393 x0.24.900 0.4x 4.6536( x 0.2)（3（3（）（），0.4,0.6）S(x)=-1.33930.6x0.8929 x - 0.44.65360.6x4.0857 x0.4x33x0.6,0.8（）（），）- 0.89290.8x- 2.5893 x - 0.64.0857 0.8x 3.3036 x0.6（x3（3（）1.9( x，x0.8,1.0）- 2.58931.0- 0 x0.83.3036 1.0 x0.8)接著，在Command Window里輸入畫圖的程序代碼，下面是畫牛頓插值以及三次樣條插值圖形的程序：x=0.2 0

10、.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38;plot(x,y)hold onfor 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 11優(yōu)質(zhì)范文.x0=0.2+0.08*kfor i=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0

11、.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,r)hold ongtext(三次樣條自然邊界 )gtext(原圖像 )gtext(4次牛頓插值 )運行上述程序可知：給出的 （ xi ， yi ）， xi =0.2+0.08i，i=0

12、 ， 1, 11, 10點， S（ x）及 P4（ x）圖形如下所示：10.94次牛頓插值0.80.7三次樣條自然邊界0.6原 圖 像0.50.40.30.20.30.40.50.60.70.80.911.11.20.2實驗二：在區(qū)間 -1,1上分別取 n10,20 用兩組等距節(jié)點對龍格函數(shù)f ( x)1作多項式插值及三次樣條插25x21優(yōu)質(zhì)范文.值，對每個 n 值，分別畫出插值函數(shù)即f ( x) 的圖形。我們先求多項式插值：在 MATLAB的 Editor中建立一個多項式的M-file,輸入如下的命令（如牛頓插值公式）：function C,D=newpoly(X,Y)n=length(X)

13、;D=zeros(n,n)D(:,1)=Yfor j=2:nfor k=j:nD(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:1C=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=polyval(C,x);plot(x,y,X

14、,Y,.);grid on;xp=-1:0.2:1;z=1./(1+25*xp.2);plot(xp,z,r)得到插值函數(shù)和f （x）圖形：優(yōu)質(zhì)范文.當 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*xp.2);plot(xp,z,r)得到插值函數(shù)和f （x）圖形：優(yōu)質(zhì)范文.下面再求三次樣條插值函數(shù)，在M

15、ATLAB的 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);for k=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-te

16、mp*U(k-1);endM(N)=U(N-1)/B(N-1);for k=N-2:-1:1優(yōu)質(zhì)范文.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-1S(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 時，我們在 Command Window中輸入以

17、下的命令：clear,clcX=-1:0.2:1;Y=1./(25*X.2+1);dx0= 0.0739644970414201;dxn= -0.0739644970414201;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=polyval(S(4,:),x4-X(4);plot(x1,y1,x2,y2,x3,y3,x4,

18、y4, X,Y,.)結(jié)果如圖：當 n=20 時，我們在 Command Window中輸入以下的命令：優(yōu)質(zhì)范文.clear,clcX=-1:0.1:1;Y=1./(25*X.2+1);dx0= 0.0739644970414201;dxn= -0.0739644970414201;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;y

19、4=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 上，兩種哪個更精確？ny k l k ( x j ) 得出。L (x) 可由公式 L (x)=8nk0三次樣條可以利用自然邊界條件。寫成矩陣：優(yōu)質(zhì)范文.20000M 0d0

20、12100M 1d102220M 2d 200323M 3d300042M 4d 4其中j =h j-1，i=h j， dj=6fxj-1,x j ,x j+1 ，n=0=0 d 0=dn=0h j-1h jh j-1h jl 0(x)=(x - 1)( x4 )( x9)( x16 )( x25)( x 36 )( x49 )( x64 )(0 1)( 04)( 09)(016)( 0 25 )(036)(049)(064 )l(x)=(x - 0)( x4)( x9)( x16)( x25)( x36)( x49 )( x64 )1(10)(14 )(19)(116)(125)(136)(

21、149)(164 )l2(x)=(x - 0)( x1)( x9)( x16)( x25)( x36)( x49)( x64)( 40)( 41)( 49)( 416)( 425)( 436)( 449)( 464)l(x)=(x - 0)( x1)( x4)( x16)( x25)( x36)( x49)( x64)3(90)( 91)( 94)( 916)( 925)( 936)( 949)(964)l 4(x)=(x - 0)( x1)( x4)( x9)( x25)( x36 )( x49)( x64)(160)(161)(164)(169)(1625)(1636)(1649)(166

22、4)l(x)=( x - 0) (x - 1)( x4)( x9)( x16)( x36)( x49)( x64)5( 250)( 251)( 254)( 259)( 2516 )( 2536)( 2549)( 2564)l(x)=(x - 0)( x1)( x4)( x9)( x16)( x25)( x49)( x64)6(360)(361)( 364)( 369)( 3616 )(3625)( 3649)( 3664 )（）x4)(x9)(x16)(x25)(x36)(x64)l 7(x)=x - 0(x - 1)(（49）1)( 494)( 499)( 4916)( 49（36)( 49

23、64)0 (4925) 49l(x)=( x - 0)(x- 1)( x4)( x9 )( x16 )( x25 )( x36 )( x49 )8( 640)( 641)( 644)( 649)( 6416 )( 6425 )( 6436 )( 6449 )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計算：可得矩陣形式的線性方程組為：2.000000000000M 00.25002.00000.7500000000M 100.37502.00000.625000000M

24、2000.41672.00000.58330000M 30000.43752.0000 0.5625000M 400000.45002.00000.550000M 5000000.45832.00000.54170M 60000000.46342.00000.5357M 7000000002.0000M 8在 MATLAB中的 Editor 中輸入程序代碼，以下是三次樣條函數(shù)的程序代碼：function tgsanci(n,s,t) %n代表元素數(shù)， s,t 代表端點的一階導。01 .00.10 .02860 .01190 .00610 .00350 .00220優(yōu)質(zhì)范文.y=0 1 2 3

25、4 5 6 7 8;x=0 1 4 9 16 25 36 49 64;n=9for j=1:1:n-1h(j)=x(j+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);endd(1)=0d(n)=0a=zeros(n,n);for j=1:1:na(j,j)=2;endr(1)=0;u(n)=0;for j=1:1:n-

26、1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*dt=ap=zeros(n-1,4); %p矩陣為 S(x) 函數(shù)的系數(shù)矩陣優(yōu)質(zhì)范文.for j=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)2/6)/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)2/6)/h(j);endp解得 :M=0;M =-0.5209;M=0.0558;M=-0.0261;M=0.0006;M5=-0.0029;M=-0.0008;M=-0.0009;0123

27、467M=0，則三次樣條函數(shù)：8（x30.0868( x0)3（x）1.0868(x，0,1）0 10 1- 0) x（4x3（3（4x）（x），x1,4）- 0.0289- 0.0031x - 10.59380.63881（x3（3（9x）（x），x 4,9）0.0019 90.0009x - 40.35350.62714（x3（3（x），x9,16S(x)=）- 0.0006160 x90.4590 16- 0.5708(x- 9)0x3（163（25x）0.5600(x - 16)，x16,25（）0 25- 0.0001x- 0.4436（x3（x3（36x）0.5470(x- 25)

28、，x 25,36）0 36- 0250.4599（x3（x3（48x）0.5404(x- 36)，x36 ,48）0 48- 0360.4633（x30( x48)3（64x）0.5333(x- 48)，x 48,64 ）0 640.4689下面進行畫圖，在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(拉格朗日插值)優(yōu)質(zhì)范文.%下面是求拉