西南交大數(shù)值分析2013上機(jī)作業(yè)

1、數(shù)值分析上機(jī)作業(yè)姓 名：XXX 學(xué) 號(hào)：XXX 專 業(yè)：XXX 聯(lián)系電話：XXX 序言長(zhǎng)期以來(lái)，在數(shù)值分析教學(xué)中，數(shù)值分析作為一門純粹的理論課程進(jìn)行授課，但結(jié)果便忽略數(shù)值分析上機(jī)的試驗(yàn)環(huán)節(jié)，忽略了將編程語(yǔ)言與教材中給出的算法結(jié)合起來(lái)，以解決實(shí)際問(wèn)題。而對(duì)于C語(yǔ)言和Fortran語(yǔ)言，只有基本熟練掌握后才具有一定得編程能力，而且此兩種語(yǔ)言編程效率較低，不適于數(shù)值分析編程。本次數(shù)值分析上機(jī)實(shí)習(xí)所用的語(yǔ)言是MATLAB語(yǔ)言。MATLAB語(yǔ)言是著名的矩陣計(jì)算數(shù)學(xué)軟件，除具備卓越的數(shù)值計(jì)算能力外，它還提供了專業(yè)水平的符號(hào)計(jì)算，文字處理，可視化建模仿真和實(shí)時(shí)控制等功能。MATLAB語(yǔ)言主要特點(diǎn)：1）高效的

2、數(shù)值計(jì)算及符號(hào)計(jì)算功能，能使用戶從繁雜的數(shù)學(xué)運(yùn)算分析中解脫出來(lái)；2）具有完備的圖形處理功能，實(shí)現(xiàn)計(jì)算結(jié)果和編程的可視化；3）友好的用戶界面及接近數(shù)學(xué)表達(dá)式的自然化語(yǔ)言，使學(xué)者易于學(xué)習(xí)和掌握；4）功能豐富的應(yīng)用工具箱(如信號(hào)處理工具箱、通信工具箱等) ，為用戶提供了大量方便實(shí)用的處理工具。編程環(huán)境:MATLAB由一系列工具組成。這些工具方便用戶使用MATLAB的函數(shù)和文件，其中許多工具采用的是圖形用戶界面。簡(jiǎn)單的編程環(huán)境提供了比較完備的調(diào)試系統(tǒng)，程序不必經(jīng)過(guò)編譯就可以直接運(yùn)行，而且能夠及時(shí)地報(bào)告出現(xiàn)的錯(cuò)誤及進(jìn)行出錯(cuò)原因分析。簡(jiǎn)單易用：MATLAB是一個(gè)高級(jí)的矩陣/陣列語(yǔ)言，它包含控制語(yǔ)句、函數(shù)、

3、數(shù)據(jù)結(jié)構(gòu)、輸入和輸出和面向?qū)ο缶幊烫攸c(diǎn)。用戶可以在命令窗口中將輸入語(yǔ)句與執(zhí)行命令同步，也可以先編寫好一個(gè)較大的復(fù)雜的應(yīng)用程序（M文件）后再一起運(yùn)行。而且這種語(yǔ)言可移植性好、可拓展性極強(qiáng)。強(qiáng)處理能力：MATLAB是一個(gè)包含大量計(jì)算算法的集合。其擁有600多個(gè)工程中要用到的數(shù)學(xué)運(yùn)算函數(shù)，可以方便的實(shí)現(xiàn)用戶所需的各種計(jì)算功能。函數(shù)中所使用的算法都是科研和工程計(jì)算中的最新研究成果，而且經(jīng)過(guò)了各種優(yōu)化和容錯(cuò)處理。圖形處理：MATLAB自產(chǎn)生之日起就具有方便的數(shù)據(jù)可視化功能，以將向量和矩陣用圖形表現(xiàn)出來(lái)，并且可以對(duì)圖形進(jìn)行標(biāo)注和打印。高層次的作圖包括二維和三維的可視化、圖象處理、動(dòng)畫和表達(dá)式作圖?？捎糜诳?/p>

