工程數(shù)值方法

1、工程數(shù)值方法學(xué)習(xí)內(nèi)容：Chapter 1 線性代數(shù)方程組的數(shù)值解法Chapter 2 插值問題與數(shù)值微分Chapter 3 數(shù)值積分方法Chapter 4 常微分方程(組)初值問題的數(shù)值方法Chapter 5 常微分方程(組)邊值問題的數(shù)值方法Chapter 6 橢圓型偏微分方程的數(shù)值方法Chapter 7 加權(quán)殘值方法參考書目：1 武漢大學(xué)、山東大學(xué)合編，計(jì)算方法，高教版，19792 林成森編，數(shù)值計(jì)算方法（上、下），科學(xué)出版社，20003 中科院研究生數(shù)學(xué)叢書，工程中的數(shù)值方法，科學(xué)出版社，20004 曾紹林編，工程數(shù)學(xué)基礎(chǔ)（研究生數(shù)學(xué)叢書），科學(xué)出版社，20015 李慶揚(yáng)編，數(shù)值分析基礎(chǔ)

2、教程，高等教育出版社，20026 李慶揚(yáng)編，數(shù)值分析（第4版），清華版，20037 關(guān)治編，數(shù)值計(jì)算方法，清華版，20048 李岳生、黃有謙編，數(shù)值逼近，高教版，19789 李榮華編，微分方程數(shù)值解法，人教版，198010 邱吉寶編著，加權(quán)殘值法的理論與應(yīng)用，宇航版，1992Chapter 1 線性代數(shù)方程組的數(shù)值解法 線性代數(shù)方程組的求解是工程實(shí)踐中最常遇到的問題。據(jù)不完全統(tǒng)計(jì)，在工程實(shí)踐中提出的計(jì)算問題中，有近一半涉及到求解線性方程組。例如：結(jié)構(gòu)有限元分析問題，大地測量問題，氣象預(yù)報(bào)問題，電力傳輸網(wǎng)分析問題，各種電路分析問題，數(shù)據(jù)擬合問題，以及非線性方程組與微分方程的數(shù)值求解問題等等。因此

3、，學(xué)習(xí)并掌握線性代數(shù)方程組求解的基本理論與方法無疑是十分必需的。 本章將介紹目前一些利用計(jì)算機(jī)求解線性代數(shù)方程組常用的、且簡單有效的數(shù)值方法。 求解線性方程組的數(shù)值方法盡管很多，但歸并起來可分為兩大類：（1） 直接法（精確法）凡經(jīng)有限次的四則運(yùn)算，若運(yùn)算中沒有舍入誤差即可求得方程組精確解的方法。如：克萊姆（Cramer）法則方法、消元法、分解法、分解法等等。（2） 迭代法（近似法）將求解方程組的問題轉(zhuǎn)化為構(gòu)造一個(gè)無限迭代的序列，在實(shí)現(xiàn)該序列過程中的每一步計(jì)算結(jié)果，均是把前一步所得的結(jié)果施行相同的計(jì)算步驟進(jìn)行修正而獲得的，而這一無限序列的極限就是原方程組的精確解答。如：簡單迭代法、賽德爾迭代法、

4、牛頓法、共軛斜量法等等。需要指出的，在一般情況下，我們使用直接法和迭代法兩類方法都不可能完全獲得原方程組的精確解答。原因很顯然：（1）實(shí)際中在使用直接法時(shí)不可能沒有數(shù)值計(jì)算的舍入誤差，故此時(shí)所謂精確方法的解并不是絕對精確的；（2）實(shí)際中在使用迭代法時(shí)，不可能將極限過程無限進(jìn)行到底，而只能進(jìn)行有限次的迭代，故獲得是滿足精度要求的近似解答。 關(guān)于這兩類方法求解的誤差分析，我們將在每類方法的介紹之后進(jìn)行簡要討論。§1.1 直接法Cramer法則與求逆方法設(shè)n元n個(gè)非齊次線性代數(shù)方程組為： (1.1)利用矩陣和向量符號，這個(gè)方程組可表為: (1.2)其中,方程組（1.1）的系數(shù)矩陣（階方陣）

5、； 方程組（1.1）的右端已知向量（階列陣）； 方程組（1.1）待求的解向量（階列陣）。§1.1.1 Cramer法則由線性代數(shù)中關(guān)于線性方程組解的定理可知：若A的行列式，則方程組（1.1）有唯一解。此時(shí)，根據(jù)Cramer法則，其解的表達(dá)式為： (1.3)其中，即將中的第j列元素依次換為右端已知向量的元素所構(gòu)成的階方陣的行列式。易見，利用Cramer法則給出的（1.3）式求解n階線性方程組（1.2）式，需要計(jì)算個(gè)階行列式，而每個(gè)n階行列式將有項(xiàng)，其中每一項(xiàng)又含有n個(gè)因子，故展開每一個(gè)n階行列式僅乘法運(yùn)算就需次（忽略加法運(yùn)算次數(shù)），而對個(gè)n階行列式，其乘法運(yùn)算的次數(shù)為：對式（1.3）其

6、除法次數(shù)為：n故求解方程組（1.2）其乘法和除法總的運(yùn)算次數(shù)為：次例如，若階，則次這是一個(gè)十分驚人的數(shù)字，即使利用超高速的電子計(jì)算機(jī)能夠勝任此計(jì)算次數(shù)，僅由于多個(gè)數(shù)的連乘亦有可能造成溢出而無法繼續(xù)運(yùn)算。因此，Cramer法則這個(gè)在理論上完善且精確的求解方法，僅在理論上和一些特殊情況下可以發(fā)揮作用，而對高階線性代數(shù)方程組的實(shí)際求解中幾乎沒有多少實(shí)用價(jià)值。為此，人們不得不研究其它一些計(jì)算簡單、且行之有效的求解方法。§1.1.2 求逆方法若，則A非奇異，即存在，則方程組（1.2）的解向量可表為： （1.4）常用的矩陣求逆方法有：（1）伴隨矩陣法；（2）初等變換法；但當(dāng)方陣A的階數(shù)n較大時(shí)，

