數(shù)值分析課程上機作業(yè)計算報告

1、西南交通大學-數(shù)值分析-報告PAGE 18/21數(shù)值分析課程上機作業(yè)計算報告班 級： 學 號： 姓 名： 專 業(yè)：大地測量學及測繪工程 指導老師： 聯(lián)系電話： 序 言 通過數(shù)值分析的理論知識的學習，此次實驗將我們學過的理論知識運用于實踐之中。本次實驗，我選用的計算機語言為MATLAB，其主要有一下幾個特點。1編程效率高MATLAB是一種面向科學與工程計算的高級語言，允許使用數(shù)學形式的語言編寫程序，且比BASIC、FORTRAN和C等語言更加接近我們書寫計算公式的思維方式，用MATLAB編寫程序猶如在演算紙上排列出公式與求解問題。因此，MATLAB語言也可通俗地稱為演算紙式科學算法語言。由于它編

2、寫簡單，所以編程效率高，易學易懂2. 用戶使用方便MATLAB語言與其他語言相比，較好的解決了上述問題，把編輯、編譯、鏈接和執(zhí)行融為一體。它能在同一畫面上進行靈活操作，快速排除輸入程序中的書寫錯誤、語法錯誤以至語義錯誤，從而加快了用戶編寫、修改和調試程序的速度，可以說在編程和調試過程中它是一種比VB還要簡單的語言。3. 方便的繪圖功能MATLAB的繪圖是十分方便的，它有一系列繪圖函數(shù)（命令），例如線性坐標、對數(shù)坐標、半對數(shù)坐標及極坐標，均只需調用不同的繪圖函數(shù)（命令），在圖上標出圖題、XY軸標注，格（柵）繪制也只需調用相應的命令，簡單易行。另外，在調用繪圖函數(shù)時調整自變量可繪出不變顏色的點、線

3、、復線或多重線。這種為科學研究著想的設計是通用的編程語言所不能及的。目 錄 TOC o 1-3 h z u HYPERLINK l _Toc311462221 1.實驗一1 HYPERLINK l _Toc311462222 1.1題目1 HYPERLINK l _Toc311462223 1.2計算思路1 HYPERLINK l _Toc311462227 1.3計算結果1 HYPERLINK l _Toc311462228 1.4總結6 HYPERLINK l _Toc311462229 2.第二題7 HYPERLINK l _Toc311462230 2.1題目7 HYPERLINK l

4、 _Toc311462231 2.2 松弛思想分析7 HYPERLINK l _Toc311462232 2.3問題的求解7 HYPERLINK l _Toc311462233 2.4總結10 HYPERLINK l _Toc311462234 3.第三題11 HYPERLINK l _Toc311462235 3.1題目11 HYPERLINK l _Toc311462236 3.2 Runge-Kutta法的基本思想11 HYPERLINK l _Toc311462237 3.3 問題的求解 PAGEREF _Toc311462237 h 11 HYPERLINK l _Toc311462

5、238 3.4問題的總結 PAGEREF _Toc311462238 h 14 HYPERLINK l _Toc311462239 總結15 HYPERLINK l _Toc311462240 附件16 HYPERLINK l _Toc311462241 實驗一程序設計16 HYPERLINK l _Toc311462242 實驗二程序設計16 HYPERLINK l _Toc311462243 實驗三程序設計171實驗一：插值問題1.1題目已知：a=-5,b=5, 以下是某函數(shù) f(x)的一些點（xk,yk）, 其中xk=a+0.1(k-1) ,k=1,.,101；（數(shù)據(jù)略）。請用插值類方法

6、給出函數(shù)f(x)的一個解決方案和具體結果。并通過實驗考慮下列問題：（1）Ln(x)的次數(shù)n越高，逼近f(x)的程度越好？（2）高次插值收斂性如何？（3）如何選擇等距插值多項式次數(shù) ？（4）若要精度增高，你有什么想法？ 比如一定用插值嗎？（5）逼近某個函數(shù)不用插值方式，有何變通之舉？（6）函數(shù)之間的誤差如何度量，逼近的標準又是什么？(7) 如何比較好的使用插值多項式呢？1.2計算思路本題我選用拉格朗日插值函數(shù)，拉格朗日插值函數(shù)的構建算法是從101組數(shù)據(jù)中等距選取n組數(shù)據(jù)來構造n階的拉格朗日插值函數(shù)，然后在-5,5區(qū)間選取m個插值點帶入插值多項式計算出m組（x,y），本題中，我分別取了n=6,8,