4、學(xué)計(jì)算和工程繪圖。程序接口：MATLAB可以利用MATLAB編譯器和C/C+數(shù)學(xué)庫(kù)和圖形庫(kù)，將自己的MATLAB程序自動(dòng)轉(zhuǎn)換為獨(dú)立于MATLAB運(yùn)行的C和C+代碼。工具箱是MATLAB函數(shù)的子程序庫(kù)，每一個(gè)工具箱都是為某一類學(xué)科專業(yè)和應(yīng)用而定制的，主要包括信號(hào)處理、控制系統(tǒng)、神經(jīng)網(wǎng)絡(luò)、模糊邏輯、小波分析和系統(tǒng)仿真等方面的應(yīng)用。此數(shù)值分析上機(jī)實(shí)習(xí)報(bào)告包含兩道大題，其中第一題包含兩個(gè)小題，為必做題，第二題為選作題，從五個(gè)題中選做兩個(gè)。每道題分為四個(gè)部分，分別為題目分析、編寫程序、計(jì)算結(jié)果和結(jié)果分析，對(duì)每個(gè)題的不同情況進(jìn)行了一定的分析，做到理論與實(shí)際結(jié)合，以解決實(shí)際問(wèn)題。目錄序言I1.必做題一11.

5、1 Gauss消元法11.1.1 第一小題11.1.2 第二小題2 第三小題21.2高斯塞德爾2 第一小題2 第二小題3 第三小題31.3 計(jì)算結(jié)果對(duì)比分析42. 必做題二52.1 步長(zhǎng)為0.0552.2 步長(zhǎng)為0.162.3 步長(zhǎng)為0.272.4 計(jì)算結(jié)果對(duì)比分析73. 選做題題三93.1 選取初值為293.2 選取初值為393.3 選取初值為493.4 結(jié)果分析104. 選做題題四114.1 復(fù)化梯形公式114.2 復(fù)化Simpson公式114.3 Romberg算法114.4 結(jié)果分析11總結(jié)12附錄131.必做題一寫出對(duì)一般的線性方程組通用的Gauss消元, Gauss-Seidel迭

6、代程序。并以下面的線性方程組為例進(jìn)行計(jì)算，討論所得到的計(jì)算結(jié)果是否與理論一致。（1）（2）（3）1.1 Gauss消元法1.1.1 第一小題1.1.1.1右端項(xiàng)為b1A=6,2,-1;1,4,-2;-3,1,4;b=-3,2,4'x=Gauss(A,b)計(jì)算結(jié)果： x = -0.7273 0.8081 0.25251.1.1.2右端項(xiàng)為b2A=6,2,-1;1,4,-2;-3,1,4;b=100,-200,345'x=Gauss(A,b)計(jì)算結(jié)果： x = 36.3636 -2.0707 114.04041.1.2 第二小題1.1.2.1右端項(xiàng)為b1A=1,0.8,0.8;0.

7、8,1,0.8;0.8,0.8,1;b=3,2,1'x=Gauss(A,b)計(jì)算結(jié)果：x = 5.7692 0.7692 -4.23081.1.2.2右端項(xiàng)為b2A=1,0.8,0.8;0.8,1,0.8;0.8,0.8,1;b=5,0,-10'x=Gauss(A,b)計(jì)算結(jié)果： x = 32.6923 7.6923 -42.30771.1.3第三小題A=1,3;-7,1;b=4,6'x=Gauss(A,b)計(jì)算結(jié)果： x = -0.6364 1.54551.2高斯塞德爾1.2.1第一小題.1右端項(xiàng)為b1A=6,2,-1;1,4,-2;-3,1,4;b=-3,2,4&#

8、39;x=Gauss_seidel(A,b)計(jì)算結(jié)果：x = -0.7273 0.8081 0.25251.2.1.2右端項(xiàng)為b2A=6,2,-1;1,4,-2;-3,1,4;b=100,-200,345'x=Gauss_seidel(A,b)計(jì)算結(jié)果： x = 36.3636 -2.0707 114.04041.2.2第二小題.1右端項(xiàng)為b1A=1,0.8,0.8;0.8,1,0.8;0.8,0.8,1;b=3,2,1'x=Gauss_seidel(A,b)計(jì)算結(jié)果 x = 5.7692 0.7692 -4.23071.2.2.2右端項(xiàng)為b2A=1,0.8,0.8;0.8,1

9、,0.8;0.8,0.8,1;b=5,0,-10'x=Gauss_seidel(A,b)計(jì)算結(jié)果： x = 32.6923 7.6923 -42.30771.2.3第三小題A=1,3;-7,1;b=4,6'x=Gauss_seidel(A,b)計(jì)算結(jié)果： x = 0.1753×10131 1.22681.3 計(jì)算結(jié)果對(duì)比分析表1.1 Gauss消元法與Gauss-Seidel迭代法對(duì)比計(jì)算結(jié)果方法（1）（2）（3）b1b2b1b2bGauss消元x1-0.727336.36365.769232.6923-0.6364x20.8081-2.07070.76927.692

