西北農(nóng)林科技大學(xué)數(shù)值分析數(shù)值法實(shí)驗(yàn)報(bào)告

1、. . . .數(shù)值法實(shí)驗(yàn)報(bào)告專業(yè)班級(jí)：信息與計(jì)算科學(xué) 121 姓名：金輝 學(xué)號(hào)：20120142801）實(shí)驗(yàn)?zāi)康谋敬螌?shí)驗(yàn)的目的是熟練數(shù)值分析第二章“插值法”的相關(guān)內(nèi)容，掌握三種插值方法：牛頓多項(xiàng)式插值，三次樣條插值，拉格朗日插值，并比較三種插值方法的優(yōu)劣。本次試驗(yàn)要求編寫(xiě)牛頓多項(xiàng)式插值，三次樣條插值，拉格朗日插值的程序編碼，并在 MATLAB軟件中去實(shí)現(xiàn)。2）實(shí)驗(yàn)題目實(shí)驗(yàn)一：已知函數(shù)在下列各點(diǎn)的值為xi 0.2 0.4 0.6 .0.8 1.0f （xi ） 0.98 0.92 0.81 0.64 0.38試用 4 次牛頓插值多項(xiàng)式 P4（x）及三次樣條函數(shù) S（x）（自然邊界條件）對(duì)數(shù)據(jù)進(jìn)行

2、插值。用圖給出 （xi ，yi ），xi =0.2+0.08i ，i=0 ，1, 11, 10 ，P4（x）及 S（x）。實(shí)驗(yàn)二：在區(qū)間-1,1 上分別取 n 10,20 用兩組等距節(jié)點(diǎn)對(duì)龍格函數(shù)f (x)121 25x作多項(xiàng)式插值及三次樣條插值，對(duì)每個(gè) n 值，分別畫(huà)出插值函數(shù)即 f (x) 的圖形。實(shí)驗(yàn)三：下列數(shù)據(jù)點(diǎn)的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函數(shù)的近似，在區(qū)間 0,64 上作圖。（1）用這 9 各點(diǎn)作 8 次多項(xiàng)式插值 L8(x).（2）用三次樣條（自然邊界條件）程序求 S（x）。從結(jié)果看在 0,64 上，那個(gè)

3、插值更精確；在區(qū)間 0,1 上，兩種哪個(gè)更精確？3）實(shí)驗(yàn)原理與理論基礎(chǔ)數(shù)值分析第二章“插值法”的相關(guān)內(nèi)容，包括：牛頓多項(xiàng)式插值，三次樣條插值，拉格朗日4）實(shí)驗(yàn)內(nèi)容實(shí)驗(yàn)一：學(xué)習(xí)資料. . . .已知函數(shù)在下列各點(diǎn)的值為xi 0.2 0.4 0.6 .0.8 1.0f （xi ） 0.98 0.92 0.81 0.64 0.38試用 4 次牛頓插值多項(xiàng)式 P4（x）及三次樣條函數(shù) S（x）（自然邊界條件）對(duì)數(shù)據(jù)進(jìn)行插值。用圖給出 （xi ，yi ），xi =0.2+0.08i ，i=0 ，1, 11, 10 ，P4（x）及 S（x）。(1) 首先我們先求牛頓插值多項(xiàng)式， 此處要用 4 次牛頓插值多

4、項(xiàng)式處理數(shù)據(jù)。已知 n 次牛頓插值多項(xiàng)式如下：Pn=f(x 0)+fx 0,x 1(x-x 0)+ fx 0,x 1,x 2(x-x 0) (x-x 1)+···+fx 0,x 1，··· xn(x-x 0) ···(x-x n-1)我們要知道牛頓插值多項(xiàng)式的系數(shù)，即均差表中得部分均差。在MATLAB的 Editor 中輸入程序代碼， 計(jì)算牛頓插值中多項(xiàng)式系數(shù)的程序如下：function varargout=newtonliu(varargin)clear,clcx=0.2 0.4 0.6 0.8 1.0

