數(shù)值分析實驗教學指導(dǎo)書

1、數(shù)值分析實實 驗驗 指指 導(dǎo)導(dǎo) 書書徐州師范大學計算機科學與技術(shù)學院2009 年 11 月前 言數(shù)值分析實驗指導(dǎo)書根據(jù)數(shù)值分析實驗教學的需要所編寫，主要包括拉格朗日插值法、最小二乘法、數(shù)值積分、常微分方程初值問題數(shù)值解法、非線性方程求解、高斯消元法、數(shù)值計算方法應(yīng)用等七個方面的內(nèi)容。前六個實驗，都給出了實驗?zāi)康?、實驗要求、實驗原理和實驗?nèi)容。實驗原理主要介紹了算法的公式，并舉例說明該方法的使用，通過這一部分，可以有效的將上課內(nèi)容和實驗內(nèi)容對應(yīng)聯(lián)系起來。實驗內(nèi)容分為算法流程圖、必做題和選做題三種形式。算法流程圖在前兩個實驗給出，然后讓學生逐步自己畫流程圖，鍛煉了學生根據(jù)公式分析、書寫流程圖的能力

2、，加深學生對算法的理解；學生通過編寫必做題部分的程序、調(diào)試并運行，提高學生的編程能力；這兩部分的內(nèi)容是要求每一位學生都必須完成的。選做題是為基礎(chǔ)較好的同學準備的，提高學生的進一步解決問題的能力。本書不僅可以作為普通院校學生的數(shù)值分析上機實習指導(dǎo)教材，也可作為廣大讀者朋友掌握數(shù)值分析算法的實用教程。由于時間倉促以及水平有限，錯誤與不足之處難免存在，敬請使用者給與批評指正！目目 錄錄實驗一實驗一 拉格朗日插值法（拉格朗日插值法（2 課時）課時）.1一、實驗?zāi)康?1二、實驗要求.1三、實驗原理.1四、實驗內(nèi)容.2實驗二實驗二 最小二乘法（最小二乘法（4 課時）課時）.5一、實驗?zāi)康?5二、實驗要求.

3、5三、實驗原理.5四、實驗內(nèi)容.6實驗三實驗三 數(shù)值積分（數(shù)值積分（2 課時）課時）.10一、實驗?zāi)康?10二、實驗要求.10三、實驗原理.10四、實驗內(nèi)容.11實驗四實驗四 常微分方程初值問題數(shù)值解法（常微分方程初值問題數(shù)值解法（4 課時）課時）.15一、實驗?zāi)康?15二、實驗要求.15三、實驗原理.15四、實驗內(nèi)容.16實驗五實驗五 非線性方程求解（非線性方程求解（2 課時）課時）.19一、實驗?zāi)康?19二、實驗要求.19三、實驗原理.19四、實驗內(nèi)容.20實驗六實驗六 高斯消元法（高斯消元法（2 課時）課時）.23一、實驗?zāi)康?23二、實驗要求.23三、實驗原理.23四、實驗內(nèi)容.25實

4、驗七實驗七 數(shù)值計算方法應(yīng)用（數(shù)值計算方法應(yīng)用（6 課時）課時）.27一、實驗?zāi)康?27二、實驗要求.27三、實驗內(nèi)容.27參考文獻參考文獻.28實驗一實驗一 拉格朗日插值法（拉格朗日插值法（2 課時）課時）一、實驗?zāi)康囊?、實驗?zāi)康?.了解拉格朗日插值法的基本原理和方法。2.掌握拉格朗日插值多項式的用法。二、實驗要求二、實驗要求1.用 C 語言編程實現(xiàn)拉格朗日插值。2.進一步加深對拉格朗日插值法的理解。三、實驗原理三、實驗原理（一）拉格朗日插值公式（一）拉格朗日插值公式經(jīng)過個點，構(gòu)造一個 n 次多項式，形如：) 1( n),( ,),(),(1100nnyxyxyxnkkknxlyxp0)()