7、10，20四種不同的階數(shù)，然后利用matlab繪圖命令查看插值函數(shù)圖像與原函數(shù)圖像的逼近效果。1.3 計算結果當階數(shù)為5階時，插值圖像與原函數(shù)圖象如下所示（紅色曲線為插值函數(shù)曲線，藍色為原函數(shù)曲線）表1（迭代數(shù)據(jù)）iXiYi1-5.0000 25.000026-2.50006.2693510 10.0000711.0000 6.26931015.000025.0000 由圖象可以看出，除了在插值點附近的擬合效果還可以，其他的點都差距比較大，總的來說效果不是很好。當階數(shù)為7階時圖象如下：表2（迭代數(shù)據(jù)）iXiYi1-5.0000 25.000017-3.400011.560033-1.8000

8、2.883249-0.2000 5.2312651.40003.0219 813.00008.9991 974.600021.1600由此圖象可以看出，利用7階插值來進行計算，在0值附近的擬合不是很好，但總得來說比5階函數(shù)較好。當階數(shù)為10階時，圖象如下：表3（迭代數(shù)據(jù)）iXiYi1-5.0000 25.000012-3.900015.210023-2.8000 7.840534-1.7000 2.5554. 893.800014.4400 1004.900024.0100利用10階函數(shù)進行插值計算時，在函數(shù)的兩端會出現(xiàn)相比于原函數(shù)偏離比較大的點，這就是runge（龍格）現(xiàn)象，但不是很明顯。當

9、階數(shù)為20階時，圖象如下：表4（迭代數(shù)據(jù)）iXiYi1-5.0000 25.00006-4.500020.250011-4.0000 16.000016-3.5000 12.2500. 914.000016.0000 964.500020.2500由此圖像可以看出，runge（龍格）現(xiàn)象已經(jīng)非常明顯了，由于在兩端的龍格現(xiàn)象太明顯，所以此圖是在去除了兩端偏離比較大的點。1.4 總結利用拉格朗日插值函數(shù)構造插值多項式是比較常用的插值計算的方法之一，其解決了不需要求解方程組的問題，它也有一定的缺點：（1）從實驗中可以看出，并不是拉格朗日的插值階數(shù)越高，其插值效果越好，雖然在一些局部擬合的效果較好，但

10、是當階數(shù)超過7階時就容易出現(xiàn)runge現(xiàn)象，由上圖中的階數(shù)n=20時就可以看出；（2）但是當插值階數(shù)太少時，插值的效果又不是很好，由上圖中階數(shù)n=5時可以看出，插值效果非常不好；（3）由上面的圖象可以分析，利用拉格朗日進行插值計算時，所選取的插值點并不是在區(qū)間里平均選取最好，根據(jù)不同的函數(shù)曲線，應當用不同的方法選點。實驗二2.1 題目松弛因子對SOR法收斂速度的影響。用SOR法求解方程組Ax=b,其中要求程序中不存系數(shù)矩陣A,分別對不同的階數(shù)取w=1.1, 1.2, .,1.9進行迭代，記錄近似解x(k)達到|x(k)-x(k-1)|10-6時所用的迭代次數(shù)k，觀察松弛因子對收斂速度的影響，并

11、觀察當w0或w2會有什么影響？2.2 松弛思想分析SOR法（松弛法），是由高斯-塞德爾迭代方法演變來的。高斯-塞德爾迭代方法的迭代格式為，而松弛法的迭代格式為，通過大量的實驗發(fā)現(xiàn)，的取值很關鍵，如果取了好的，迭代方程的收斂速度會加快，根據(jù)理論基礎可以知道，當02時松弛法才能夠收斂。2.3問題的求解分別取不同階數(shù)系數(shù)矩陣和不同松弛因子進行計算，結果如下圖所示：當取階數(shù)為5時，=1. 1時：迭代次數(shù)：迭代過程：x = 0.8250 0.7769 0.7636 0.7600 1.0340 x = 0.9561 0.9453 0.9426 1.0176 1.0014x = 0.9893 0.9868