5、;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:nfor j=i:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-i+1);endendfor i=1:nfprintf('%4.2f',x(i);for j=1:ifprintf('%10.5f',FF(i,j);en

6、d學(xué)習(xí)資料. . . .fprintf('n');end由 MATLAB計(jì)算得：xi f(x i ) 一階差商 二階差商 三階差商 四階差商0.20 0.9800.40 0.920 -0.300000.60 0.810 -0.55000 -0.625000.80 0.640 -0.85000 -0.75000 -0.208331.0 0.380 -1.30000 -1.12500 -0.62500 -0.52083所以有四次插值牛頓多項(xiàng)式為：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

7、-0.6)-0.52083 (x-0.2)(x-0.4)(x-0.6)(x-0.8 ）(2) 接下來(lái)我們求三次樣條插值函數(shù)。用三次樣條插值函數(shù)由上題分析知 , 要求各點(diǎn)的 M值：M2 0 0 0 0 00M0.500 2 0.500 0 0 - 3.75001000.500 020.5000.500 200.500MM23- 4.50000 0 0 0 2 0M4三次樣條插值函數(shù)計(jì)算的程序如下：function tgsanci(n,s,t) %n 代表元素?cái)?shù)，s,t 代表端點(diǎn)的一階導(dǎo)。x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38;n=5for

8、 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-1學(xué)習(xí)資料. . . .f(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-1a(j+1,j)=u(j+1);a(j,j+1)=r(j)

9、;endb=inv(a)m=b*d'p=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);endpT 得到 m=(0 -1.6071 -1.0714 -3.1071 0)即 M0=0 ;M1= -1.6071;M 2= -1.0714; M 3= -3.1071; M 4=0則根據(jù)三次樣條函數(shù)定義，可得：（0 0.43x）1

10、.3393（x0.2 ）4.9 0（00.4x 4.6536( x 0 .2)，x） 0.2,0.4S(x)=- 1.3393（0.6-0.892 （9 0.83x）3x -）0 x.8929 （2.589 （3 x3-0.4 ）3- 0.6）4.6536 0.6 （4.85 （7 0.8x）x）4.0 857 x （3.303 （6 x0.4 ），x0.6 x）， 0.4,0.6 0.6,0.8-2.1.（ 03x）-（0 x30.8）3.3036 （1.0 x）1.9 ( x0. 8) x， 0.8,1. 0學(xué)習(xí)資料. . . .接著，在 Command Window里輸入畫(huà)圖的程序代碼，

11、下面是畫(huà)牛頓插值以及三次樣條插值圖形的程序：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 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 11x0=0.2+0.08*kfor i=1:1:4y0(i)=0.98-0.3*(x

12、(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,'r')hold ongtext(' 三次樣條自然邊界

13、 ')gtext(' 原圖像')gtext('4 次牛頓插值 ')運(yùn)行上述程序可知：給出的 （xi ，yi ），xi =0.2+0.08i ，i=0 ，1, 11, 10 點(diǎn)，S（x）及 P4（x）圖形如下所示：學(xué)習(xí)資料. . . .10.94 次 牛 頓 插 值0.8三 次 樣 條 自 然 邊 界0.70.6原 圖 像0.50.40.30.20.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2實(shí)驗(yàn)二：在區(qū)間-1,1 上分別取 n 10,20 用兩組等距節(jié)點(diǎn)對(duì)龍格函數(shù)f (x)121 25x作多項(xiàng)式插值及三次樣條插值，對(duì)每個(gè)

14、 n 值，分別畫(huà)出插值函數(shù)即 f (x) 的圖形。我們先求多項(xiàng)式插值：在 MATLAB的 Editor 中建立一個(gè)多項(xiàng)式的 M-file, 輸入如下的命令 （如牛頓插值公式）：function C,D=newpoly(X,Y)n=length(X);D=zeros(n,n)D(:,1)=Y'for 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);學(xué)習(xí)資料. . .

15、 .end當(dāng) n=10 時(shí)，我們?cè)?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,Y,'.');grid on;xp=-1:0.2:1;z=1./(1+25*xp.2);plot(xp,z,'r')得到插值函數(shù)和 f （x）圖形：當(dāng) n=20 時(shí)，我們?cè)?Command Window中輸入以下的命令：clear,clf,hold on;X=-1:0.1:1;Y=1.