7、求逆非常麻煩且計(jì)算量很大。故對高階線性代數(shù)方程組來說，求逆方法也是一種中看不中用的方法。§1.2 直接法Gauss（高斯）消去法盡管這是一種較為古老的方法，但至今仍不失為最常用和最有效的方法之一?；舅枷耄和ㄟ^逐次消元處理，將原方程組化為等價(jià)的三角形方程組進(jìn)行求解。§1.2.1 三角形方程組所謂三角形方程組無非是以下兩種形式的方程組：(1)下三角形式(2.1)矩陣記為 (2.1')其中系數(shù)矩陣的元素滿足關(guān)系：主對角線以上的元素均為零。（2）上三角形式 （2.2）矩陣記為 （2.2'）其中系數(shù)矩陣中主對角線以下的元素均為零，即：對三角形方程組的求解是十分簡單的

8、。顯然對于下三角方程組（2.1），其求解步驟如下： 從第一個(gè)方程中解得：； 將代入第二個(gè)方程中，從中解得：； 將代入第三個(gè)方程中，從中解得：；······如此逐個(gè)方程求解的過程向前遞推下去，直到第n步。 將代入第n個(gè)方程中，從中解得：。對上述過程其完整的求解計(jì)算格式可歸結(jié)為： （2.3）這一求解過程稱為前推過程。同理，對于上三角方程組（2.2），其求解步驟亦可如法泡制，即： 從第n個(gè)方程中解得：； 將代入第n-1個(gè)方程中，從中解得：； 將、代入第n-2個(gè)方程中，從中解得：；·····&

9、#183; 將代入第1個(gè)方程中，從中解得。對上述過程其完整的求解計(jì)算格式可歸納為： （2.4）這一求解過程稱為回代過程。通常是將（2.4）式合并為一式表為： （2.4a）式中約定：當(dāng)足標(biāo)j的取值大于其上界n時(shí)，和式。§1.2.2 Gauss消去法高斯消去法的求解過程就是首先利用矩陣的初等行變換方法將原方程組逐次消元，使之化為等價(jià)的具有上三角形式的方程組，然后再按上三角方程求解的計(jì)算格式（2.4a）式求出原方程組的解。整個(gè)求解過程可分為消元和回代兩個(gè)過程。1. 簡例 為了便于說明高斯消去法的求解過程，以如下4階方程組的求解為例： (2.5)其展開式為： （2.5）其增廣矩陣為： （2.

10、5.0）（1） 消元過程利用若干輪的初等行變換處理，設(shè)法將中的A部分化為上三角形式。若主元素，則可保留方程組中的第一個(gè)方程，并利用將其余三個(gè)方程中的第一個(gè)未知量消去，具體做法是取數(shù)：再對增廣矩陣（2.5.0）中的第i行進(jìn)行如下初等變換： （其中r表示行或排row）這樣原方程的增廣矩陣（2.5.0）被變換為如下形式： （2.5.1）其中被改變的元素為：至此，完成了第1輪初等行變換處理。若增廣矩陣（2.5.1）中的主元素，則可保留第1個(gè)和第2個(gè)方程，利用同上的處理方法消去第3個(gè)和第4個(gè)方程中的第二個(gè)未知量，即取數(shù)：再對增廣矩陣（2.5.1）進(jìn)行如下初等行變換：得到新的增廣矩陣如下： （2.5.2）

11、其中的被改變的元素為：至此，完成了第2輪初等行變換處理。若增廣矩陣（2.5.2)中主元素，則保留第13個(gè)方程，將第4個(gè)方程中的第三個(gè)未知量消去，即取數(shù)：再對增廣矩陣（2.5.2）進(jìn)行如下初等行變換：得到新的增廣矩陣為： （2.5.3）其中被改變的元素為：至此，完成了第3輪初等行變換處理。此時(shí)，原4階方程組通過3輪初等行變換之后，將原方程組（2.5）的增廣矩陣中的A部分已化為上三角形式。這一演變過程即為消元過程。（2） 回代過程由線性代數(shù)方程組理論可知：經(jīng)初等行變換之后得到的新的增廣陣（2.5.3）所對應(yīng)的方程組與原方程組（2.5）為等價(jià)方程組，即（2.5.3）的解即為原方程組（2.5）的解。由

12、前關(guān)于上三角方程組求解的計(jì)算格式（2.4a），可依次逆序求得原方程組的解為：逐步求解出的過程即為回代過程。2. Gauss法的求解步驟與計(jì)算格式通過以上4階方程組的求解過程，我們不難歸納出對于n階方程組Gauss消去法的求解步驟與計(jì)算格式。設(shè)n階線性方程組為： （1） 消元過程：通過（n-1）輪初等行變換處理將原方程組變?yōu)槿缦碌葍r(jià)的上三角方程組，即： （為上三角方陣）具體做法是：將記為 在第k輪的消元過程中，增廣陣中的前k行元素與增廣陣中前k行元素保持不變，對其他改變元素和的計(jì)算格式為： （2.6）（2） 回代過程：求解由消元過程形成的上三角方程組，即由（2.4a）式的合并表達(dá)式，可得解的計(jì)算

13、格式為： （2.7）式中約定：當(dāng)足標(biāo)j的取值大于其上界n時(shí)，和式。3. n階方程組高斯消去法的計(jì)算量（乘除法次數(shù)） 在消元過程中（公式（2.6）的計(jì)算中），即在第次消元執(zhí)行過程中需作（）次除法，次乘法，故在整個(gè)消元過程中共需乘、除法次數(shù)為：。在回代過程中（即公式（2.7）的計(jì)算中），共需作次乘、除法。則高斯消元法總的乘除法次數(shù)為：其中是二次方以下的余項(xiàng)，可見主要運(yùn)算量是花費(fèi)在消去過程中。 例對階，高斯消元法的乘除法運(yùn)算次數(shù)為，而Cramer法的乘除法運(yùn)算次數(shù)為。顯見，高斯消元法的計(jì)算量減少的非常顯著。§1.3 主元素消去法§1.3.1主元素選取的必要性 在前述的Gauss消