12、1.0069 1.0005 1.0000 x = 0.9974 1.0025 1.0001 1.0000 1.0000 x = 1.0010 1.0001 1.0000 1.0000 1.0000 x = 0.9999 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000當階數(shù)n=6，=1. 3時：迭代次數(shù)：迭代過程：x = 0.9750 0.9669 0.9642 0.9634 0.9631 1.2880 x = 0.9967 0.99

13、73 0.9979 0.9983 1.1041 0.9474x = 1.0001 1.0002 1.0001 1.0344 0.9629 1.0037x = 1.0000 1.0000 1.0111 0.9812 1.0062 1.0009x = 1.0000 1.0036 0.9917 1.0050 1.0000 0.9997x = 1.0012 0.9966 1.0030 0.9995 0.9997 1.0000 x = 0.9985 1.0015 0.9994 0.9999 1.0000 1.0000 x = 1.0009 0.9997 1.0000 1.0001 1.0000 1.00

14、00 x = 0.9996 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0001 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000當階數(shù)n=5時，=-1時：請輸入階數(shù)n：5請

15、輸入w的值：-1num= 1, error=2.53e+000num= 2, error=3.93e+000num= 3, error=7.35e+000num= 4, error=1.80e+001num= 5, error=4.84e+001num= 6, error=1.32e+002num= 7, error=3.63e+002num= 10; error=1.00e+003num= 11, error=2.77e+003num= 12, error=7.68e+003num= 13, error=2.12e+004num= 14, error=5.87e+004num= 15, er

16、ror=1.62e+005num= 16, error=4.47e+005num= 17, error=1.23e+006num= 20, error=3.38e+006num= 21, error=9.27e+006num= 22, error=2.54e+007num= 23, error=6.96e+007num= 24, error=1.90e+008num= 25, error=5.20e+008num= 26, error=1.42e+009num= 27, error=3.87e+009num= 30, error=1.06e+010num= 31, error=2.88e+01

17、0num= 32, error=7.84e+010num= 33, error=2.13e+011num= 34, error=5.81e+011num= 35, error=1.58e+012num= 36, error=4.29e+012x = 1.0e+012 * -0.3019 0.8370 -1.5384 2.1285 -1.9856由數(shù)據(jù)中可以分析得出，當=-1時,迭代是不收斂的。當階數(shù)n=5時，=2.5時請輸入階數(shù)n：5請輸入w的值：2.5num= 1, error=1.38e+001num= 2, error=2.11e+001num= 3, error=2.58e+001nu

18、m= 4, error=2.80e+001num= 5, error=4.63e+001num= 6, error=6.61e+001num= 7, error=8.74e+001num= 10, error=1.38e+002num= 11, error=1.72e+002num= 12, error=2.63e+002num= 13, error=4.26e+002num= 14, error=8.10e+002num= 15, error=1.51e+003num= 16, error=2.52e+003num= 17, error=3.58e+003num= 20, error=4.4

19、6e+003num= 21, error=6.34e+003num= 22, error=9.29e+003num= 23, error=1.25e+004num= 24, error=1.83e+004num= 25, error=2.82e+004num= 26, error=3.74e+004num= 27, error=6.20e+004num= 30, error=1.18e+005num= 31, error=1.96e+005num= 32, error=2.94e+005num= 33, error=4.09e+005num= 34, error=5.63e+005num= 3

20、5, error=7.87e+005num= 36, error=1.24e+006x = 1.0e+005 * 0.4613 1.5993 -1.4788 1.3274 -2.6027由數(shù)據(jù)中可以分析得出，當=2.5時,迭代也是不收斂的。2.4 總結（1）在本題中，迭代次數(shù)和矩陣的未知數(shù)有以及的取值有一定的聯(lián)系，不一定是未知數(shù)多迭代次數(shù)就多，也有可能在取相同時，未知數(shù)多迭代次數(shù)反而較少，這和SOR法的迭代原理有關。（2）當0或2.0，SOR法不收斂，這是必要條件。實驗三3.1 題目用Runge-Kutta 4階算法對初值問題y/=-20*y，y(0)=1按不同步長求解，用于觀察穩(wěn)定區(qū)間的作用