5、(使得 成立。kknyxp)(),2 , 1 , 0(nk其中kknkkkkkknkkkxxxxxxxxxxxxxxxxxxxxxl10)()()()()()()(110110為插值基函數(shù)。（二）例子（二）例子已知滿足作的二次拉格朗日插值多)(xf15)225(,13)169(,12)144(fff)(xf項式，并求的近似值。)175(f解：設(shè)，則的二次拉格朗15,13,12225,169,144210210yyyxxx，)(xf日插值基函數(shù)為：4536)169)(144()()()(1400)225)(144()()()(2025)225)(169()()()(1202102 2101201

6、 2010210 xxxxxxxxxxxlxxxxxxxxxxxlxxxxxxxxxxxl因此的二次拉格朗日插值多項式為；)(xf)()()()(2211002xlyxlyxlyxL且。73158230.13)175(15)175(13)175(12)175()175( 210 2lllLf四、實驗內(nèi)容四、實驗內(nèi)容（一）算法流程圖（一）算法流程圖如果程序的數(shù)據(jù)輸入項（函數(shù)參數(shù)）為：插值節(jié)點及函數(shù)值，及待求點 x 的值；輸出為待求點 x 對應(yīng)的函數(shù)值，則程序流程圖如圖 1-1。輸入 xi,yi,及 n,x,k=0y = 0P =1j=0,1nk=j?P = P*(x-xj)/(xk-xj)y =

7、 y+P*ykk=k+1輸出x,ykn圖 1-1 拉格朗日插值法算法流程圖（二）編程作業(yè)（二）編程作業(yè)編寫拉格朗日插值法通用子程序，并用以下函數(shù)表來上機求,。)15. 0(f)31. 0(fx0.00.10.1950.30.4010.5f (x)0.398940.396950.391420.381380.368120.35206代碼如下： / * C 語言源代碼 */#include main() static float Lx10,Ly10;int n,i,j;float x,y,p;printf(enter n=);scanf(%d,&n); / printf(enter xin)

8、; for(i=1;i=n;i+) scanf(%f,&Lxi);printf(enter yin); for(i=1;i=n;i+) scanf(%f,&Lyi);printf(enter x=);scanf(%f,&x); for( ) ; for( ) if(i!=j) ; / ； printf(y=%fn,y);請完成這個程序，并在兩處注釋處寫上正確的注釋。在執(zhí)行程序時，如果求的值，那么屏幕上應(yīng)出現(xiàn)如下內(nèi)容：)15. 0(fenter n=6enter xi00.10.1950.30.4010.5enter yi0.398940.396950.391420.38

9、1380.368120.35206enter x=0.15y=0.394473（三）選做題（三）選做題參考教材牛頓插值公式，編程實現(xiàn)用牛頓插值公式求上述條件下對應(yīng)節(jié)點的函數(shù)值。實驗二實驗二 最小二乘法（最小二乘法（4 課時）課時）一、實驗?zāi)康囊弧嶒災(zāi)康?.了解最小二乘擬合的基本原理和方法，注意與插值方法的區(qū)別。2.掌握最小二乘法。二、實驗要求二、實驗要求1.掌握用 C 語言作最小二乘多項式擬合的方法。2.進一步加深對最小二乘法的理解。三、實驗原理三、實驗原理（一）最小二乘多項式擬合（一）最小二乘多項式擬合已知數(shù)據(jù)對，求多項式), 2 , 1)(,njyxjj)()(0nmxaxPmiii使得

10、為最小，這就是一個最小二乘問題。 njmijijinyxaaaa12010),(（二）法方程組為（二）法方程組為 . a.aa . 111n10121111112111miinimiiimiiminiminiminiminimiimiiminimiiyxyxyxxxxxxxxm（三）最小二乘法計算步驟（三）最小二乘法計算步驟用線性函數(shù)為例，擬合給定數(shù)據(jù)。bxaxP)(), 2 , 1)(,miyxii算法描述：步驟 1：輸入值，及。m), 2 , 1)(,miyxii步驟 2：建立法方程組。AYAXAT步驟 3：解法方程組。步驟 4：輸出。bxaxP)(四、實驗內(nèi)容四、實驗內(nèi)容（一）算法流程圖

11、（一）算法流程圖1.算法整體流程圖（如圖 2-1）輸入xi,yi(i=1,2,m)及m,n生成中間矩陣C生成法方程組的系數(shù)矩陣A=CTC生成法方程組的右端向量b=CTY解法方程中得ai(i=0,1,n)輸出ai(i=0,1,n)圖 2-1 最小二乘法算法整體流程圖2.“生成中間矩陣 C”算法流程圖（如圖 2-2）i=1,2,mCi11j=2,3,n+1Cijxi*Cij-1圖 2-2 最小二乘法中“生成中間矩陣 C”算法流程圖3.中間矩陣 C 的重要作用miinimiiimiiimiiminiminiminiminiminimiimiimiiminimiimiinnnmnmnnyxyxyxyY

12、xxxxxxxxxxxmAxxxxxxxxxYxxxxxxxxxmmm11211T121211111131211121T22222Tm212222Cb,.CC,.1.11Cyyy ,.1.1.1C1221122111則（二）編程作業(yè)（二）編程作業(yè)測得銅導(dǎo)線在溫度()時的電阻如下表，求電阻 R 與溫度 T 的近似函數(shù)關(guān)Ti)(Ri系。i0123456()Ti19.125.030.136.040.045.150.0)(Ri76.3077.8079.2580.8082.3583.9085.10【提示提示】在進行程序?qū)崿F(xiàn)時，務(wù)必注意中間矩陣的作用，以及非奇次線性方程組求解問題！為了實驗的順利完成，此處