14、元法中，未知量是按其在方程式中的自然順序逐次消去的；用來消去其他方程式中未知量的方程亦是按自然順序選取的，故這種方法又稱為自然順序消元法。但實(shí)際計(jì)算結(jié)果表明：這種按自然順序消元法的方法在某些情況下可能具有如下兩個(gè)嚴(yán)重的缺陷：（1） 對，在消元過程中有可能出現(xiàn)主元素為零的情況，假設(shè)在第k主步的消元過程中，主元素，則有數(shù)（除法溢出），盡管，但此時(shí)消元過程將無法進(jìn)行下去，即導(dǎo)致消元過程就此失敗。（2） 若有主元素，但其絕對值已很小，此時(shí)消元過程仍可以繼續(xù)下去，但用絕對值很小的主元素作為除數(shù)，將會(huì)導(dǎo)致其他元素?cái)?shù)量級的嚴(yán)重增長和舍入誤差的擴(kuò)散，使得最終獲得的解答極不準(zhǔn)確。這里給出一個(gè)簡單的例子。例1.

15、二元線性方程組其精確解為：現(xiàn)取三位浮點(diǎn)十進(jìn)制數(shù)求解（即計(jì)算中保留小數(shù)后三位有效數(shù)字）：1 按自然順序消元，即有：由此得方程的近似解為：（其結(jié)果已完全失真）2 行序交換后再消元（即原兩個(gè)方程的位置對調(diào)），即有：由此得方程的近似解為：（與精確解相當(dāng)接近）?，F(xiàn)分析考察計(jì)算過程1產(chǎn)生解失真結(jié)果的原因： 令：（解1的誤差）將和分別代入1中的第一個(gè)方程中，則分別有： 上兩式相減得：即有：，這表明的誤差傳播給誤差時(shí)被反向擴(kuò)散了倍。上述例子表明：在消元過程中，適當(dāng)選取主元素是十分必要的。大量實(shí)踐表明，若A為對稱正定陣，順序消元法可以保證計(jì)算過程對于舍入誤差的數(shù)值穩(wěn)定性；但對于一般方陣A，則必須進(jìn)行選主元處理后

16、，消元法才能得到滿意的解答。§1.3.2 主元素消去法 該法僅是在高斯消元法中的每一消去主步之前增加一個(gè)主元素的選取過程，其余過程則完全相同。 最常用和最有效的主元素選取方法有兩種：1. 列主元素消去法 未知數(shù)仍按自然順序消去，但從各方程要消去的那個(gè)未知量的系數(shù)中選出模最大者作為下一步消元過程的主元素。 如擬消去未知數(shù)時(shí)，則在要消去的各方程中找出系數(shù)中模最大者，假如它是第（）個(gè)方程中的系數(shù)，則將第個(gè)方程與第k個(gè)方程交換位置，則使所選的主元素處于原來主元素的位置上，再按順序消去法計(jì)算公式進(jìn)行處理。例2 方程組 其增廣陣為：回代解得 2. 全主元素選取法 未知數(shù)不再按自然順序消去，而是在

17、每輪消去步之前，從所要消去未知數(shù)的所有方程中尋找出未知量系數(shù)模最大者作為當(dāng)前的主元素，并再通過行或列的調(diào)換將其置換至原來主元素所在位置處。 在第k步主步消元之前，須先從第k至n個(gè)方程中的第k至n個(gè)未知數(shù)的系數(shù)中找出模最大者，假如它是，則將第k個(gè)與第個(gè)方程交換位置，再將各方程中與的兩項(xiàng)位置交換（相當(dāng)于系數(shù)陣A及右端項(xiàng)中的第k行與行交換，A中的第k列與列交換）。此時(shí)所選主元素處于原來主元素的位置，然后再按順序消元法公式進(jìn)行處理。 例： 兩種主元消去法的優(yōu)缺點(diǎn)：（1） 列主元消去法程序簡單，在消去過程中可保持矩陣A的某些特征（帶狀，特征塊型），并有足夠的計(jì)算精度。（2） 全主元消去法計(jì)算精度高，但程

18、序稍復(fù)雜，在消去過程中將無法保持矩陣A的某些特征。§1.4 矩陣的分解法前面介紹了高斯消元法，其本質(zhì)就是將線性方程組的系數(shù)矩陣A分解為下三角形方陣和上三角形方陣的乘積。本節(jié)首先從矩陣運(yùn)算的角度來分析消元法，然后給出利用矩陣分解法求解方程組的步驟和計(jì)算公式。§1.4.1 消元法與矩陣分解設(shè)n階方程組有解，即（矩陣A非奇異）。 由前知，Gauss消去法的求解過程為：記為 ，or 其中：為上三角形n階方陣。即通過n-1輪消元主步處理將原方程組變成等價(jià)的上三角形方程組，即，而每一輪消元主步相當(dāng)于對原方程組的增廣陣進(jìn)行了若干次初等行變換。按初等行變換的作用，即為用相應(yīng)的初等方陣去左乘

19、被變換的矩陣。由此，容易驗(yàn)證： 第1輪消元主步過程相當(dāng)于對增廣陣進(jìn)行了如下的矩陣乘法運(yùn)算： ， ，其中，初等變換矩陣，是第1列對角線以下元素非零，其它元素與單位陣相同的下三角方陣。類似地，第2輪消元主步過程可等價(jià)表示為如下矩陣乘法：其中，初等變換矩陣是第2列對角線以下元素非零，其它元素與單位陣相同的下三角方陣。······第（n-1）輪消元主步等價(jià)為如下矩陣乘法：其中，初等變換矩陣是第(n-1)列對角線以下元素非零，其它元素與單位陣相同的下三角方陣。此時(shí)，即為上三角形方陣，為右端項(xiàng)由上述遞推關(guān)系，顯然有： （4.1）其中為第i輪消元主