21、，推薦兩種步長h=0.1,0.2。注：此方程的精確解為：y=e-20 x3.2 Runge-Kutta法的基本思想Runge-Kutta法的基本思想是通過f (x, y)某些點函數(shù)值的適當線性組合替換Euler 法中的f(xk,yk)，可能使得方法的精確度更高。由于有了這種思想，因此標準的四階Runge-Kutta的算法如下(h為求解步長）：K1=hf(xk,yk)如果yk某種方法第k步的近似值，y(xk)是其準確值，其絕對誤差為k，即有k=y(xk)-yk。假定第k步之后的計算不再有舍入誤差，只是由k引起的擾動m，都有mk，則稱此方法是絕對穩(wěn)定的。數(shù)值解法的穩(wěn)定性一般與步長h以及f(x,y)

22、有關。3.3 問題的求解給定初值x為0，我選定迭代終點為1，,分別選取h=0.1、h=0.2進行計算，并對迭代過程和迭代圖像進行分析。計算結果如下：當h=0.1時請輸入a的值：0請輸入b的值：1請輸入步長h：0.1x 1=0.000000 y 1=1.000000 yi 1= 1.000000 error 1= 0.000000 x 2=0.100000 y 2=0.333333 yi 2= 0.135335 error 2= 0.197998x 3=0.200000 y 3=0.111111 yi 3= 0.018316 error 3= 0.092795x 4=0.300000 y 4=0

23、.037037 yi 4= 0.002479 error 4= 0.034558x 5=0.400000 y 5=0.012346 yi 5= 0.000335 error 5= 0.012010 x 6=0.500000 y 6=0.004115 yi 6= 0.000045 error 6= 0.004070 x 7=0.600000 y 7=0.001372 yi 7= 0.000006 error 7= 0.001366x 8=0.700000 y 8=0.000457 yi 8= 0.000001 error 8= 0.000456x 9=0.800000 y 9=0.000152

24、yi 9= 0.000000 error 9= 0.000152當h=0.2時由圖象和數(shù)據(jù)中可以看出，當步長為0.2時，迭代是不收斂的。3.4問題的總結（1）用標準的四階 Runge-Kutta 算法計算常微分方程時，不同的步長算法的穩(wěn)定性有影響。不夠穩(wěn)定的步長下面的計算，誤差會越來越大，結果失真嚴重。（2）一般情況下，步長越小，標準的四階 Runge-Kutta 算法的穩(wěn)定性越高，精度也越高?？偨Y通過此次數(shù)值分析上機實習，認識到了計算機編程在數(shù)值分析課程中的重要性。也讓我意識到了學習編程語言的重要性。要想要將數(shù)學應用于實際工程中，特別是對于工科的學生來講，這是一門非常重要的實踐課程。在此次報

25、告中，首次接觸了Matlab這門軟件，之所以選這門，是因為很多工程中對數(shù)據(jù)的處理需要使用，這對我本身的專業(yè)是非常重要的。本學期學習的數(shù)值分析方法，讓我在對數(shù)學的理解上又有了新的認識。在一連串看似雜亂無章的數(shù)據(jù)中，用數(shù)值分析進行處理，有了很多的收獲。在學習的過程中，意識到了這門課對我專業(yè)的重要性，讓我對自己的專業(yè)有了更好的發(fā)揮。對我寫論文以及解決實際工程問題有很大的幫助。最后，感謝老師給我這樣的機會去接觸這門語言，雖然只了解了皮毛，可是仍然收獲頗多。由于初次接觸這門軟件，在報告中仍難免會有不完善甚至錯誤的地方，望諒解！附錄實驗一程序function yy=L(xi) f=fopen(xk.txt,r); x=fscanf(f,%f); g=fopen(yk.txt,r); y=fscanf(g,%f); n=input(請輸入插值函數(shù)階數(shù)：); s=0; xi=-3.5:0.1:3.5; for i=1:fix(101/(n-1):101 t=ones(1,length(xi); for j=1:fix(101/(n-1):101 if j=i t=t.*(xi-x(j)/(x(i)-x(j); end end s=s+t*y(i); end yy=s; plot(x,y,b,xi,