16、/(1+25*X.2);C,D=newpoly(X,Y);x=-1:0.01:1;學(xué)習(xí)資料. . . .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）圖形：下面再求三次樣條插值函數(shù)， 在 MATLAB的 Editor 中建立一個(gè)多項(xiàng)式的 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

17、);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);學(xué)習(xí)資料. . . .for k=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);endM(N)=U(N-1)/B(N-1);for k=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-

18、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當(dāng) n=10 時(shí)，我們?cè)?Command Window中輸入以下的命令：clear,clcX=-1:0.2:1;Y=1./(25*X.2+1);dx0= 0.0739644970414201;dxn= -0.0739644970414201;S=csfit(X,Y,dx0,d

19、xn)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,y4, X,Y,'.')結(jié)果如圖：學(xué)習(xí)資料. . . .當(dāng) n=20 時(shí)，我們?cè)?Command Window中輸入以下的命令：clear,clcX=-1:0.1:1;Y=1./(25*X.2+1)

20、;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,y4, X,Y,'.')結(jié)果如圖：學(xué)習(xí)資料. . . .實(shí)驗(yàn)三：

21、下列數(shù)據(jù)點(diǎn)的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函數(shù)的近似，在區(qū)間 0,64 上作圖。（1）用這 9 各點(diǎn)作 8 次多項(xiàng)式插值 L8(x).（2）用三次樣條（自然邊界條件）程序求 S（x）。從結(jié)果看在 0,64 上，那個(gè)插值更精確；在區(qū)間 0,1 上，兩種哪個(gè)更精確？nL8(x) 可由公式 Ln(x)=y k l ( xj ) 得出。kk 0三次樣條可以利用自然邊界條件。寫(xiě)成矩陣：20 0 0 M d0002 0 M01 1 1d10 2 M d02 2 2 20 0 2 M d33 3 320 0 0 M d4 4 4其中

22、j =hhj-1j-1hj， i=hhj-1jhj，dj=6fx j-1 ,x j ,x j+1 ， n= 0=0 d 0=dn=0(x - 1)( x 4 )( x 9 )( x 16 )( x 25 )( x 36 )( x 49 )( x 64 )l 0(x)=(0 1)( 0 4 )( 0 9)( 0 16 )( 0 25 )( 0 36 )( 0 49 )( 0 64 )l 1(x)=(x - 0)( x 4)( x 9)( x 16 )( x 25)( x 36 )( x 49 )( x(1 0)(1 4)(1 9)( 1 16)( 1 25)(1 36 )(1 49 )(164)

23、64)l 2(x)=(x - 0)( x 1)( x 9)( x 16 )( x 25 )( x 36 )( x 49)( x(4 0)( 4 1)( 4 9)( 4 16)( 4 25)( 4 36)( 4 49)( 464)64)l 3(x)=(x - 0)( x 1)( x 4)( x 16 )( x 25 )( x 36 )( x 49 )( x(9 0)(9 1)( 9 4)( 9 16)( 9 25)( 9 36)( 9 49)( 964)64)l 4(x)=(16(x - 0)( x 1)( x 4)( x 9)( x 25 )( x 36 )( x 49 )( x64)0)(1

24、6 1)(16 4 )(16 9)(16 25)(16 36)(16 49)( 1664)l 5(x)=( 25(x - 0) (x - 1)( x 4)( x 9)( x 16 )( x 36 )( x 49 )( x0)( 25 1)( 25 4)( 25 9)( 25 16 )( 25 36 )( 25 49 )(64)2564)l 6(x)=(36(x - 0)( x 1)( x 4)( x 9)( x 16)( x 25 )( x 49 )( x0)( 36 1)( 36 4)( 36 9)( 36 16)( 36 25)( 36 4964)( 36)64)l 7(x)=（49（x