20、步對應(yīng)的初等變換矩陣，其表達(dá)式為： （4.2）它是第i列對角線以下元素非零，其它元素與單位陣相同的下三角形矩陣。按逆陣的定義，容易得到下三角形矩陣的逆陣為： （4.3）所以，從(4.1)中第一式中可解得： （4.4）其中 為三角方陣 （由矩陣乘法規(guī)則容易驗(yàn)證） （4.5）即L為下三角方陣。由（4.4）式可知： （4.4'）這表明原n階方程組的系數(shù)陣A經(jīng)n-1次高斯消元處理之后即被分解為下三角矩陣L與上三角矩陣的乘積。這一結(jié)果揭示了高斯消元過程與矩陣分解之間的本質(zhì)關(guān)系，即為矩陣的三角形分解定理。定理4.1(三角形分解定理) 若n階方陣A的各階主子式，則A定可分解為下三角陣與上三角陣的乘積

21、，即：。將，代入方程組中，則原方程組可表為： （4.5）若令： （上三角方程組） （4.6）從而（4.5）可表為： （下三角方程組） （4.7）即原方程組的求解問題亦被相應(yīng)分解為相繼求解（4.7）下三角形系數(shù)的方程組和（4.6）上三角形系數(shù)的方程組問題。關(guān)于（4.7）和（4.6）的計(jì)算格式見前公式（2.3）和（2.4）。為表明其正確性，過程結(jié)果如下：即先從（4.7）中求出，即，代入（4.6）中，則有，從中解出：。這樣處理的好處是顯然的，它的突出優(yōu)點(diǎn)是：在求解具有相同系數(shù)矩陣A而右端項(xiàng)不同的方程組時(shí)（在結(jié)構(gòu)有限元中，對應(yīng)同一結(jié)構(gòu)受到不同的荷載工況作用），對系數(shù)矩陣A無需再作重新分解處理。如對另一

22、方程組：，順序求解 即可。§1.4.2 直接求解法由定理4.1，對方程組的系數(shù)矩陣進(jìn)行如下三角形分解：按矩陣乘法規(guī)則，A中的元素等于L中的第i行與中的第j列相乘的結(jié)果，即：當(dāng)時(shí)（下三角元素），有： （4.8.1）其中約定：當(dāng)足標(biāo)大于上界時(shí)，上和式恒為零。當(dāng)時(shí)（上三角元素），有： （4.8.2）分別從（4.8.1）式和（4.8.2）式中可獲得分解的計(jì)算公式為： （緊湊格式） （4.9） 在上式中，由矩陣乘法規(guī)則顯然有： 1° L的第1列元素即等于A中的第1列元素，即有： （4.9.1） 2° 的第1行元素為： （4.9.2） 不難獲得在分解公式（4.9）中計(jì)算和的各

23、元素的順序?yàn)椋?第1步，由（4.9.1）式求出中的第1列元素； 第2步，由（4.9.2）式求出中的第1行元素； 第3步，由（4.9）第1式求出中的第2列元素； 第4步，由（4.9）第2式求出中的第2行元素； 如此交替進(jìn)行，最終可求出三角陣和中的所有元素，即實(shí)現(xiàn)了對矩陣A的三角形分解。在三角形矩陣的分解公式(4.9)中，無須求出順序消去法中A陣的中間結(jié)果，而可直接由原矩陣A的元素算出和U的元素，故該法稱之為直接分解法。相應(yīng)的分解公式(4.9)亦稱為緊湊格式。當(dāng)?shù)娜切畏纸馔瓿珊?，則原方程組可轉(zhuǎn)化為逐次求解如下兩三角形方程組，即 (4.7) (4.6) 由計(jì)算格式（2.3）式，對下三角形方程組（4

24、.7）其解的計(jì)算格式為： （4.10）當(dāng)求出之后，由計(jì)算格式（2.4），對上三角形方程組（4.6）其解的計(jì)算格式為： （4.11）可以證明：上述的直接求解算法與高斯消元法是等價(jià)的，計(jì)算公式（4.10）和（4.11）分別對應(yīng)于消元法中的消元過程和回代過程。§1.5 對稱矩陣的分解法§1.5.1 平方根法（即法）設(shè)n階方程組，若系數(shù)矩陣A為對稱正定陣，即（對稱），（正定），則由矩陣的三角形分解定理4.1，A可被分解為： （5.1） （5.2）為下三角陣，則為上三角陣。將按矩陣乘法規(guī)則展開，由出發(fā)，按直接分解法類似的步驟，（5.1）式兩端元素之間關(guān)系為：其中約定：當(dāng)時(shí)，和式； 當(dāng)

25、時(shí)，和式。再由以上兩式，很容易獲得矩陣元素的計(jì)算式為： （5.3）此式即為對稱正定矩陣的平方根分解法。在此法中，由于下三角陣的第k列元素即為上三角陣的第k行元素，因此，我們只需順序地計(jì)算出陣的第1列至第n列元素即完成了對A陣的三角形分解，其分解的計(jì)算量僅為直接分解法的一半左右。若分解，則對原方程組的求解可分解為如下兩三角形方程組的順序求解:它們的求解公式為： 即將（4.11）式中的元素?fù)Q為即可。由式（5.3）可見，在平方根分解法中，需要進(jìn)行n次開方運(yùn)算，而在計(jì)算機(jī)上作開方運(yùn)算需要通過調(diào)用函數(shù)利用子程序來實(shí)現(xiàn)，這一則增大了運(yùn)算量，二則可能降低了計(jì)算精度。為了避免開方運(yùn)算的出現(xiàn)，人們在法的基礎(chǔ)上又

26、研究提出了一種改進(jìn)的分解法，即Cholesky（喬累斯基）法。§1.5.2 Cholesky（喬累斯基）法（法） 對，若A對稱正定，則A可被分解為： （5.4）其中 （5.5）為下三角陣，其主對角元素均為1，其余元素亦與前（5.2）式L陣略有不同（見后）；對角陣 （5.6）事實(shí)上，我們只需在前（5.1）的分解表達(dá)式中，對陣L稍作數(shù)學(xué)上的處理即可獲得（5.4）式結(jié)果。具體作法如下： 從（5.2）表出的下三角陣L的第j列元素中提出因子，由矩陣乘法，則有： （5.7）其中： （5.8）將（5.7）的關(guān)系式代入（5.1）式中，則有：（注意：此式L陣與（5.2）中L陣有所不同）仿照前面的處理，