10、31.5455x30.2525114.0404-4.2308-42.3077Gauss-Seidelx1-0.727336.36365.769232.69230.1753×10131x20.8081-2.07070.76927.69231.2268x30.2525114.0404-4.2307-42.3077表1.2 高斯賽德爾迭代法斂散性分析題目序列號(hào)特征值譜半徑斂散性（1）、（2）1= 02= 0.0469 + 0.3504i3= 0.0469 - 0.3504i（B）=0.3536（B）<1收斂（3）、（4）1=02= 0.7040 + 0.1280i3= 0.7040

11、- 0.1280i（B）= 0.7155（B）<1收斂（5）1=02=-21（B）=21（B）>1發(fā)散小結(jié)：從表1.2可知，對(duì)矩陣方程組，利用高斯賽德爾迭代法計(jì)算，當(dāng)矩陣的譜半徑（B）<1時(shí)，結(jié)果收斂；當(dāng)（B）>1時(shí)，結(jié)果發(fā)散。而且（B）值越小，收斂速度越快。計(jì)算結(jié)果收斂與否只取決于迭代矩陣的譜半徑，與初始量及方程組的右端項(xiàng)無(wú)關(guān)，只與迭代矩陣有關(guān)。前兩個(gè)小題高斯-塞德爾迭代方法在b不同時(shí)其計(jì)算結(jié)果與理論值很接近，但第三小題結(jié)果卻與理論值有很大差異，原因是此種情況下高斯-塞德爾迭代法不收斂。通過(guò)計(jì)算可知： =Max=21>1所以此時(shí)高斯-塞德爾迭代算法計(jì)算結(jié)果與理論

12、值差異很大。2.必做題二給定初值問(wèn)題 （精確解為），用Runge-Kutta 4階算法按步長(zhǎng)求解，分析其中遇到的現(xiàn)象及問(wèn)題。2.1 步長(zhǎng)為0.05fun=(x,y)-20*y;x,y=runge_kutta(fun,0,2,1,0.05)計(jì)算結(jié)果：N =40x = Columns 1 through 6 0 0.0500 0.1000 0.1500 0.2000 0.2500 Columns 7 through 12 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 Columns 13 through 18 0.6000 0.6500 0.7000 0.75

13、00 0.8000 0.8500 Columns 19 through 24 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 Columns 25 through 30 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 Columns 31 through 36 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 Columns 37 through 41 1.8000 1.8500 1.9000 1.9500 2.0000y = Columns 1 through 6 1.0000 0.37

14、50 0.1406 0.0527 0.0198 0.0074 Columns 7 through 12 0.0028 0.0010 0.0004 0.0001 0.0001 0.0000 Columns 13 through 18 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 19 through 24 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 25 through 30 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 31 thr

15、ough 36 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 37 through 41 0.0000 0.0000 0.0000 0.0000 0.00002.2 步長(zhǎng)為0.1fun=(x,y)-20*y;x,y=runge_kutta(fun,0,2,1,0.1)計(jì)算結(jié)果：N =20x = Columns 1 through 6 0 0.1000 0.2000 0.3000 0.4000 0.5000 Columns 7 through 12 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 Colum

16、ns 13 through 18 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 Columns 19 through 21 1.8000 1.9000 2.0000y = Columns 1 through 6 1.0000 0.3333 0.1111 0.0370 0.0123 0.0041 Columns 7 through 12 0.0014 0.0005 0.0002 0.0001 0.0000 0.0000 Columns 13 through 18 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Column

17、s 19 through 21 0.0000 0.0000 0.00002.3 步長(zhǎng)為0.2fun=(x,y)-20*y;x,y=runge_kutta(fun,0,2,1,0.2)計(jì)算結(jié)果：N =10x = Columns 1 through 6 0 0.2000 0.4000 0.6000 0.8000 1.0000 Columns 7 through 11 1.2000 1.4000 1.6000 1.8000 2.0000y = Columns 1 through 5 1 5 25 125 625 Columns 6 through 10 3125 15625 78125 390625