13、給出解非奇次線性方程組的高斯消元法的函數(shù)。請認真閱讀并理解。float gs(float a2020,float b20,int n )int i,j,k,l; float s; k=1; while(k!=n+1) if(akk!=0) for(i=k+1;i=n+1;i+) aik=aik/akk; bi=bi-aik*bk; for(j=k+1;j=1;k-) s=0; for(l=k+1;l=n+1;l+) s=s+akl*bl; bk=(bk-s)/akk; return 0;實驗主程序如下（請加上必要的注釋） 。實驗三實驗三 數(shù)值積分（數(shù)值積分（2 課時）課時）一、實驗?zāi)康囊?、實?/p>

14、目的1.了解數(shù)值積分的基本原理和方法。2.掌握復(fù)合梯形公式。3.了解求積公式外推思想、Romberg 公式及 Romberg 積分法。二、實驗要求二、實驗要求1.編寫定步長復(fù)合梯形公式。2.編寫變步長復(fù)合梯形公式。3.進一步加深對數(shù)值積分的理解。三、實驗原理三、實驗原理（一）定步長復(fù)合梯形公式（一）定步長復(fù)合梯形公式1.公式將積分區(qū)間等分，分點為，其中稱為積分nba,niihaxi, 1 , 0，nabh步長。)()(2)(2)(11bfxfafhTdxxfnkknba2.例子用復(fù)合梯形公式求積分的近似值。 （取 8 位小數(shù)，精確解為dxx10214）14159265. 3138988. 3)

15、1 ()87()43()85()21()83()41()81( 2)0(1618fffffffffT（二）變步長復(fù)合梯形公式（二）變步長復(fù)合梯形公式1.公式nknnnabkafnabTT122) 12(221遞推公式：), 3 , 2 , 1(2) 12(221)()(21121221kabiafabTTbfafabTkkkikk2.例子用遞推公式求積分的近似值，使誤差不超過。dxx1021461013117647. 3)43()41(41211 . 3)21(212132421) 1 ()0(2124121ffTTfTTffT四、實驗內(nèi)容四、實驗內(nèi)容（一）算法流程圖（一）算法流程圖1.定步長

16、復(fù)合梯形算法流程圖2.變步長復(fù)合梯形算法流程圖（二）編程作業(yè)（二）編程作業(yè)求 的近似值。10214x（1）編寫定步長復(fù)合梯形程序求解上式；（2）編寫變步長復(fù)合梯形程序求解上式，使誤差不超過 10-6?！咎崾尽空埜鶕?jù)前面的算法流程圖進行編寫程序。1.定步長復(fù)合梯形程序求解2.變步長復(fù)合梯形程序求解實驗四實驗四 常微分方程初值問題數(shù)值解法（常微分方程初值問題數(shù)值解法（4 課時）課時）一、實驗?zāi)康囊弧嶒災(zāi)康?.掌握 Euler 法和改進的 Euler 法公式的用法。2.進一步加深對微分方程數(shù)值解的理解。二、實驗要求二、實驗要求1.編寫歐拉法程序。2.編寫改進的歐拉法程序，學會用改進的歐拉公式來求解