27、可推得D陣和L陣元素的分解計(jì)算公式為： （5.9）此處約定：當(dāng)，時(shí)，和式，。對稱方陣的這種分解即為喬累斯基分解方法，顯見在此分解公式中不再含有開方運(yùn)算。利用喬累斯基分解法解方程的過程如下： （5.10） （5.11）對三角形方程（5.10）和（5.11）的求解過程和計(jì)算格式同前。§1.6 求解三對角方程組的追趕法在利用樣條函數(shù)進(jìn)行插值和利用有限差分法求解場問題的數(shù)值計(jì)算中，常常面臨著求解下列形式的三對角方程組：用矩陣表示 （6.1）其中，系數(shù)方陣：稱為三對角陣。其特點(diǎn)是：所有非零元素都集中在主對角及其相鄰的兩條次對角線上，而其余元素全為零。針對方程組系數(shù)陣的這種帶狀特點(diǎn)，人們給出高斯

28、消去法的一種簡化形式追趕法，其求解步驟：消元過程：第一步：將中的第1個(gè)方程中的系數(shù)化為1，即有：其中，第二步：注意到在的剩下的方程中，只有第2個(gè)方程含有，則利用式將方程中第2個(gè)方程中的消去，并將的系數(shù)化為1，即有其中，如此一步步順序加工方程組中的每個(gè)方程，設(shè)經(jīng)第步后，方程中的第個(gè)方程被化為：其中，當(dāng)經(jīng)第步之后，方程中的第個(gè)方程被化為：其中，將與中第個(gè)方程（最后一個(gè)方程）聯(lián)立，可解得變量為：其中，至此，通過上述消元過程，使原方程組化為如下形式更為簡單的等價(jià)方程組：即 回代過程：此方程組為二對角形的，其系數(shù)陣中的非零元素均集中分布在主對角線以上的兩條對角線上，且主對角元素均為1。顯然，通過自下而上

29、的逐步回代，即可從方程中依次求得解，其求解公式為：其中，按式計(jì)算，按式計(jì)算。上述算法即為追趕法，它的消元過程和回代過程分別為追的過程和趕的過程。具體地說，追的過程就是按算式和順序地求出系數(shù)和；趕的過程就是按算式逆序地求出解?？梢?，追趕法的消元原理與高斯法是相同的，但它針對三對角方程組的具體特點(diǎn)進(jìn)行了相應(yīng)的處理，即將系數(shù)陣中大量的零元素撇開，從而大大節(jié)省了計(jì)算量和存貯量，提高了算法的有效性。§1.7 簡單迭代方法（Jacobi雅可比法）迭代法的基本思想：將求解方程組歸結(jié)為重復(fù)計(jì)算一組彼此獨(dú)立的線性表達(dá)式，產(chǎn)生一個(gè)解向量序列，并使該向量序列收斂至某個(gè)極限向量，則即為方程組的準(zhǔn)確解。迭代法

30、的核心是建立相應(yīng)的迭代計(jì)算格式。對于求解線性方程組常用的方法有：簡單迭代法、Seidel（塞德爾）迭代法和松弛法等等。它們之間的差異僅在于迭代格式有所不同。本節(jié)介紹簡單迭代方法。§1.7.1求解方法和迭代格式以3階方程組為例說明迭代法的計(jì)算過程：設(shè)對角主元，將第個(gè)方程的第個(gè)未知量用其余項(xiàng)表出，則被化為等價(jià)形式：其中：式即為原方程組的等價(jià)方程，同時(shí)它又是簡單迭代的計(jì)算格式?，F(xiàn)任取一組方程組的近似解（通常稱為初值或零次近似解）：代入式的右端，即完成了一次迭代，可求出一組新的近似解（一次近似解），即有：將一次近似解再代入式的右端，即完成了第二次迭代，求得二次近似解為：。如法重復(fù)上述過程，經(jīng)

31、次迭代后，得到方程組的次近似解為：式即為簡單迭代法的迭代計(jì)算公式，若引入矩陣符號：令 則迭代公式可以矩陣表為：（7.4）式中：B迭代矩陣；常數(shù)列向量。（7.4）式屬于一階線性定常迭代形式。一階：新解僅是其前面一次近似解的函數(shù)；線性：新解為其前面近似解的線性函數(shù)；定常：無論為何值，計(jì)算的規(guī)則不變（即均為定常矩陣），且與迭代次數(shù)無關(guān)。按照3階方程組簡單迭代法求解過程和計(jì)算格式，容易寫出n階方程組的簡單迭代公式為：等價(jià)變形為（從中得的形式解表達(dá)式）：其中：將縮寫為：或統(tǒng)一表為：從上式可直接寫出簡單迭代法的迭代計(jì)算格式（表達(dá)式）為： 的矩陣形式完全同（7.4）式，即：其中，例1.（華中計(jì)算方法P208

32、）用簡單迭代法求解：依次從三個(gè)方程中分離出，建立迭代計(jì)算格式：（a）取初值，迭代結(jié)果如下：000010.720.830.8420.9711.071.1531.0571.15711.248281.09981.199411.2997891.099941.199941.29992已趨向于本問題的準(zhǔn)確解：。可見由迭代格式（a）所獲得的解序列最終將收斂于準(zhǔn)確解。§1.7.2 迭代的收斂條件簡單迭代法雖然求解過程簡單，但并不是隨便構(gòu)造一個(gè)迭代計(jì)算格式就能實(shí)現(xiàn)的，其中有一個(gè)非常重要的理論問題在迭代中必須予以關(guān)注，這就是所謂的收斂性問題?，F(xiàn)再一次考察例1，與前處理不同，先將分別從原方程組（a）的第2

33、，3和1個(gè)方程中分離出來，并建立相應(yīng)的迭代格式，即：（a）仍取初值迭代結(jié)果： 若繼續(xù)迭代下去，所得解向量的值會(huì)越來越大，即迭代格式（a）是發(fā)散的（不收斂）。這里首先給出迭代收斂性的有關(guān)數(shù)學(xué)描述。收斂性定義：對迭代格式（7.10）當(dāng)序列極限存在，則稱該迭代格式為收斂的，且即為原方程組的解。在以上兩例中，為何對同一方程組，構(gòu)造的兩個(gè)迭代計(jì)算表達(dá)式有所不同，會(huì)導(dǎo)致迭代結(jié)果收斂和發(fā)散兩種截然不同的結(jié)果？其原因何在？進(jìn)一步，欲確保迭代格式具有收斂性的條件是什么？收斂條件討論：將準(zhǔn)確解代入迭代公式中，必有兩端相等，即由式-，得：顯然有：令：即將矩陣中各行向量一階范數(shù)中的最大者記為。（注：假如向量為：，則其