18、 1953125 Column 11 97656252.4計(jì)算結(jié)果對(duì)比分析表2.1 步長(zhǎng)為0.05時(shí)的Runge-Kutta四階算法K001.00001010.05000.37500.3678790.0071220.10000.14060.1353350.0052630.15000.05270.0497870.0029140.20000.01980.0183160.0014850.25000.00740.0067380.0006660.30000.00280.0024790.0003270.35000.00100.0009128.8×10-580.40000.00040.000335

19、6.5×10-590.45000.00010.0001232.34×10-5100.50000.00014.54×10-55.5×10-5110.55000.00001.67×10-51.67×10-5表2.2 步長(zhǎng)為0.1時(shí)的Runge-Kutta四階算法K00 1.00001010.10000.33330.1353350.19796520.20000.11110.0183160.09278430.30000.03700.0024790.03452140.40000.01230.0003350.01196550.50000.0041

20、4.54×10-50.00405560.60000.00146.14×10-60.00139470.70000.00058.32×10-70.00049980.80000.00021.13×10-70.000290.90000.00011.52×10-81×10-4101.00000.00002.06×10-92.06×10-9表2.3 步長(zhǎng)為0.2時(shí)的Runge-Kutta四階算法K0011010.200050.0183164.98168420.4000250.00033524.9996630.60001256.

21、14×10-612540.80006251.13×10-762551.000031252.06×10-9312561.2000156253.78×10-111562571.4000781256.91×10-137812581.60003906251.27×10-1439062591.800019531252.32×10-161953125102.000097656254.25×10-189765625小結(jié)：（1）用標(biāo)準(zhǔn)的四階Runge-Kutta算法計(jì)算常微分方程時(shí)，不同的步長(zhǎng)對(duì)算法的穩(wěn)定性有影響。不夠穩(wěn)定的步長(zhǎng)下

22、面的計(jì)算，誤差會(huì)越來(lái)越大，結(jié)果失真嚴(yán)重。（2）一般情況下，步長(zhǎng)越小，標(biāo)準(zhǔn)的四階Runge-Kutta算法的穩(wěn)定性越高，精度也越高。3. 選做題題三利用牛頓法求方程于區(qū)間的根，考慮不同初值下牛頓法的收斂情況。3.1 選取初值為2計(jì)算結(jié)果： 8 3.097944090589694Elapsed time is 3.721679 seconds.3.2 選取初值為3計(jì)算結(jié)果： 6 2.982091433074349Elapsed time is 3.250174 seconds.3.3 選取初值為4計(jì)算結(jié)果： 7 4.348876621999810Elapsed time is 3.548942 s

23、econds.3.4 結(jié)果分析表3 不同初值下牛頓法的收斂情況234計(jì)算結(jié)果迭代次數(shù)867所用時(shí)間（秒）3.7216793.2501743.548942迭代過(guò)程n12.98209143307434924.348876621999810343.0979440905896945678由上表可知，當(dāng)初值離根越近，牛頓法收斂的速度越快，迭代次數(shù)越少，計(jì)算機(jī)的計(jì)算時(shí)間越少；相反，初值離根越遠(yuǎn)，收斂速度越慢，迭代次數(shù)越多，計(jì)算機(jī)的計(jì)算時(shí)間越多，占用的內(nèi)存也就越高。4. 選做題題四已知，利用復(fù)化梯形公式、復(fù)化Simpson公式和Romberg算法求的近似值；并觀察實(shí)際計(jì)算結(jié)果，比較它們的收斂速度。4.1 復(fù)

24、化梯形公式計(jì)算結(jié)果：4.2 復(fù)化Simpson公式計(jì)算結(jié)果：4.3 Romberg算法計(jì)算結(jié)果：4.4 結(jié)果分析的精確值為4根據(jù)以三種的計(jì)算結(jié)果可知，收斂速度最快的是Romberg算法，復(fù)化Simpson算法次之，復(fù)化梯形算法最慢，原因如下表所示。表4 Romberg算法、復(fù)化Simpson算法、復(fù)化梯形算法比較求積公式代數(shù)精度誤差量級(jí)Tn(f)1O(h2)Sn(f)3O(h4)Rn(f)7O(h8)其中上表中h為步長(zhǎng)，O(hm)表示同階無(wú)窮小。因此它們的收斂速度依次為Romberg算法>復(fù)化Simpson算法>復(fù)化梯形算法.總結(jié)通過(guò)本次的編程作業(yè)，我收獲很多。逐漸學(xué)會(huì)了MATLA

