數(shù)值分析實(shí)驗(yàn)

1、數(shù)值分析計(jì)算實(shí)習(xí)報(bào)告冊(cè) 專業(yè) 信息與計(jì)算科學(xué)學(xué)號(hào) 20112984姓名 張妃_ 20122013年第一學(xué)期實(shí)驗(yàn)一 數(shù)值計(jì)算的工具 Matlab1.解釋下MATLAB程序的輸出結(jié)果程序：t=0.1n=1:10e=n/10-n*t答案：t = 0.1000 小數(shù)點(diǎn)后輸出四位n = 1 2 3 4 5 6 7 8 9 10 n取1到10e =1.0e-015 * 0 0 -0.0555 0 0 -0.1110 -0.1110 0 0 0 e等于-0.0555×2.下面MATLAB程序的的功能是什么？程序：x=1;while 1+x>1,x=x/2,pause(0.02),endx=1

2、;while x+x>x,x=2*x,pause(0.02),endx=1;while x+x>x,x=x/2,pause(0.02),end答案：當(dāng)x=1時(shí)，對(duì)x無(wú)限二分，到1+x>1結(jié)束，且循環(huán)之間間隔0.02秒當(dāng)x=1時(shí)，對(duì)無(wú)限循環(huán)，到x+x>x結(jié)束，且循環(huán)之間間隔0.02秒當(dāng)x=1時(shí)，對(duì)x無(wú)限二分，到x+x>x結(jié)束，且循環(huán)之間間隔0.02秒3.考慮下面二次代數(shù)方程的求解問(wèn)題 公式是熟知的，與之等價(jià)地有，對(duì)于，應(yīng)當(dāng)如何選擇算法。答案：clc,cleara=1;b=100000000;c=1;x1=(-b+sqrt(b*b-4*a*c)/2*ax2=(-b+s

3、qrt(b*b-4*a*c)/2*ax3=2*c/(-b+sqrt(b*b-4*a*c)x4=2*c/(-b+sqrt(b*b-4*a*c)d=x1*x1+100000000*x1+1e=x2*x2+100000000*x2+1f=x3*x3+100000000*x3+1g=x4*x4+100000000*x4+1x1 =-7.4506e-009=-7.4506×x2 = -7.4506e-009x3 =-134217728x4 = -134217728d =0.25494194030762e =0.25494194030762f = 4.592625709481985e+015g

4、= 4.592625709481985e+015因?yàn)閐，e比f(wàn)，g更接近0所以公式更好。4.函數(shù)有冪級(jí)數(shù)展開(kāi)利用冪級(jí)數(shù)計(jì)算的MATLAB程序?yàn)閒unction s=powersin(x)s=0;t=x;n=1;while s+t=s; s=s+t； t=-x2/（(n+1)*(n+2)）*t； n=n+2；end修改后程序function s=powersin(x)s=0;x=pi/2;t=x;k=0;n=1;while s+t=s; s=s+t; t=-x2/(n+1)*(n+2)*t; n=n+2; k=k+1;endK結(jié)果X=21*pi/2k =60ans = 0.9999X=11*pi

5、/2k = 37ans = -1.0000X=pi/2k =11ans = 1.0000(a) 解釋上述程序的終止準(zhǔn)則。t=0時(shí)，終止(b)對(duì)于計(jì)算的進(jìn)度是多少？分別計(jì)算多少項(xiàng)？ ans =1.0000 ans =-1.0000 ans = 0.9999 11,37,605.考慮調(diào)和級(jí)數(shù)，它是微積分中的發(fā)散級(jí)數(shù)，在計(jì)算機(jī)上計(jì)算該級(jí)數(shù)的部分和，會(huì)得到怎么樣的結(jié)果，為什么？clc;clears=0;n=20; t=1/i; s=s+t;endS結(jié)果s = 18.5697for i=n-1:n+100000s = 8.5952s = 18.5697S會(huì)隨著i的項(xiàng)數(shù)的增大而增大，到一定程度后不變。因?yàn)?/p>

6、在計(jì)算機(jī)中如果一個(gè)數(shù)很小，計(jì)算機(jī)會(huì)默認(rèn)為0.6.指數(shù)函數(shù)的級(jí)數(shù)展開(kāi)，如果對(duì)于，用上述的級(jí)數(shù)近似計(jì)算指數(shù)函數(shù)的值，這樣的算法結(jié)果是否會(huì)好，為什么？1.clc,clears=1;t=1;x=-1;for i=1:100000 t=t*i; tt=xi/t; s=s+tt;ends X=1;s=2.7183 x=-1; s= 0.36792.function s=powersin(x)s=0;m=1;t=1;n=1;while s+t=s; s=s+t; m=m*n; t=-x(n)/m; n=n+2;endx=0,f=1;x=1,y= -0.4107;x=-1,f=2,41077.考慮數(shù)列，它的統(tǒng)

7、計(jì)平均定義為 ，標(biāo)準(zhǔn)差數(shù)學(xué)上等價(jià)于作為標(biāo)準(zhǔn)差的兩種算法，你將如何評(píng)價(jià)他們的得與失。1. clc,cleara=input('請(qǐng)輸入數(shù)據(jù)')n=length(a);sm=sum(a);ave=sm/n;b1=0;for i=1:n b1=(a(i)-ave)2+b1;endb1b2=(1/(n-1)*b1)(1/2)結(jié)果b1=133.9368b2=11.57312.clc,cleara=input('請(qǐng)輸入數(shù)據(jù)')n=length(a);sm=sum(a);ave=sm/n;b3=0;for i=1:n b3=(a(i)2+b3;endb3b1=(1/(n-1)*

8、(b3-n*(ave2)(1/2)結(jié)果b3=189545;b1=11.8630實(shí)驗(yàn)二 插值法計(jì)算實(shí)習(xí)題1. 已知函數(shù)在下列個(gè)點(diǎn)的值為0.20.40.60.81.00.980.920.810.640.38試用4次插值多項(xiàng)式及三次樣條插值(自然邊界條件)對(duì)數(shù)據(jù)進(jìn)行插值。用圖給出，及。一對(duì)于多項(xiàng)式插值：程序1：（第一種編程）（1）建立function函數(shù)functionc, d=poly1(x, y) n=length(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) /

9、(x(k)-x(k-j+1); endend c =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（2)建立一個(gè)主程序xunan1.mclcclearpoly1(0.2,0.4,0.6,0.8,1.0, 0.98, 0.92, 0.81, 0.64, 0.38)結(jié)果：ans = -0.5208 0.8333 -1.1042 0.1917 0.9800得出所求的牛頓多項(xiàng)式為： 程序2:(第二種編程）x=0.2,0.4,0.6,0.8,1;y=0.98,0.92,0.8

10、1,0.64,0.38;a=polyfit(x,y,length(x)-1);%插值poly2sym(a)%輸出插值多項(xiàng)式結(jié)果：ans = 畫(huà)圖：程序1：x=0.2 ;y=-0.5208*x4+0.8333*x3-1.1042*x2+0.1917*x+0.9800x=0.2+0.08i ;y=-0.5208*x4+0.8333*x3-1.1042*x2+0.1917*x+0.9800x=0.2+0.08*11i;y=-0.5208*x4+0.8333*x3-1.1042*x2+0.1917*x+0.9800x=0.2+0.08*10i ;y=-0.5208*x4+0.8333*x3-1.104

11、2*x2+0.1917*x+0.9800結(jié)果：y =0.9800y = 0.9847 - 0.0135iy =1.2324 - 0.4306iy = 1.2334 - 0.3466i程序2：x=0.2 0.2+0.08i 0.2+0.08*11i 0.2+0.08*10i;y=0.9800 0.9847-0.0135i 1.2324-0.4306i 1.2334-0.3466iplot(x,y,'-r+')結(jié)果：二對(duì)于樣條插值：程序1：x = 0.2:0.2:1.0;y= 0.98, 0.92, 0.81, 0.64, 0.38;xx = linspace(0.2,1.0,21

12、)cs = spline(x,y,xx)結(jié)果：xx = Columns 1 through 12 0.2000 0.2400 0.2800 0.3200 0.3600 0.4000 0.4400 0.4800 0.5200 0.5600 0.6000 0.6400 Columns 13 through 21 0.6800 0.7200 0.7600 0.8000 0.8400 0.8800 0.9200 0.9600 1.0000cs = Columns 1 through 12 0.9800 0.9718 0.9617 0.9497 0.9358 0.9200 0.9022 0.8823 0

13、.8603 0.8362 0.8100 0.7815Columns 13 through 21 0.7506 0.7168 0.6801 0.6400 0.5963 0.5488 0.4970 0.4409 0.380clearclcx = 0.2:0.2:1.0;y= 0.98, 0.92, 0.81, 0.64, 0.38;cs = spline(x,0 y 0)xx = linspace(0.2,1.0,21);plot(x,y,'o',xx,ppval(cs,xx),'-')畫(huà)圖：程序：clearclcx = 0.2:0.2:1.0;y= 0.98, 0

14、.92, 0.81, 0.64, 0.38;cs = spline(x,0 y 0)xx =0.2 0.2+0.08i 0.2+0.08*11i 0.2+0.08*10i;plot(x,y,'o',xx,ppval(cs,xx),'-')2. 在區(qū)間-1,1上分別取用兩組等距節(jié)點(diǎn)對(duì)龍格函數(shù)作多項(xiàng)式插值及三次樣 條插值，對(duì)每個(gè)值，分別畫(huà)出插值函數(shù)及的圖形。一當(dāng)時(shí)的程序：結(jié)果：a=-1;b=1;n=10;for i=1:11x(i)=a+(i-1)*(b-a)/n;y(i)=1/(1+25*x(i)*x(i);endplot(x,y,'r-')hol

15、d onz=a:(b-a)/(2*n):b;n=length(x);for j=2:n for i=n:-1:j y(i)=(y(i)-y(i-1)/(x(i)-x(i-j+1); endendu=y(n);m=length(z);for j=1:m for i=n-1:-1:1 u=y(i)+u*(z(j)-x(i); v(j)=u; end u=y(n);endplot(z,v,'b')hold off二當(dāng)時(shí)的程序：a=-1;b=1;n=20for i=1:21x(i)=a+(i-1)*(b-a)/n;y(i)=1/(1+25*x(i)*x(i);endplot(x,y,&

16、#39;r-')hold onz=a:(b-a)/(2*n):b;n=length(x);for j=2:n for i=n:-1:j y(i)=(y(i)-y(i-1)/(x(i)-x(i-j+1); endendu=y(n);m=length(z);for j=1:m for i=n-1:-1:1 u=y(i)+u*(z(j)-x(i); v(j)=u; end u=y(n);endplot(z,v,'b')hold off結(jié)果：3. 下列數(shù)據(jù)點(diǎn)的插值x01491625364964y012345678可以得到平方根函數(shù)的近似值，在區(qū)間0,64上作圖。（1） 用這九個(gè)

17、點(diǎn)作8次多項(xiàng)式插值 .（2） 用三次樣條（第一邊界條件）程序求從得到結(jié)果看在0,64上，哪個(gè)插值更精確；在區(qū)間0,10上，兩種插值哪個(gè)更精確？一多項(xiàng)式插值作圖：程序：x=0 1 4 9 16 25 36 49 64;y=0:1:8;p=polyfit(x,y,8)X=0:1:64;Y=polyval(p,X);plot(x,y,'r-',X,Y,'y-')title('8次多項(xiàng)式')xlabel('x')ylabel('y')1.建立function函數(shù)function f=nan(x,y,x0)syms t;if

18、(length(x)=length(y) n=length(x);else disp('x!=y') return;endf=0;for(i=1:n) l=y(i); for(j=1:i-1) l=l*(t-x(j)/(x(i)-x(j); end f=f+l; if(i=n) if(nargin=3) f=subs(f,'t',x0); else f=collect(f); f=vpa(f,6); end endend2.建立一個(gè)主程序xunan1.mclearclcx=0 1 4 9 16 25 36 49 64;y=0:1:8;f=nan(x,y)結(jié)果：f

19、 = .764716e-12*t8+.535301e-10*t7+.161447e-7*t6+.171358e-5*t5+.136120e-3*t4+.615224e-2*t3+.132896*t2+.860814*t（2）三次樣條插值：程序：%得到三次樣條擬合函數(shù)%得到三次樣條擬合函數(shù)syms x sx1=0 1 4 9 16 25 36 49 64;y1=0 1 2 3 4 5 6 7 8;x2=0:1:64;y2=spline(x1,y1,x2);p=polyfit(x2,y2,3) S=p(1)+p(2)*x+p(3)*x2+p(4)*x3 plot(x2,y2)hold onplot

20、(x1,y1,'r-')結(jié)果：p = 0.0000 -0.0042 0.2527 0.8876 S =所以在0,64區(qū)間上樣條插值更精確。二當(dāng)0,10區(qū)間上x(chóng)=0 1 4 9 16 25 36 49 64;y=0:1:8;p=polyfit(x,y,8)X=0:1:10;Y=polyval(p,X);plot(x,y,'r-',X,Y,'y-')title('8次多項(xiàng)式')xlabel('x')ylabel('y')結(jié)果：程序：syms x sx1=0 1 4 9 16 25 36 49 64;y1

21、=0 1 2 3 4 5 6 7 8;x2=0:1:10;y2=spline(x1,y1,x2);p=polyfit(x2,y2,3) S=p(1)+p(2)*x+p(3)*x2+p(4)*x3 plot(x2,y2)hold onplot(x1,y1,'r-')title('三次樣條')結(jié)果：由圖可知在0,10區(qū)間上多項(xiàng)式插值更精確。實(shí)驗(yàn)三 函數(shù)逼近與快速傅立葉變換 1. 對(duì)于給函數(shù)在區(qū)間-1,1上去 ，試求3次曲線擬合，試畫(huà)出擬合曲線并打印方程，與實(shí)驗(yàn)二，第二章的計(jì)算實(shí)習(xí)題2的結(jié)果比較。clc,clearn=11;for i=1:n x(i)=-1+0.2*

22、(i-1); y(i)=1/(1+25*x(i)*x(i);endA=polyfit(x,y,3)z=polyval(A,x)plot(x,y,'k+',x,z,'r')結(jié)果：A = 0.0000 -0.5752 -0.0000 0.4841f(x)= -0.5752*X2+0.48412. 由實(shí)驗(yàn)給出數(shù)據(jù)表 00.10.20.30.50.81 1.00.410.500.610.912.022.46試求3次，4次多項(xiàng)式的曲線擬合，再根據(jù)數(shù)據(jù)曲線圖形，求一個(gè)另外函數(shù)的擬合曲線，用圖示數(shù)據(jù)曲線及相應(yīng)的三種擬合曲線。clc,clearx=0.0 0.1 0.2 0.3

23、 0.5 0.8 1.0y=1.0 0.41 0.50 0.61 0.91 2.02 2.46A=polyfit(x,y,3)z=polyval(A,x)plot(x,y,'k+',x,z,'r')結(jié)果：A = -6.6221 12.8147 -4.6591 0.9266f(x)=-6.6221*x3+12.8147*X2+-4.6591*x+0.92664次:clc,clearx=0.0 0.1 0.2 0.3 0.5 0.8 1.0y=1.0 0.41 0.50 0.61 0.91 2.02 2.46A=polyfit(x,y,4)z=polyval(A,x

24、)plot(x,y,'k+',x,z,'r')結(jié)果：A = 2.8853 -12.3348 16.2747 -5.2987 0.9427f(x)=2.8853*x4-12.3348*x3+16.2747*X2+-5.2987*x+0.9427根據(jù)數(shù)據(jù)曲線形狀，我覺(jué)得它像是一個(gè)開(kāi)口向上的二次函數(shù)的圖像.clc,clearx=0.0 0.1 0.2 0.3 0.5 0.8 1.0y=1.0 0.41 0.50 0.61 0.91 2.02 2.46A=polyfit(x,y,2)z=polyval(A,x)plot(x,y,'k+',x,z,'

25、;r')結(jié)果：A =3.1316 -1.2400 0.7356f(x)= 3.1316*X2+-1.2400*x+0.7356實(shí)驗(yàn)五 數(shù)值微分 精度及步長(zhǎng)的關(guān)系實(shí)驗(yàn)?zāi)康模簲?shù)值計(jì)算中誤差是不可避免的，希望同學(xué)能通過(guò)本實(shí)驗(yàn)初步認(rèn)識(shí)數(shù)值分析中極端重要的兩個(gè)概念：截?cái)嗾`差和舍入誤差，并認(rèn)真體會(huì)誤差對(duì)計(jì)算結(jié)果的影響.問(wèn)題提出：設(shè)一元函數(shù) ，則 在 的導(dǎo)數(shù)定義為實(shí)驗(yàn)內(nèi)容：計(jì)算在的導(dǎo)數(shù)值可以用如下的差分給出其近似： （1） （2）這也就給出了計(jì)算函數(shù)導(dǎo)數(shù)的兩個(gè)簡(jiǎn)單算法.請(qǐng)給出幾個(gè)計(jì)算高階導(dǎo)數(shù)的近似算法，并完成如下工作：1. 對(duì)同樣的h，比較（1）和（2）的計(jì)算結(jié)果.2. 針對(duì)計(jì)算高階導(dǎo)數(shù)的算法，比較

26、h取不同值時(shí)的計(jì)算結(jié)果.實(shí)驗(yàn)要求：選擇有代表性的函數(shù) （最好選擇多個(gè)），利用MATLAB提供的繪圖工具畫(huà)出該函數(shù)在某個(gè)區(qū)間的導(dǎo)數(shù)曲線 .再將數(shù)值計(jì)算的結(jié)果用MATLAB畫(huà)出來(lái)，認(rèn)真思考實(shí)驗(yàn)的結(jié)果.同學(xué)們還可以圍繞這一問(wèn)題設(shè)計(jì)其它的實(shí)驗(yàn)內(nèi)容，特別認(rèn)真的體會(huì)導(dǎo)數(shù)的結(jié)果。對(duì)于問(wèn)題1：我們分別用算法（1）、（2）計(jì)算函數(shù)在上導(dǎo)數(shù),并與在的值進(jìn)行對(duì)比，畫(huà)出圖像（如下圖1所示）：由圖像可知：其差距并不明顯.于是我們便計(jì)算出分別用算法（1）、（2）計(jì)算函數(shù)在上導(dǎo)數(shù)的誤差，并畫(huà)出其圖像，（如下圖2所示，其中藍(lán)色表示用算法（1）計(jì)算得到的誤差，紅色表示用算法（2）計(jì)算得到的誤差）： 圖1圖2由圖像可知：對(duì)于計(jì)算

27、函數(shù)導(dǎo)數(shù)，(2)是比(1)更好的算法.對(duì)于問(wèn)題1實(shí)驗(yàn)結(jié)論的分析：不論采用怎樣的算法，計(jì)算結(jié)果通常都會(huì)有誤差。比如算法(1)，由Taylor公式，知： 所以有 利用上式來(lái)計(jì)算，其截?cái)嗾`差為： 所以誤差是存在的，并且當(dāng)步長(zhǎng)h越來(lái)越小時(shí)，(1)的近似程度也越來(lái)越好。類似地可以分析(2)的截?cái)嗾`差為：上述截?cái)嗾`差的分析表明(2)是比(1)更好的算法，因?yàn)閷?duì)步長(zhǎng)h(<<1)，(2)比(1)更接近于.關(guān)于問(wèn)題1的程序如下：function y=dsh1(fu,x,h)y=(feval(fu,x+h)-feval(fu,x)/h;y;function y=dsh2(fu,x,h)y=(feval

28、(fu,x+h)-feval(fu,x-h)/(2*h);y;clear,clc;x=-pi:0.1:pi;y=cos(x);y1=dsh1('sin',x,0.01);y2=dsh2('sin',x,0.01);plot(x,y,'b*',x,y1,'ro',x,y2,'go')pause(1);wucha1=abs(y1-y);wucha2=abs(y2-y);plot(x,wucha1,'b*',x,wucha2,'ro')對(duì)于問(wèn)題2： 我們用如下算法來(lái)計(jì)算二階導(dǎo)數(shù)：該算法程序

29、如下：function y=dsh3(fu,x,h)y=(feval(fu,x+2*h)-feval(fu,x)-feval(fu,x+h)+feval(fu,x-h)/(2*h*h);y;當(dāng)h分別為0.01,0.001及0.000001時(shí)，通過(guò)該算法計(jì)算出此時(shí)sin(x)在上的二階導(dǎo)數(shù)，并畫(huà)出圖像，與在上的圖像進(jìn)行對(duì)比，其圖像如下圖所示：圖3有圖像可看出其差別并不大。于是我們做出誤差誰(shuí)的變化圖像，如下圖所示：圖4由圖像我們不難發(fā)現(xiàn)：當(dāng)h小到一定程度后，計(jì)算結(jié)果的誤差不但不再減少，反而會(huì)變大。對(duì)于問(wèn)題2實(shí)驗(yàn)結(jié)論的分析：事實(shí)上當(dāng)步長(zhǎng)h過(guò)小時(shí)，計(jì)算結(jié)果的誤差變大就是由于舍入誤差的緣故。關(guān)于問(wèn)題2的

30、程序如下：clear,clc;x=-pi:0.1:pi;y=-sin(x);h=0.01;y1=dsh3('sin',x,h);h=0.001;y2=dsh3('sin',x,h);h=0.000001;y3=dsh3('sin',x,h);plot(x,y,'b*',x,y1,'r*',x,y2,'go',x,y3,'yo')pause(1);h=;for i=7:-1:2 a=10(-i):(10(-i+1)-10(-i)/9:10(-i+1); h=h,a;endwucha=;for i=1:6