34、一階范數(shù)為：） （7.14a） （7.14b）即準(zhǔn)確解與第次近似解和第次近似解之間的最大誤差。由（7.14a）和（7.14b）式，則式可表為： 此式即為第次迭代誤差與第次迭代誤差之間的關(guān)系不等式。由此可知，若，即使，必有當(dāng)，即迭代解向量的計(jì)算誤差將隨迭代進(jìn)程趨于零。這表明此迭代過程（格式）是收斂的。由此，可得迭代收斂性的判別定理如下：定理7.1 在迭代格式中，若，則該迭代格式對任意給定的初值均是收斂的。這一定理在實(shí)際中的應(yīng)用是特別簡單和方便的。這里，我們用此定理判別例1中兩種迭代格式的收斂性：對迭代式（a）：（收斂）對迭代式（a）：（不收斂）這一判別結(jié)果與前面的數(shù)值計(jì)算結(jié)果是完全一致的。

35、67;1.8 Gauss-Seidel（高斯-塞德爾）迭代法與松弛法簡單迭代法雖然具有迭代格式特別簡單的優(yōu)點(diǎn)，但其迭代的收斂速度較慢。為了加快迭代收斂速度，人們在簡單迭代法的基礎(chǔ)上研究提出了幾種改進(jìn)的迭代方案，主要有Gauss-Seidel迭代法和松弛法。§1.8.1 Gauss-Seidel法在簡單迭代格式式中將第一式求得的階近似解代替第二式右端的第階近似解去求得近似解；再將、代替第三式右端的、去求出，逐次這樣處理直到第n個(gè)方程。第輪迭代的計(jì)算式為：它是在簡單迭代法的計(jì)算過程中，逐次以新值代替老值再進(jìn)行迭代。式即是所謂的Gauss-Seidel迭代法的計(jì)算公式。為了寫出其一般迭代公

36、式，先不妨詳細(xì)寫出中的第個(gè)表達(dá)式如下：將其中新值部分和老值部分分別用和式形式表為，則有： 若將的關(guān)系式代入，則式可表為： （8.2）式或（8.2）式即為Gauss-Seidel迭代法的計(jì)算公式。由于新值比其老值更準(zhǔn)確一些，顯然用這些已知的新值去代替老值進(jìn)行迭代，自然迭代收斂于其準(zhǔn)確解的速度要比簡單迭代法來的快。為了進(jìn)行比較，仍以上節(jié)例1為例。例3 方程組同上節(jié)例1，將例1的迭代格式改寫為Gauss-Seidel迭代格式如下：其中迭代矩陣與簡單迭代法中的完全相同。取初值，迭代結(jié)果如下：000010.720000.902001.1644021.043081.167191.2820561.09999

37、1.199991.30000此結(jié)果與前簡單迭代結(jié)果相比可見，Seidel法的第6輪迭代結(jié)果的精度已高出簡單迭代法第9次迭代結(jié)果的精度，即Seidel法的迭代收斂速率明顯加快。顯然，在Seidel迭代公式中，其迭代矩陣與簡單迭代法中的完全相同。故關(guān)于簡單迭代法的收斂性及其判別定理7.1亦完全適用于Seidel迭代法。§1.8.2 松弛法松弛法是另一種線性加速的方法。其基本作法：將前一步的結(jié)果與Seidel方法的迭代值適當(dāng)進(jìn)行線性組合，構(gòu)造一個(gè)收斂速度較快的近似解序列，其迭代計(jì)算過程的每一步分為迭代和加速兩個(gè)步驟。其迭計(jì)算公式為：迭代：加速：其中，系數(shù)稱為松弛因子（即新解與老解的線性組合

38、系數(shù)）。我們可將迭代公式和加速公式合并在一起，則松弛法的迭代計(jì)算格式為：這就是松弛法的計(jì)算公式。顯然，當(dāng)取松弛因子時(shí)，式就變成了Gauss-Seidel法迭代公式。因此，適當(dāng)選取因子可得到高于Gauss-Seidel法的收斂速度。其加速作用的簡單解釋為：若迭代收斂，則式的計(jì)算結(jié)果定比更精確些，因此在加速公式中應(yīng)加大在線性組合中的比重，盡可能擴(kuò)大它的影響效果。為了達(dá)到這一目的，通常應(yīng)取松弛因子。但的值并不是越大越好，在實(shí)際計(jì)算中，可根據(jù)系數(shù)矩陣的性質(zhì)，結(jié)合經(jīng)驗(yàn)通過多次試算來選定最佳的松弛因子之值。數(shù)學(xué)上可以證明：為保證迭代格式式收斂，須要求：。當(dāng)松弛因子，稱為超松弛法；當(dāng)松弛因子，即為Gauss

39、-Seidel法；當(dāng)松弛因子，稱為低松弛法。在松弛法中，關(guān)于“松弛”一詞的數(shù)學(xué)作用，以及松弛因子取值范圍的證明可參見計(jì)算方法（武大、山大P69）。因涉及到矩陣譜半徑和譜范數(shù)等概念及其有關(guān)的定理，故在此不展開討論。Chapter 2插值與數(shù)值微分§2.1 插值問題概述§2.1.1 背景工程上：在機(jī)械制造、產(chǎn)品幾何造型中，常常給定一批離散的樣點(diǎn)，要求利用數(shù)學(xué)的方法，自動(dòng)構(gòu)造一條光滑連續(xù)的曲線（面），使該曲線（面）通過這些樣點(diǎn)，以滿足幾何造型設(shè)計(jì)要求，或者據(jù)此進(jìn)行機(jī)械加工或放樣。例如：賦形拋物面天線（非標(biāo)準(zhǔn)的拋物面），汽車、飛機(jī)的外殼等等。數(shù)學(xué)上：在利用數(shù)值方法求解有關(guān)連續(xù)函數(shù)問