17、常微分方程初值問題。三、實驗原理三、實驗原理（一）常微分方程初值問題（一）常微分方程初值問題 (式 1)的數(shù)值解法，這也是科學與工程計算經(jīng)常遇到的問題。（二）歐拉法（二）歐拉法求初值問題(式 1)的一種最簡單方法是將節(jié)點的導(dǎo)數(shù)用差商代替，于是(式 1)的方程可近似寫成 (式 2)從出發(fā)，由(式 2)求得再將代入(式 2)右端，得到的近似，一般寫成稱為解初值問題的 Euler 法。（三）改進歐拉法（三）改進歐拉法先用 Euler 法計算出的近似，將隱式梯形公式改為稱為改進 Euler 法，它實際上是顯式方法。即（四）例子（四）例子求初值問題的數(shù)值解，取步長。 （精確解為1)0(102yxyxyy

18、，1 . 0h）2/1)21 ()(xxy解：（1）利用歐拉法9 ,2 , 1 , 0,1/2 . 01 . 101iyyxyyiiii（2）利用改進歐拉法9 ,2 , 1 , 0,11 . 0) 1 . 0(21 . 0/2)(05. 001121211iyKyxKyKyxyKKKyyiiiiiiii計算結(jié)果如下：ixiEuler 方法 yi改進 Euler 法 yi精確解 y(xi)01234567891000.10.20.30.40.50.60.70.80.9111.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.71

19、77791.78477011.0959091.1840961.2662011.3433601.4164021.4859561.5525151.6164761.6781681.73786911.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051四、實驗內(nèi)容四、實驗內(nèi)容（一）算法流程圖（一）算法流程圖1.請根據(jù)歐拉公式，畫出其算法流程圖。2.請根據(jù)改進歐拉公式，畫出其算法流程圖。（二）編程作業(yè)（二）編程作業(yè)編寫 Euler 法和改進的 Euler 法程序。求微分方程在區(qū)間 10 等1y(0)0

20、,1xy,dxdy份的近似解。1.Euler 法程序2.改進的 Euler 法程序在編寫改進的 Euler 法程序時，有關(guān)輸入輸出部分，可參照以下屏幕上應(yīng)出現(xiàn)內(nèi)容。please input a,b,h and a00 1 0.1 1x=0.000000, y=1.000000 x=0.100000, y=0.905000 x=0.200000, y=0.819025x=0.300000, y=0.741218x=0.400000, y=0.670802x=0.500000, y=0.607076x=0.600000, y=0.549404x=0.700000, y=0.497210 x=0.8

21、00000, y=0.449975x=0.900000, y=0.407228x=1.000000, y=0.368541（三）選做題（三）選做題使用梯形公式編寫程序，解決上述例子問題。 y x第第1章章 O x* x1 x0實驗五實驗五 非線性方程求解（非線性方程求解（2 課時）課時）一、實驗?zāi)康囊?、實驗?zāi)康?.了解求解非線性方程的解的常見方法。2.編寫牛頓迭代法程序求解非線性方程。二、實驗要求二、實驗要求1.設(shè)計牛頓迭代法算法，編寫程序上機調(diào)試。2.進一步加深對迭代法的理解。三、實驗原理三、實驗原理（一）牛頓迭代法（一）牛頓迭代法又稱為牛頓-雷夫生方法（Newton-Raphson met

22、hod） ，是一種在實數(shù)域和復(fù)數(shù)域上通過迭代計算求出非線性方程的數(shù)值解方法。方法的基本思路是利用一個根的猜測值x0做初始近似值，使用函數(shù) f(x)在 x0處的泰勒級數(shù)展式的前兩項做為函數(shù) f(x)的近似表達式。由于該表達式是一個線性函數(shù)，通過線性表達式替代方程 f(x)= 0 中的 f(x)求得近似解 x1。即將方程 f(x)= 0 在 x0處局部線性化計算出近似解 x1，重復(fù)這一過程，將方程 f(x)= 0 在 x1處局部線性化計算出 x2，求得近似解 x2，。詳細敘述如下：假設(shè)方程的解x*在 x0附近（x0是方程解 x*的近似） ，函數(shù) f(x)在點 x0處的局部線化表達式為)()()()