25、B編程語(yǔ)言的初步使用，在編程過(guò)程中加強(qiáng)了對(duì)高斯-塞德爾方法、初值問(wèn)題、數(shù)值積分收斂速度等方面的理解，更加培養(yǎng)了對(duì)于數(shù)值分析學(xué)習(xí)的樂(lè)趣，對(duì)我今后的學(xué)習(xí)工作都將會(huì)產(chǎn)生積極的影響。編程要重視最基本的原理，即這些方法的最基本的步驟，寫出合理的算法，能夠有效地提高編程速度。同一個(gè)題目，可能有程序方法實(shí)現(xiàn)，各種方法都有各自的有點(diǎn)，要合理選擇合適的方法。在大型計(jì)算中，計(jì)算速度非常重要，因此，我們要想辦法優(yōu)化程序，提高運(yùn)速度。對(duì)于實(shí)際的工程而言，比如土木、機(jī)械行業(yè)，在建立合適的力學(xué)模型后，絕大多數(shù)情況下，解析解難以求出，運(yùn)用編程的方法求出數(shù)值解，對(duì)于解決工程實(shí)際問(wèn)題，意義重大。因此，我們要重視課本上的基本知識(shí)

26、，基本迭代原理，加強(qiáng)編程訓(xùn)練，以利于以后的工作學(xué)習(xí)中實(shí)際問(wèn)題的解決。附錄程序一1.Gauss法調(diào)用M文件程序function x=Gauss(A,b)%線性方程組的高斯消去法%A為系數(shù)矩陣，b為方程組的右端項(xiàng)n,m=size(A);nb=length(b);if n=merror('系數(shù)矩陣必須是方的');endif m=nberror('b的維數(shù)與方程組的行數(shù)不匹配');endfor k=1:n-1a_max=0;for i=k:nif abs(A(i,k)>a_maxa_max=abs(A(i,k);r=i;endendif a_max<1e-1

27、5error('系數(shù)矩陣奇異');endif r>kfor j=k:n;z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%消元過(guò)程for i=k+1:nm=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endend%回帶過(guò)程if abs(A(n,n)<1e-15error('系數(shù)矩陣奇異，無(wú)法求解方程組');endx=zeros(size(b);for k=n:-1:1for j=k+1

28、:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end2. Gauss-seidelfunction x=Gauss_seidel(A,b)%線性方程組的Gauss-seidel迭代法%A為系數(shù)矩陣，b為方程組的右端項(xiàng)n,m=size(A);nb=length(b);if n=merror('系數(shù)矩陣必須是方的');endif m=nberror('b的維數(shù)與方程組的行數(shù)不匹配');endx=zeros(nb,1);k=1;B=-tril(A)triu(A,1);f=tril(A)b;while k<100y=B*x

29、+f;if norm(y-x,inf)<1e-5break;endx=y;k=k+1;end程序二Runge-Kutta調(diào)用M文件程序function x,y=runge_kutta(fun,a,b,y0,h)%四階runge_kutta算法%fun一階微分方程的函數(shù)%a,b為區(qū)間端點(diǎn)；y0為初始值%h為步長(zhǎng)，N為區(qū)間的等分?jǐn)?shù)N=(b-a)/hx=zeros(1,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=y0;for n=1:Nx(n+1)=x(n)+h;k1=h*feval(fun,x(n),y(n);k2=h*feval(fun,x(n)+1/2*h,y(n)

30、+1/2*k1);k3=h*feval(fun,x(n)+1/2*h,y(n)+1/2*k2);k4=h*feval(fun,x(n)+h,y(n)+k3);y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end程序三牛頓法調(diào)用M文件程序function x,iter,X=newton(fun,x0,eps,maxiter)%x非線性方程的近似根%iter迭代次數(shù)%X每步迭代的結(jié)果%fun輸入函數(shù)%x0初始迭代點(diǎn)%eps計(jì)算精度，默認(rèn)值為1e-6%maxiter最大迭代次數(shù)，默認(rèn)值為1e4tic %保存當(dāng)前時(shí)間if nargin<2,error('輸入?yún)?shù)至少2個(gè)'),endif nargin<3|isempty(eps),eps=1e-6;endif nargin<4|isempty(maxiter),maxiter=1e4;endk=0;err=1;while abs(err)>epsk=k+1;fx0,dfx0=feval(fun,x0);if dfx0=0error('fx在x0處的導(dǎo)數(shù)為0，終止運(yùn)算')endx1=x0-fx0/dfx0;err=x1-x0;x0=x1;X(k)=x1;endif k>=maxite