40、題時(shí)，往往需要根據(jù)給定離散樣點(diǎn)的數(shù)據(jù)插補(bǔ)中間點(diǎn)。如：用差商代替微商；結(jié)構(gòu)有限元中求非節(jié)點(diǎn)處的位移和單元非形心處的應(yīng)力等等。與此相關(guān)的數(shù)學(xué)問題就是插值問題。簡單通俗地講，所謂插值就是在給定的樣點(diǎn)間再插進(jìn)一些所需的中間值。§2.1.2 數(shù)學(xué)描述給定一組樣點(diǎn)：，設(shè)法構(gòu)造一個(gè)光滑連續(xù)的函數(shù)y=F（x），使該函數(shù)曲線（面）通過所有樣點(diǎn)，即滿足： （*）如圖示： 滿足插值條件（*）的函數(shù)F(x)即為所求的插值函數(shù)。構(gòu)造F（x）的基本要求：（1）函數(shù)形式應(yīng)為便于計(jì)算的初等函數(shù)，如多項(xiàng)式函數(shù)、分段多項(xiàng)式函數(shù)、樣條函數(shù)等等；（2）F(x)中的自由度（即待定參數(shù)的個(gè)數(shù)）應(yīng)與插值條件（*）的個(gè)數(shù)相等。為以

41、下討論方便起見，統(tǒng)一約定如下：樣點(diǎn)的節(jié)點(diǎn)；插值點(diǎn)；某個(gè)原函數(shù)（通常未知）；函數(shù)在節(jié)點(diǎn)處的函數(shù)值和一階導(dǎo)數(shù)值（已知）；待求的插值函數(shù)（待求）；原函數(shù)與插值函數(shù)在點(diǎn)處的插值余項(xiàng)（絕對誤差）；為含有點(diǎn)（）的最大區(qū)間。若插值點(diǎn)，則該插值問題稱為內(nèi)插；若插值點(diǎn)，則該插值問題稱為外推。§2.2 Lagrange（拉格朗日）插值若給定樣點(diǎn)，使（過點(diǎn)），則稱為Lagrange（拉格朗日）插值函數(shù)。下面首先介紹拉氏插值中兩種最簡單的情況，即線性插值和拋物線插值，然后在此基礎(chǔ)上給出高次拉氏插值函數(shù)的表達(dá)式。§2.2.1 線性插值（兩點(diǎn)一次插值）。求：過兩樣點(diǎn)（x0，f0），（x1，f1）的線性

42、插值函數(shù)y=F（x）。由解析幾何關(guān)系，過、兩點(diǎn)的線性函數(shù)（直線方程）可表為： (點(diǎn)斜式) （2.1a）上式經(jīng)整理，則有： (對稱式) （2.1b）為今后表述方便且統(tǒng)一起見，令： (2.2)這是兩個(gè)的一次多項(xiàng)式函數(shù)。顯然有： （2.3），在節(jié)點(diǎn)、處構(gòu)成了正交規(guī)范基，故稱之為線性插值函數(shù)的基函數(shù)。從而（2.1b）式給出的線性插值函數(shù)可由基函數(shù)表為： （2.1c）此式表明：插值函數(shù)可表為基函數(shù)的線性組合，即曲線是由基函數(shù)張起構(gòu)成的。利用微分中值定理，可推得線性插值函數(shù)的插值余項(xiàng)為： （2.4）式中：為微分中值點(diǎn)；=為的二次多項(xiàng)式?？梢姡壕€性插值函數(shù)對于原函數(shù)具有二階精度，即對于原函數(shù)的一階導(dǎo)數(shù)值亦有

43、逼近性。當(dāng)原函數(shù)為線性函數(shù)，則插值余項(xiàng)，即此時(shí)插值函數(shù)就是原函數(shù)。§2.2.2 拋物線插值（三點(diǎn)二次插值）求過三樣點(diǎn)（x0，f0），（x1，f1），（x2，f2）的二次多項(xiàng)式（拋物線）插值函數(shù)。由前知，所求插值函數(shù)可以表為其對應(yīng)基函數(shù)的線性組合形式，故不妨假設(shè)： （2.5）其中，為待求的基函數(shù)，它們應(yīng)滿足基函數(shù)的正交規(guī)范化條件，即有： （2.6）這樣的基函數(shù)不難直接構(gòu)造，現(xiàn)以為例說明之。欲使?jié)M足，顯然在的表達(dá)式中應(yīng)含有和兩個(gè)因子，故令： （2.7）其中為待求系數(shù)。由條件，解得：代入（2.7），則有 （2.8a）這是一個(gè)二次的多項(xiàng)式函數(shù)。同理可得： （2.8b） （2.8c）將（2.8

44、）式統(tǒng)一表為： （2.8）將（2.8ac）代入（2.5）式中，即得所求插值函數(shù)F(x)。易驗(yàn)證F(x)滿足插值條件：利用微分中值定理，二次插值函數(shù)的插值余項(xiàng)為： （2.9）即具有三階精度，直到二階導(dǎo)數(shù)都有逼近性。 例1 利用100,121,144的平方根求（精確值=10.7238*）之值。解：已知樣點(diǎn)為：（x0=100，f0=10），（x1=121，f1=11），（x2=144，f2=12），插值點(diǎn)：x=115，求插值函數(shù)值F（115）。代入插值公式（2.5）中，即有：可見誤差甚小。§2.2.3 n次多項(xiàng)式插值求過（n+1）個(gè)樣點(diǎn)（xi，fi）（i=0,1,2,n）的n次多項(xiàng)式插值函