23、(000 xfxxxfxf由此得一次方程0)()()(000 xfxxxf求解，得)()(0001xfxfxx如圖 5-1 所示，x1比 x0更接近于 x*。該方法的幾何意義是：用曲線上某點（x0，y0）的切線代替曲線，以該切線與 x 軸的交點（x1，0）作為曲線與 x 軸的交點（x*，0）的近似（所以牛頓迭代法又稱為切線法） 。設(shè) xn是方程解 x*的近似，迭代格式 ( n = 0，1，2，)()(1nnnnxfxfxx就是著名的牛頓迭代公式，通過迭代計算實現(xiàn)逐次逼近方程的解。牛頓迭代法的圖 5-1 牛頓迭代法示意圖最大優(yōu)點是收斂速度快，具有二階收斂。以著名的平方根算法為例，說明二階收斂速度

24、的意義。（二）例子（二）例子已知，求等價于求方程 f(x) = x2 2 = 0 的解。由于。應(yīng)用牛4 . 12 2xxf2)(頓迭代法，得迭代計算格式， （n = 0，1，2，）)/2(211nnnxxx取 x0= 1.4 為初值，迭代計算 3 次的數(shù)據(jù)列表如表 5-1：表 5-1 牛頓迭代法數(shù)值實驗迭代次數(shù)近似值15 位有效數(shù)誤差01.41.41421356237310-1.42e-00211.414285714285711.414213562373107.21e-00521.414213564213561.414213562373101.84e-00931.414213562373091

25、.41421356237310-2.22e-016四、實驗內(nèi)容四、實驗內(nèi)容（一）算法流程圖（一）算法流程圖牛頓迭代法算法流程圖（二）編程作業(yè)（二）編程作業(yè)編寫 Newton 迭代法通用子程序。實現(xiàn)方程 f(x)=x6-x-1=0 的滿足精度要求的解。要求求解過程中用一個變量 I 控制三種狀態(tài)，其中：i=0 表示求解滿足給定精度的近似解；i=1 表示 f(x0)=0，計算中斷；i=2 表示迭代 n 次后精度要求仍不滿足。（說明：用教材 P139 例 8）（三）選做題（三）選做題1.用簡化牛頓法計算上述例題。2.用牛頓下山法計算上述例題。實驗六實驗六 高斯消元法（高斯消元法（2 課時）課時）一、實

26、驗?zāi)康囊?、實驗?zāi)康?.掌握高斯選主元消去法公式的用法，適用范圍及精確度。2.通過高斯選主元消去法求矩陣方程的解，驗證高斯消去法。二、實驗要求二、實驗要求1.寫出高斯選主元消去法解線性方程組算法，編寫程序上機調(diào)試出結(jié)果。2.進一步加深對高斯消去法的理解。三、實驗原理三、實驗原理（一）高斯消元法（一）高斯消元法Gauss 消去法就是將方程組nnnnnnnnnnbxaxaxabxaxaxabxaxaxa.22112222212111212111通過(n-1)步消元，將上述方程組轉(zhuǎn)化為上三角方程組再回代求此方程組的解。（二）選列主元素高斯消去法（二）選列主元素高斯消去法給定線性方程組 Ax=b，記 A

27、(1)=A，b(1)=b，列主元 Gauss 消去法的具體過程如下：首先在增廣矩陣 B(1)=(A(1),b(1)的第一列元素中，取.max1)1(11)1(1rraakinik為主元素，然后進行第一步消元得增廣矩陣 B(2)=(A(2),b(2)。 再在矩陣 B(2)=(A(2),b(2)的第二列元素中，取.max2)2(22)2(2rraakinik為主元素，然后進行第二步消元得增廣矩陣 B(3)=(A(3),b(3)。按此方法繼續(xù)進行下去, 經(jīng)過n-1 步選主元和消元運算,得到增廣矩陣 B(n)=(A(n),b(n).則方程組 A(n)x=b(n)是與原方程組等價的上三角形方程組,可進行

28、回代求解. 易證，只要|A|0,列主元 Gauss 消去法就可順利進行。（三）例子（三）例子采用 4 位十進制浮點計算，分別用順序 Gauss 消去法和列主元 Gauss 消去法求解線性方程組：9812 . 4120032001 .1291. 58334. 06781. 0167. 001. 0012. 0321321321xxxxxxxxx方程組具有四位有效數(shù)字的精確解為 x1*=17.46，x2*=-45.76，x3*=5.546解（1）用順序 Gauss 消去法求解，消元過程為553231065171011750041.44010. 8101000. 006781. 01670. 001

29、00. 00120. 01017981044531467041.44010. 8101000. 006781. 01670. 00100. 00120. 00 .981200. 41200320010.12910. 58334. 0000. 16781. 0167. 00100. 00120. 0回代得： x3=5.546，x2=100.0，x1=-104.0 （2）用列主元 Gauss 消去法求解，消元過程為選主元31rr 6744. 01670. 01055. 0079.11909. 54584. 000 .981200. 4120032006781. 01670. 00100. 00120. 010.12910. 58334. 0000. 10 .981200. 4120032000 .981200. 41200320010.12910. 58334. 0000. 16781. 01670. 00100. 00120. 025329. 00961. 00079.11909. 54584. 000 .981200. 4120032005329. 00961. 00079