25、0 (x - 1)( x 4)( x 9)( x 16 )( x 25)( x 36 )( x- ）0）(49 1)( 49 4)( 49 9)( 49 16 )( 49 25) 49 36 )(（64)4964)l 8(x)=(64( x - 0)(x - 1)( x 4)( x 9 )( x 16 )( x 25 )( x 36 )( x0 )( 64 1)( 64 4 )( 64 9 )( 64 16 )( 64 25 )( 64 36)(49)6449 )L8(x)= l 1(x)+2 l 2(x)+3 l 3(x)+4 l 4(x)+5 l 5(x)+6 l 6(x)+7 l 7(x

26、)+8 l 8(x)學(xué)習(xí)資料. . . .求三次樣條插值函數(shù)由 MATLAB計(jì)算：可得矩陣形式的線性方程組為：2.0 000 0 0 0 0 0 0 0 0 0M00.2 500 2 .0 000 0 .7 500 0 0 0 0 0 0 M 1.010 0 .3 750 2 .0 000 0.6250 0 0 0 0 0 M 0.12M0 0 0.4167 2.0000 0.5833 0 0 0 0 0.02863M0 0 0 0.4375 2.0 000 0.5625 0 0 0 0.01194M0 0 0 0 0.4 500 2.0000 0 .5 500 0 0 0.00615M0 0

27、 0 0 0 0.4583 2 .0 000 0.5 417 0 0.00356M0 0 0 0 0 0 0 .4 634 2.0000 0.5 3 57 0.00227M0 0 0 0 0 0 0 0 2.0 000 08在 MATLAB中的 Editor 中輸入程序代碼，以下是三次樣條函數(shù)的程序代碼：function tgsanci(n,s,t) %n 代表元素?cái)?shù)， s,t 代表端點(diǎn)的一階導(dǎo)。y=0 1 2 3 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

28、(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:n學(xué)習(xí)資料. . . .a(j,j)=2;endr(1)=0;u(n)=0;for j=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

29、(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解得:M0=0;M1=-0.5209;M 2=0.0558;M3=-0.0261;M 4=0.0006;M5=-0.0029;M 6=-0.0008;M 7=-0.0009;M8=0，則三次樣條函數(shù)：（0 13x）0.868 ( x 0)3（0 1 x） 1.0868(x - 0)， x 0,1-0.289 （4 x3

30、）- 0.0031（x-31）0.5938 （4 x 0.6388） （x1），x1, 40.19 9（3x）0.9（x -34）0.3535 （9 x） 0.6271 x（4），x4,9S(x)=0-0.6 16（0 253x）-3x）（0x0.1（x39）316）-0.4590 （160.4436 25（x）- 0.5708(xx 0.5600(x）-9)，- 16)，xx9,1 6 16,25（0 363x）-（0x325）0.4599 （36 x 0.5470(x）-25)，x25,36 （0 483x）-（0x336）0.463 （3 48 x）0.5404 (x-36)，x36,4

31、8（0 643x）0( x 48)30.4689 （64x 0.5333(x）- 48) x ，48,64下面進(jìn)行畫(huà)圖，在 Command Window中輸入畫(huà)圖的程序代碼：%畫(huà)圖形比較那個(gè)插值更精確的函數(shù)：x0=0 1 4 9 16 25 36 49 64;y0=0 1 2 3 4 5 6 7 8;學(xué)習(xí)資料. . . .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&#

32、39;);hold on;plot(x0,y0,':b');hold on%axis(0 2 0 1); % 看0 1 區(qū)間的圖形時(shí)加上這條指令gtext(' 三次樣條插值 ')gtext(' 原圖像')gtext(' 拉格朗日插值 ')%下面是求拉格朗日插值的函數(shù)function y=lagr1(x0,y0,x)n=length(x0); m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s;end拉格朗日插值函數(shù)與三次樣條插值函數(shù)如圖中所示， 綠色實(shí)