45、數(shù)。設(shè)： （2.10）由前討論可知，該n次多項(xiàng)式函數(shù)共用（n+1）項(xiàng)，它可由其（n+1）個(gè)基函數(shù)的線性組合形式表出，即有： （2.11）其中（n+1）個(gè)基函數(shù)必須滿足正交規(guī)范條件，即： （2.12）由此條件仿前處理，可推得基函數(shù)的表達(dá)式為： （2.13）將（2.13）式代入（2.11）式中，得n次多項(xiàng)式插值函數(shù)的計(jì)算表達(dá)式為： （2.14）可見，拉氏插值公式在邏輯結(jié)構(gòu)上表現(xiàn)為二重循環(huán)，其中內(nèi)循環(huán)是由累積獲得基函數(shù)，外循環(huán)由節(jié)點(diǎn)的函數(shù)值與基函數(shù)之積的累加獲得插值點(diǎn)處的函數(shù)值。利用微分中值定理，可得n次拉氏插值的余項(xiàng)為： （2.15）該插值函數(shù)具有（n+1）階精度。例2 已知原函數(shù)的四個(gè)樣點(diǎn)值為：

46、0.10.150.250.30.9048730.8607080.7788010.740818求插值點(diǎn)處的f（x）的近似值。（精確值）解：因?yàn)榻o定四個(gè)樣點(diǎn)，故采用四點(diǎn)三次拉氏插值公式。由（2.11）式，即有： （a）由（2.14）式，其中基函數(shù)為： （b）將和x=0.2分別代入（b）中，得基函數(shù)在x=0.2處的值為：由（a）式x=0.2處的插值函數(shù)值為： 由插值余項(xiàng)公式（2.15）式，則有：從而在插值點(diǎn)處的余項(xiàng)的估計(jì)為：由于為單減函數(shù)，故。從而得的估計(jì)區(qū)間為：在處插值余項(xiàng)的精確值為：落入誤差的估計(jì)區(qū)間內(nèi)。§ 2.3 Hermite（埃爾米特）插值在許多情況下，不僅要求插值函數(shù)過所有的樣

47、點(diǎn)，甚至還要求與在樣點(diǎn)處的導(dǎo)數(shù)值彼此相等，即要求既“過點(diǎn)”還要求“相切”。這就是所謂的埃爾米特（Hermite）插值問題。其數(shù)學(xué)表述為：若使：， 且使： 則稱為Hermite插值函數(shù)。本節(jié)首先介紹埃爾米特插值中兩種最簡單情況，即一點(diǎn)一次帶導(dǎo)數(shù)插值和兩點(diǎn)三次帶導(dǎo)數(shù)插值，然后再給出一般高階埃氏插值的計(jì)算公式。§2.3.1 一點(diǎn)一次帶導(dǎo)數(shù)插值（切線插值）給定一樣點(diǎn)為，其中。求一次插值多項(xiàng)式函數(shù)F(x)，滿足：， 根據(jù)以上兩個(gè)插值條件，由解析幾何，構(gòu)造為： （3.1）顯然有： 利用微分中值定理，其插值余項(xiàng)為： （3.2）該插值函數(shù)具有一階導(dǎo)數(shù)的逼近性。§2.3.2 兩點(diǎn)三次帶導(dǎo)數(shù)插

48、值給定兩樣點(diǎn)：，求既過點(diǎn)又相切的多項(xiàng)式插值函數(shù)。即插值條件為： （*）共有四條，故設(shè)F(x)為三次多項(xiàng)式函數(shù)（共有四項(xiàng)），即 （3.3）將（*）式代入上式，得以a0，a1，a2，a3為未知量的線性方程組為：聯(lián)立解得a0，a1，a2，a3，代入插值函數(shù)（3.3）式中，并經(jīng)整理得的統(tǒng)一表達(dá)式為： （3.4）其中：即為基函數(shù)，它們表達(dá)式為： （3.5）其中：即為兩點(diǎn)線性插值的基函數(shù)；為節(jié)點(diǎn)步長。易驗(yàn)證： 這是一組四維正交規(guī)范基，這表明：所求插值函數(shù)（3.4）式即為基函數(shù)的線性組合。這與拉氏插值的結(jié)果是一致的，差別僅在于基函數(shù)有所不同。由微分中值定理，（3.4）式的插值余項(xiàng)為： （3.7）=該插值函數(shù)

49、具有三階導(dǎo)數(shù)的逼近性。對于基函數(shù)（3.5）的另一種表示因: x的二次多項(xiàng)式 (3.6)則有： 從而基函數(shù)亦可表為： （3.5a）此種形式便于推廣到高次埃爾米特插值多項(xiàng)式函數(shù)的基函數(shù)表達(dá)式之中。§2.3.3 高階埃氏插值給定(n+1)個(gè)樣點(diǎn)，求既過點(diǎn)又相切的多項(xiàng)式插值函數(shù)。插值條件為： （3.8）共有個(gè)，故插值函數(shù)應(yīng)為次多項(xiàng)式函數(shù)。仿照前兩點(diǎn)三次帶導(dǎo)數(shù)的埃氏函數(shù)的產(chǎn)生方法，這個(gè)次多項(xiàng)式插值函數(shù)可表為： （3.9）其基函數(shù)，（），可按（3.5a）式的格式表為： （3.10）其中 （3.11） （3.12）其插值余項(xiàng)為： （3.13）具有階導(dǎo)數(shù)的逼近性。埃氏插值的幾何意義是：插值函數(shù)()與原函數(shù)(x)在節(jié)點(diǎn)處即“等值”又“公切”。埃氏插值還可推廣到包含更高階導(dǎo)數(shù)以及各節(jié)點(diǎn)包含的導(dǎo)數(shù)階不均等的情況。§1.4 關(guān)于多項(xiàng)式插值的穩(wěn)定性根據(jù)拉氏和埃氏插值的余項(xiàng)公式(2.15)和(3.13)，似乎表明插值函數(shù)的次數(shù)愈高則精度愈高，與原函數(shù)的貼近效果就愈好。其實(shí)并非如此，現(xiàn)分析如下：在插值過程中總的誤差來自于以下兩個(gè)方面：1. 原函數(shù)f（x）被插值函數(shù)F（x）代替所引起的原理性誤差，即余項(xiàng)E(x)=f(x)-F(x)（截?cái)嗾`差），其誤差分析即為余項(xiàng)公式。2. 節(jié)點(diǎn)數(shù)據(jù)fi本身的誤差（實(shí)驗(yàn)誤差、