計算方法 第3章 線性方程組數(shù)值解法

第三章 線性方程組的數(shù)值解法 線性方程組是應用最為廣泛的數(shù)學模型，很多復雜問題中都含有線性方程組子問題，因此討論線性方程組問題的求解很有必要，本章將討論線性方程組的數(shù)值解法。 線性方程組的一般形式： (3.1)如果記： 這里A稱為系數(shù)矩陣，x稱為解向量，b稱為右端項，則得線性方程組的矩陣形式： (3.2) 求解線性方程組問題(3.1)或(3.2)的數(shù)值方法可分為兩類：直接解法和迭代解法，其直接解法是通過有限次初等運算，求得其解，雖然直接解法的推導過程都是無誤差的，但是由于計算機的運算都是有舍入誤差的，所求解其實是一個有誤差的近似解；迭代解法則是從某個初始近似解出發(fā)，按照一個確定的迭代公式得到一個更好的近似解，反復迭代，直到求得一個滿足精度要求的近似解。 本章首先討論線性方程組的直接解法然后再介紹迭代解法。3.1 消去法1.順序Gauss消去法首先回顧一下線性代數(shù)中所講的線性方程組消去法過程，然后歸納出消去法的數(shù)值算法，請看如下的例子：求解線性方程組解：求解線性方程組的第一階段稱為消元過程，其方法是：第2個方程減去第1個方程的倍，第3個方程減去第1個方程的2倍，得第3個方程減去第2個方程的倍，得這一過程就是消元過程，即把方程化為等價的上三角方程(對角線下變?yōu)?)。第一個階段完成后，進入第二個階段，稱為回代過程，其方法是：先由第3個方程解出，將代入第2個方程解出，再將和代入第1個方程解出也就解出所有的未知量。如下就是所求的解： 歸納以上求解方法，求解線性方程組包括兩個過程，消元過程和回代過程。首先給出回代過程的算法，回代過程其實是一個特殊形式的方程組的求解方法，就是一個上三角線性方程組的求解方法，如： (3.3)第1步：根據(jù)第n個方程解出，得 第2步：根據(jù)已求出的和第n-1個方程求，得 假如已經(jīng)求出，代入第k個方程可以求出，得： (3.4)回代過程就是對實施這一公式，注意必須從后向前計算方可，所以此過程叫做回代過程，(3.4)式也是求解下三角方程的一般公式，可看作一個獨立的算法。 算法3.1(上三角線性方程組的回代算法)[初始化]設置上三角方程系數(shù)矩陣A，右端項向量b[回代過程]對于循環(huán)對于循環(huán)[算法結束]再看Gauss消去法的消元過程，對于線性方程組的消去法本質(zhì)上可看作將其增廣矩陣用初等行變換化為梯形矩陣的過程。為清楚的表示每次消元前后系數(shù)矩陣和右端項的狀態(tài)，通常以和表示系數(shù)矩陣和右端項，其元素分別記作，其上標表示在第k次消元前的狀態(tài)，其初始增廣矩陣為： (3.5)首先對第1列消元，也就是對增廣矩陣做初等行變換將元素約化為0，希望消元之后形如： (3.6*)為此，對，令消元因子，將第i行減去第1行的倍，則可將(3.5)式變?yōu)?3.6*)式，之后的論述中，在不引起混淆的情況下，將(3.6*)式記為： (3.6)第2列消元，對，令，將第i行減去第2行的倍，則有對第3、4、…、n-1列消元，最后得到上三角方程組的增廣矩陣為以上的消元過程可用一個統(tǒng)一公式表示為： (3.7)通常在編程求解線性方程組時，消元因子存儲于矩陣對角線以下相應的位置上，因為此時這些位置的元素將被消元為0。算法3.2(順序Gauss消去法)[初始化]置系數(shù)矩陣A，右端項向量b；[消元過程]對于循環(huán)對于循環(huán)[回代過程]執(zhí)行算法3.1，即[算法結束]定理3.1順序Gauss消去法能夠順利進行的充要條件是系數(shù)矩陣A的順序主子式，且，特別有。證：從前述的算法可知，消元過程能夠順利進行的充要條件是步驟1.1)中對角線元素，以下只要證明順序主子式為即可。事實上，順序消元過程所做行變換并不改變順序主子式的值，而在完成第k次消元過程后，增廣矩陣被約化成：顯然第k個順序主子式為：2．順序Gauss算法的計算量通常我們在度量一個有限步終止的算法的計算量時，主要考慮整個計算過程中使用乘除法運算的總次數(shù)，在以后統(tǒng)計算法的計算量時，我們僅統(tǒng)計算法中乘除運算的次數(shù)。消元過程的計算量 回代過程的計算量 總的乘除法運算的計算量為總計算量=可以估計出總的加減法運算的計算量也約為列主元Gauss消去法前節(jié)所講順序Gauss消去法，有以下兩個問題：在消元過程中當某個元素對角線為零時，算法無法進行；如果某個對角線元素的絕對值非常小，將會引起很大的誤差。因此順序Gauss消去法不具實用性，但由此方法改進的列主元消去法是最為廣泛使用的方法，以下我們通過一個實例體會一下順序Gauss消去法的問題和稍作改進之后的效果：用順序Gauss消去法求解線性方程組(運算過程中保留4位有效數(shù)字)解：寫出方程組的增廣矩陣消元過程：，第2行減去第1行的倍，得回代過程得到計算解為：，此方程的準確解為，顯然此計算解已經(jīng)沒有任何意義。如果將兩個方程交換一下位置再求解，增廣矩陣：消元過程：，第2行減去第1行的倍，得回代過程得到計算解為：，與精確解完全相同！ 第1種解法產(chǎn)生非常大的誤差的原因是的絕對值太小，做為分母將產(chǎn)生太大的誤差而導致最后的解面目全非；用第2種方法是通過交換方程的位置使做分母的的絕對值盡可能的大，正如第一章誤差分析所述，避免絕對值太小的數(shù)做分母，可以提高計算結果的精度。這一過程稱為選主元(所選的對角線元素稱為主元)，一般的做法如下。 設進行第k步計算時增廣矩陣如下：為了在第k次消元過程中使消元因子中的分母絕對值盡可能的大，在第k列中取第k到第n個元素中絕對值最大者： (3.8)在選定主元之后，考慮如下情況如果，則線性方程組的行列式的值為0，其解不定(無解或有無窮多組解，停止計算)；如果，則的絕對值最大，以作為主元；如果則交換矩陣的第行，交換之后第k行對角線元素絕對值變?yōu)樽畲?。而交換增廣矩陣的兩行相當于交換兩個方程的位置，顯然不改變原方程的解；另一個需要指出的是，由于誤差的原因判斷的方法采?。簩τ谥付ǖ某浞中〉?，如果，則認為。增加選主元機制的Gauss消去法即為列主元Gauss消去法，列主元Gauss消去法的計算過程包括：選主元，交換主元行，消元。以下我們寫出適合編程的Gauss算法。算法3.3(列主元Gauss消去法)[初始化]置系數(shù)矩陣A，右端項向量b，取充分小的數(shù)；[消元過程]對于循環(huán)[選主元]如果，停止計算，算法失敗如果交換方程，[消元]對于循環(huán)1.4.1)1.4.2)1.4.3)[回代過程]執(zhí)行算法3.1，即計算[算法結束]關于列主元Gauss消去法的計算量，從列主元Gauss消去法的算法可以看出，僅比順序Gauss消去法多找主元部分，這一部分的操作主要是比較運算及必要時交換兩個方程，其他的部分算法完全一樣，因此影響不大。全主元Gauss消去法全主元Gauss消去法與列主元Gauss消去法的區(qū)別在于取主元的方法是如果交換增廣矩陣的第行，如果交換增廣矩陣的第列(本質(zhì)上是交換未知量與，計算結束時應該交換回來)。顯然全主元消去法取主元的方法范圍更大，主元的絕對值會更大，計算過程也就更可靠更穩(wěn)定，但計算量稍大。實際計算效果表明，列主元Gauss消去法已經(jīng)很穩(wěn)定，全主元法的改進甚微，是沒必要的，所以通常使用列主元Gauss消去法。Gauss-Jordan消去法列主元Gauss消去法在消元過程中是將對角線以下的元素約化為0，將系數(shù)矩陣化為上三角矩陣，通過回代得到方程的解；Gauss-Jordan消去法是將系數(shù)矩陣約化成單位矩陣，無需回代就直接得到了方程組的解，也稱為無回代過程的消去法。設線性方程組的增廣矩陣為： 第1步計算，設，第1行除以，第行減去新第1行的倍，則增廣矩陣變?yōu)槿缦滦问?第2步計算，設，第2行除以，第行減去新第2行的倍，則增廣矩陣變?yōu)槿缦滦问皆O第k-1步計算結束時增廣矩陣為如下形式 第k步計算，設，第k行除以，第k)行減去新第k行的倍，則增廣矩陣變?yōu)橥瓿蒼-1步計算之后，增廣矩陣變?yōu)轱@然最后一列就是方程組的解。 注意到在求解過程中，總是假設對角線上的主元，若算法將無法進行下去，如同列主元消去法，在消元過程中也可用選列主元的策略來避免這種情況發(fā)生，除非系數(shù)矩陣奇異。 算法3.4(求解線性方程組的Gauss-Jordan消去法)[初始化]置系數(shù)矩陣A，右端項向量b，取充分小的數(shù)；[消元過程]對于循環(huán)[選主元]如果，停止計算，算法失?。∪绻粨Q方程 [規(guī)格化主元行][消元]即對于循環(huán)1.5.1)1.5.2)[算法結束]6.Gauss-Jordan消去法求逆矩陣及算法設計 首先通過一個求逆過程來寫出相應算法。 例3.3用Gauss-Jordan消去法求下矩陣的逆陣： 解：在矩陣A的右邊添加單位矩陣進行初等行變換此時對應初始狀態(tài)的單位矩陣就是所求的逆矩陣：通常編程實現(xiàn)求逆矩陣時并不添加右邊的單位矩陣，而是對原始矩陣直接進行變換，在完成計算時原矩陣就變成了逆矩陣。其實我們注意一下求逆的過程，有如下兩個特征：首先、從初始狀態(tài)到最終狀態(tài)，增廣矩陣的每個狀態(tài)都有n個列是單位列向量，這也就意味著每次變換，都會將原矩陣的一列變?yōu)閱挝涣邢蛄浚鴮⒃瓎挝痪仃嚨囊粋€列向量變?yōu)槟婢仃嚨囊涣?，每次變換時可讓這一對列向量占用共同的列，也就是可將右邊新產(chǎn)生的一列直接放在左邊將要變?yōu)閱挝涣邢蛄康奈恢蒙?。其次、在求逆過程中如果沒有行交換，則增廣矩陣中單位列向量的消失與新增，在左右兩部分都是按自然順序遞增的。從例3.3可以看出，因交換主元行，右邊的單位列向量順序隨主元行而變，因此在完成計算后需要換回原位置，這正說明原矩陣的行交換對應逆矩陣的列交換這一性質(zhì)。以下過程是通常編程時的實際計算過程：整個計算有個消元過程，但是最后一次消元不須選主元，所以主元行標分別是{2,3}，按逆序分別交換列，則得到正確的逆矩陣，此過程直觀性太差，僅作為調(diào)試程序時核對之用。 算法3.5(矩陣求逆)[初始化]輸入矩陣A，并以向量p記錄主元行下標(共n-1個主元行下標)，取充分小的數(shù)；[行變換]對循環(huán)[選主元]，如果，矩陣不可逆，算法失?。环駝t計算[交換主元行]如果交換行[規(guī)格化主元行]，[消元]對循環(huán)1.5.1)1.5.2)[交換列]對循環(huán)如果，交換列[算法結束]以下是計算機求逆驗證矩陣 (3.9) ， (3.10)3.2 LU分解 1.消去法的矩陣運算表示形式 從矩陣計算的角度講Gauss消去法是通過一系列初等行變換將矩陣約化成上三角矩陣；Guass-Jordan消去法是通過一系列初等行變換將矩陣約化成單位矩陣。矩陣的初等行變換本質(zhì)上是矩陣左乘一個初等矩陣，因此消去法也可以用矩陣乘積形式進行描述。 Gauss消去法的矩陣形式： Gauss-Jordan消去法的矩陣形式： 其中(3.9)式和(3.10)式中的就是我們前述消去法的消元因子，如果記： (3.11)則Gauss消去法的整個過程可以表示為：注意到形如(3.11)式的單位下三角矩陣具有以下兩個性質(zhì)：逆陣仍然是單位下三角陣；任意兩個單位下三角陣的乘積仍然是單位下三角陣，所以上式也可以表示為：因此Gauss消去法的本質(zhì)上是將系數(shù)矩陣A進行如下的矩陣分解： (3.12)其中是單位下三角矩陣，是上三角矩陣，此種分解也稱為杜里特爾(Doolittle)分解，如果是下三角矩陣，是單位上三角矩陣，則稱為克洛特(Crout)分解，以下內(nèi)容討論的是杜里特爾分解方法。LU分解上述(3.12)式也可以寫成： 實現(xiàn)矩陣的LU分解算法既可以仿照消去法設計算法也可以直接從該式推導出計算公式，以下我們討論直接LU分解方法。首先用L的第1行乘以U的各列，則有：，由此可以得到U的第1行：再讓L的各行乘以U的第1列得到，由此又可以得到L的第1列假如我們已經(jīng)得到了L的前k-1列和U的前k-1行：按如下方法確定U的第k行和L的第k列，將L的第k行與U的第各列相乘，有 從該式可以得到U的第k行 將L的第i(i>k)各行與U的第k列相乘，有 從該式可以得到L的第k列 則有如下結果：按如此方法，從執(zhí)行到n-1就完成了矩陣的LU分解。綜上所述，我們得到直接LU分解的計算公式： (3.13)定理3.2設n階矩陣A的順序主子式均不為0，則A存在LU分解，且分解是唯一的。事實上，如果分解成功，A的第k個順序主子式的值就是，因此結論是顯然的。本定理表明LU分解可以正常進行的前提是順序主子式不為0，表現(xiàn)在分解過程中就是所有。 通?？梢詫⑾氯蔷仃嘗和上三角矩陣U放在一個矩陣內(nèi)，即： (3.14)例3.4求矩陣A的LU分解 解：設對于k=1，利用(3.13)式，先計算U的第1行 再根據(jù)的計算公式計算L的第1列 ， ， 當前的分解為 對于k=2，根據(jù)的計算公式計算U的第2行 ， ， 再根據(jù)的計算公式計算L的第2列 ， 當前的分解為對于k=3，根據(jù)的計算公式計算 最終得到分解為：將下三角矩陣L和上三角矩陣U放在一個矩陣內(nèi)，如：用LU分解求解線性方程組考慮線性方程組如果系數(shù)矩陣A存在LU分解A=LU，代入原方程組令，則可先求解然后再求解即可得到原線性方程組的解。此處兩個方程組，一個是下三角方程，一個是上三角方程。上三角方程的求解可以使用算法3.1，現(xiàn)在考慮下三角方程的求解。將下三角方程展開求解上三角線性方程組是按的次序求解所以稱為回代過程，求解此下三角線性方程組時次序正好相反是所以稱為順代過程，計算公式如下： 算法3.6(求解線性方程組的直接LU分解算法)[初始化]輸入系數(shù)矩陣A，右端項b[LU分解]對于循環(huán)[計算U的第k行][計算L的第k列][解下三角方程][解上三角方程][算法結束]例3.5用LU分解方法求解下列線性方程組 解：先對系數(shù)矩陣LU分解得到：解下三角方稱：解上三角方程：前述的LU分解相當于順序Gauss消去法，如同順序Gauss消去法不具實用性一樣，算法3.6也不具實用性，應該使用更穩(wěn)定的列主元策略。由線性代數(shù)知識可知，交換兩行相當于系數(shù)矩陣左乘一個交換矩陣，將所有交換矩陣的乘積記為P，就得到如下定理：定理3.3如果矩陣A非奇異(矩陣可逆)，則存在交換矩陣P，單位下三角陣L和上三角陣U，使得 (3.15) 用列主元實現(xiàn)矩陣的直接LU分解遠不如對算法3.3稍作改進更容易，以下我們給出這一改進的算法。 算法3.7(列主元LU分解算法)[初始化]設置系數(shù)矩陣A，向量存放主元行下標，取充分小的數(shù)；[分解過程]對于循環(huán)[選主元]如果，停止計算，算法失敗！如果交換第行對于循環(huán)1.4.1)1.4.2)[算法結束]本算法的幾點說明：向量p記錄主元行的下標，必要時可以生成交換矩陣P，完成矩陣或方程的交換功能，在求解方程時根據(jù)主元行下標對右端項作相應的交換；交換主元行時，要對n個元素進行交換，這是因為此時矩陣L也要做相應的交換；每次循環(huán)都對k行k列以后的元素進行更新，步驟1.4.1)是對的最后計算結果，對的計算則隱含在步驟1.4.2)中；完成算法之后，下三角矩陣L除主對角上的1之外其他元素置于矩陣A對應的位置上，上三對角矩陣U置于矩陣A的對應位置上。3.3 三對角方程組的追趕法本節(jié)考慮如下形式的線性方程組 (3.16)稱上述方程組為三對角方程組，其系數(shù)矩陣稱為三對角矩陣。其矩陣形式為： (3.17)其中系數(shù)矩陣A是三對角矩陣。在一些實際問題中，例如解常微分方程邊值問題、偏微分方程的差分方程、解熱傳導方程、三次樣條函數(shù)、數(shù)值微分等，都會遇到此類線性方程組，求解此類方程組，用Gauss消去法顯然會浪費內(nèi)存和計算量。以下我們通過LU分解，導出求解三對角方程更為簡便有效的算法，稱其為追趕法。 對三對角矩陣進行求解LU分解：右端兩矩陣相乘我們可以得到：從而得到以下三對角陣LU分解的遞推公式： (3.18)在得到三對角陣的LU分解之后，先解下三角方程：由此得到y(tǒng)的計算公式： 再求解上三角方程：由此得到x的計算公式： 定理3.4設三對角線性方程組中系數(shù)矩陣A滿足：則系數(shù)矩陣A非奇異，且有：證明：注意到結論1)與系數(shù)矩陣A非奇異是等價的。下面我們用歸納法來證明結論1)和結論2)。對于，，顯然結論成立；先假設結論對成立，則對有從而可以得到結論1)和結論2)均成立。考慮結論3)，根據(jù)結論2)和的計算公式有：算法3.8(求解三對角方程的追趕法)[初始化]輸入數(shù)據(jù)，令。[矩陣分解][解下三角方程][解上三角方程][算法結束]例3.6求解三對角方程(計算過程保留小數(shù)點以后3位數(shù)字)解：首先對三對角矩陣進行LU分解，結果為：先解下三角矩陣再解上三角矩陣用追趕法求解三對角方程組的計算量，也就是乘除法運算次數(shù)僅為，所以對于三對角方程組通常不使用Guass消去法而是使用追趕法也就很容易理解了。有時還會遇到如下形式稱為廣義三對角方程組(如第4章的三次樣條中第三類邊界條件形成的三彎矩方程)： (3.19)對(3.19)式也可以設計出計算量的算法，仿照Gauss消去法，可用如下方法求解。 首先消元第一列：顯然僅需要消元，在完成對應的兩個方程的消元之后會增加兩個位置的非零元素，我們?nèi)匀挥帽硎?，成為如下形式：一般情況，對于，每列僅兩個位置的元素應該消元，同時處在的位置上會新增兩個非零元素，我們?nèi)杂洖楸敬伪幌脑胤枺斖瓿闪说南^程時是如下形式 (3.20)顯然只要再完成(3.20)式對的消元就可以將其變?yōu)槿缦滦问? (3.21)剩下的問題就是一個簡單的回代算法。算法3.8A求解廣義三對角方程算法[初始化]輸入數(shù)據(jù)對于執(zhí)行，，，，對于做，，，，對于，[回代過程][算法結束]顯然容易算出該算法的計算量約為次乘除法運算。3.4 平方根法平方根法本節(jié)我們討論對稱正定矩陣的一個更為有效的分解方法—平方根法或Cholesky分解，首先回顧線性代數(shù)中對稱正定矩陣的有關結論。 定義3.1設矩陣滿足，則稱矩陣A是對稱矩陣，如果對任意非零向量均有 則稱A是正定矩陣，如果有則稱A是半正定矩陣。 定理3.5設A是對稱正定矩陣，則有以下結論成立：所有特征值為正；所有順序主子式為正；如果對稱正定矩陣可逆，則逆矩陣也是對稱正定矩陣。定理3.6設，且的所有順序主子式均大于0，則存在唯一的單位下三角陣和對角陣，使得 (3.22)證明：因為的順序主子式不為零，所以存在唯一的LU分解A=LU，注意到上三角矩陣U可以分解為：將其帶入則有，由對稱性得再根據(jù)LU分解的唯一性，則有得證。 注意到如果矩陣A對稱正定的，則有，令：將(3.22)式可以改寫成如下形式：則有如下結論。定理3.7(Cholesky分解)設是正定對稱矩陣，則存在唯一的非奇異下三角陣，使得 (3.23)其中的對角元素均為正數(shù)，稱為Cholesky(喬萊斯基)分解或平方根分解。 至此我們僅給出了Cholesky分解的存在性條件，以下討論分解算法。如同LU分解，也考慮直接的分解方法，從下式推導出計算公式： (3.24)我們寫出下三角部分的表達式為假設已得到L的第1列到第k-1列，現(xiàn)要計算第k列，則有 (3.25) 例3.7求以下對稱矩陣的Cholesky分解，并求其行列式的值 解： 2．改進平方根法 平方根算法是穩(wěn)定的，缺點是算法中有開平方運算。一次開平方運算的CPU時間遠遠大于乘除運算(開平方運算一般是迭代法實現(xiàn)的)，因此通常我們應該盡量避免使用平方根運算，改進平方根法就是無開方運算的平方根法。 其實，所謂無開方運算的平方根法就是基于上節(jié)所述的(3.22)式進行求解，將(3.22)展開以作為中間矩陣，則有，直接對按矩陣乘法規(guī)則推導出計算公式為： (3.26)如果求解線性方程組Ax=b，以下只要分別求解下三角方程和上三角方程，即用如下計算公式： (3.27)在編程實現(xiàn)改進平方根算法時，對角陣可以存放在矩陣的對角線上，中間矩陣可以借用矩陣的上三角部分，改進平方根算法如下： 算法3.9改進平方根矩陣分解算法[初始化]輸入系數(shù)矩陣和右端項[分解]對于循環(huán)[計算L]對[計算D][算法結束]注：(1)在編程時，如果希望將下三角矩陣L和中間矩陣T存放在原矩陣A內(nèi)，只需要將算法3.9中的所有改為即可，其它不須作任何修改；(2)如果利用改進平方根算法求解線性方程組，可以在算法3.9的基礎上按照(3.27)公式繼續(xù)計算即可以得到方程組的解。例3.8用改進平方根算法分解以下矩陣解：按算法3.9步驟1)可求得如下結果：即：迭代法以上所述Guass消去法、LU分解求解線性方程組的方法均屬于直接方法，還有一類求解線性方程組的方法屬于迭代法，迭代法在求解大型稀疏問題、病態(tài)方程組等特殊情況時是一個非常有效的方法。求解線性方程組的直接方法和迭代法的主要區(qū)別是：直接法可以用有限次初等運算求出其解，比如Gauss消去法可以用大約次乘除法運算和次加減法運算求出方程組的解；而迭代法思想如同非線性方程的迭代法一樣，先給出一個初始近似解，根據(jù)一個迭代公式得到一個更好的近似解，通過反復執(zhí)行迭代過程逐步逼近到精確解，當近似解達到指定的精度時結束計算。如同求解非線性方程的迭代法，求解線性方程組迭代法的也是先把所給的線性方程組，化成同解的線性方程組 (3.28)這里稱B為迭代矩陣，給定初始向量并按迭代格式 (3.29)進行迭代計算。如果按上述迭代格式所得到的向量序列收斂于某一向量，則就是方程組的解，并稱此迭代法收斂；否則稱為不收斂或發(fā)散。本課程僅介紹Jacobi迭代法和Gauss-Seidel迭代法，并簡單討論其收斂條件。Jacobi迭代法設有線性方程組寫成其緊湊形式為 假設系數(shù)矩陣非奇異且主對角線元素，將方程組的非對角線項移至等式的右端并兩端除以對角線元素，則方程組的等價形式為： (3.30)按此式建立迭代格式： (3.31)則稱為Jacobi迭代格式，其相應的迭代方法稱為Jacobi迭代法。為便于今后的收斂性分析，現(xiàn)將迭代公式改寫成矩陣形式。記則，方程組改寫成：則有：相應的矩陣形式的迭代公式為：簡記為： (3.32)其中稱為Jacobi迭代矩陣，。在非線性方程的迭代法中，迭代過程的終止準則是兩相鄰迭代近似解之差的絕對值小于給定精度，而對于線性方程組的迭代法，其終止準則是兩次相鄰迭代近似解之間的某種范數(shù)小于給定精度，即 例3.9用Jacobi迭代法求解方程組 解：等價形式的方程組為：對應的迭代格式為：取初始點，終止精度為，迭代過程如下表：序號00.0000.0000.00010.444-0.1251.00020.417-0.8610.90930.918-0.7971.12740.851-1.0590.96651.043-0.9411.05960.954-1.0480.97171.035-0.9701.02780.977-1.0260.98191.019-0.9831.014100.987-1.0130.990111.010-0.9901.008120.993-1.0070.994131.005-0.9951.004140.996-1.0040.997151.003-0.9971.002160.998-1.0020.998171.002-0.9981.001180.999-1.0010.999191.001-0.9991.001200.999-1.0010.999211.000-1.0001.000經(jīng)21次迭代，收斂到線性方程組的精確解，收斂速度較慢。 算法3.10求解線性方程組的Jacobi迭代法[初始化]置初始近似解，最大迭代次數(shù)，精度；[迭代]當時循環(huán)執(zhí)行如下步驟對于循環(huán)如果，則停止計算，取為方程組的近似解；令，繼續(xù)執(zhí)行下一次循環(huán)。[算法結束]Guass-Seidel迭代法前述的Jacobi迭代的優(yōu)點是算法簡單，缺點是收斂速度較慢。我們稍加觀察就會注意到，如果迭代是收斂的，則有理由認為一個新的近似解應該優(yōu)于前一個近似解，但在迭代過程中并不是充分利用了最新的計算結果。比如在計算時仍然在使用，其實此時已經(jīng)計算出來，而近似值應該比更好，如果用參與新的計算結果會更好些。Guass-Seidel迭代法就是基于此思想而設計的，因此我們將(3.31)式改寫為： (3.33)稱其為Guass-Seidel迭代格式，相應的迭代法稱為Guass-Seidel迭代法。以下我們給出Gauss-Seidel迭代格式的矩陣形式，沿用前述的矩陣記號，(3.33)式可以寫成等價的形式：迭代格式的矩陣形式為： (3.34)這里稱為Gauss-Seidel迭代矩陣，。 例3.10請用Gauss-Seidel迭代法求解例3.9中的線性方程組，初始點和終止精度仍按原例所取。 解：首先寫出Gauss-Seidel迭代格式為：如同例3.9，取初始點，終止精度為，迭代過程如下表：序號00.0000.0000.00010.444-0.2360.94020.497-0.8371.09730.881-1.0311.04341.016-1.0311.00451.020-1.0080.99661.006-0.9990.99871.000-0.9991.00080.999-1.0001.000顯然僅經(jīng)過8次迭代，已經(jīng)達到了終止精度，迭代速度遠勝過Jacobi迭代法。 算法3.11(求解線性方程組的Gauss-Seidel迭代法)[初始化]置初始近似解，最大迭代次數(shù)，精度；[迭代]當時循環(huán)執(zhí)行如下步驟對于循環(huán)如果，則停止計算，取為方程組的近似解；令，繼續(xù)執(zhí)行下一次循環(huán)；[算法結束]讀者如果編程實現(xiàn)Jacobi迭代法和Gauss-Seidel迭代法會發(fā)現(xiàn)，Gauss-Seidel迭代法更易編程且使用更少的存儲。一般認為，Gauss-Seidel迭代法的收斂速度要優(yōu)于Jacobi迭代法，數(shù)值實驗表明對絕大多數(shù)線性方程組，確實如此，但并無法保證對每一個例子均是如此，讀者可以把下列兩個矩陣作為系數(shù)矩陣構成線性方程組，隨便取非0的右端項，前述的兩種迭代法，各有一個收斂而另一個不收斂。， 超松弛(SOR)迭代解線性方程超松馳迭代法是目前解大型方程組的一種最常用方法，該算法是對Gauss-Seidel迭代法的一種加速方法，其思想是根據(jù)兩次相鄰迭代的近似解構造一個新的更好的近似解。在Gauss-Seidel迭代過程中，對Gauss-Seidel迭代法進行如下修正： (3.35)該迭代法稱為逐次超松弛(SOR)迭代法，其中系數(shù)稱為松弛因子，當時稱為低松弛迭代，取為就是Gauss-Seidel迭代，當稱為超松弛迭代。 例3.11用超松弛迭代法求解方程組 解：取初始點，精度均取。松弛因子的迭代結果序號01.0001.0001.00015.2503.813-5.04723.1413.883-5.02933.0883.927-5.01843.0553.954-5.01153.0343.971-5.00763.0213.982-5.00473.0133.989-5.00383.0083.993-5.00293.0053.996-5.001103.0033.997-5.001113.0023.998-5.000123.0013.999-5.000松弛因子的迭代結果序號01.0001.0001.00016.3133.520-6.65022.6223.959-4.60033.1334.010-5.09742.9574.007-4.97353.0044.003-5.00662.9964.001-4.99873.0004.000-5.000 算法3.12求解線性方程組的逐次超松弛(SOR)迭代法[初始化]置初始近似解、松弛因子、迭代終止，最大迭代次數(shù)N，精度，；[迭代]當時循環(huán)執(zhí)行如下步驟對于循環(huán)如果，則停止計算，取為方程組的近似解；令，繼續(xù)執(zhí)行下一次循環(huán)[算法結束]迭代法的收斂性我們使用任何迭代算法，都要清楚該算法必須具有收斂性、了解收斂速度如何？使用發(fā)散的迭代算法或者是收斂速度太慢的迭代算法是沒有實際意義的，因此我們要研究一下前述的Jacobi迭代法和Gauss-Seidel迭代法的收斂性、收斂速度問題。 定義3.2設是上的任一向量范數(shù)。若存在使則稱序列收斂于。無論在任何一種范數(shù)下，序列收斂于，根據(jù)向量范數(shù)的等價性，必可推出在其它范數(shù)下也收斂于，即收斂性與范數(shù)的選擇無關。 定理3.8(迭代法收斂的充分條件)對于迭代形式的線性方程組，如果則對于任意的初始向量及任意的右端向量，解此線性程組的迭代格式是收斂的，且 (3.36) (3.37) 證明：由于，根據(jù)定理1.5矩陣非奇異，線性方程組有唯一解，即有，與迭代公式相減兩邊取范數(shù)，并利用矩陣范數(shù)與矩陣范數(shù)的相容性則有反復使用該不等式，則有注意到對上式兩端取極限也就是迭代格式收斂。同樣根據(jù)前邊分析還有從而有(3.36)式成立，將(3.36)式再綜合下式可推得(3.37)式成立，即 本定理說明迭代法的收斂性與右端項無關，也與所取初始點無關，但初始點的選取與收斂的快慢有關。例3.12設線性方程組的系數(shù)矩陣為試判斷使用Jacobi迭代法和Gauss-Seidel迭代法求解該方程組的收斂性。 解：分解系數(shù)矩陣為對于Jacobi迭代法，有對迭代矩陣取1-范數(shù)，則有，根據(jù)定理3.8可知Jacobi迭代法收斂。對于Gauss-Seidel迭代法，有對迭代矩陣取-范數(shù)，則有，同樣是根據(jù)定理3.8可得Gauss-Seidel迭代法收斂。 定理3.9(迭代法收斂的充要條件)對于迭代形式的線性方程組，迭代法收斂的充要條件是譜半徑小于1，即。 證明：如果方程組的解為，則應有顯然的充要條件是，根據(jù)定理1.7可知，極限成立的充要條件是，證畢。 與定理3.8相比，定理3.9的條件不容易驗證，以下我們再給出一個容易驗證的收斂條件，首先給出如下定義。 定義3.3如果矩陣滿足則稱矩陣A為對角占優(yōu)矩陣，如果嚴格不等式成立，則稱A是嚴格對角占優(yōu)的。 定理3.10在求解線性方程組的迭代法中：如果矩陣A嚴格對角占優(yōu)，則Jaocbi迭代法和Gauss-Seidel迭代法均是收斂的，當松弛因子滿足時，SOR迭代法是收斂的；如果矩陣A是對稱正定的，則Gauss-Seidel迭代法是收斂的，當松弛因子時，SOR迭代法是收斂的。在以上給出的判斷迭代法收斂條件中，主要有矩陣對角占優(yōu)、矩陣范數(shù)小于1、譜半徑小于1、矩陣對稱正定這四個條件，顯然前兩個判斷方法簡單，而后兩個因為要計算特征值或順序主子式則不便于使用。3.6 線性方程組的誤差分析 在本章的前幾節(jié)，我們給出了一些線性方程組的求解方法，一個還沒有涉及到而且很重要的問題是，這些方法所得到的解可靠嗎？如果是近似解其近似程度如何？在何種情況下所求的解是可靠的？何種方法求出的解更為可靠？這就是本節(jié)所要研究的問題。病態(tài)線性方程組在求解出方程組的計算解后，直觀上會認為可用如下方法度量其計算解的準確性，計算其解的殘向量：如果的值很小，可能認為比較精確。但有時并非如此，看下述線性方程組： (3.38)此方程組的精確解是，但將帶入方程組計算殘量的-范數(shù)值為此時殘向量確實已經(jīng)很小，但是無論如何也不可以看作是的近似解，其實是下列方程組的精確解而這兩個線性方程組的差別僅僅是元素相差，這種差別可以當做有一個微小的擾動。定義3.4如果線性方程組的系數(shù)矩陣或右端項有一個微小的擾動將導致擾動前后兩方程組的解有一個很大的變化，我們稱此線性方程組是病態(tài)的，其系數(shù)矩陣稱為病態(tài)矩陣；反之稱此方程組是良態(tài)的。眾所周知，計算機求解線性方程組的每步計算都會產(chǎn)生舍入誤差，而當前的舍入誤差又會傳播到之后的求解過程中去，導致誤差不斷的積累，而這些計算過程中的誤差及誤差積累勢必影響到最終計算解誤差的大小。目前對這類誤差的分析采用一種稱之為向后誤差分析的方法，其思想是在誤差分析過程中，把舍入誤差導致最終解的誤差看作是初始數(shù)據(jù)有一個擾動，把有舍入誤差的計算過程看作是精確的運算。之后的討論就是基于如此思想的誤差分析方法，為了度量方程組的病態(tài)或良性程度，首先給出關于矩陣條件數(shù)的概念。矩陣的條件數(shù)定義3.5設是非奇異矩陣，對于某一矩陣范數(shù)，稱 (3.39)為矩陣A求解線性方程組的條件數(shù)。特別如果取矩陣的2-范數(shù)，則可記作：如果矩陣A是對稱矩陣，分別是A的最大、最小特征值，則有 考慮(3.38)式的方程組，系數(shù)矩陣為，而的兩個特征值分別為，，則其系數(shù)矩陣A在2-范數(shù)下的條件數(shù)為，顯然該條件數(shù)相比較于系數(shù)矩陣的元素已經(jīng)非常大了。 通常我們使用(3.39)式定義的條件數(shù)來度量線性方程組的病態(tài)程度，其中的矩陣范數(shù)可取為2-范數(shù)、1-范數(shù)、范數(shù)中的任一范數(shù)，條件數(shù)愈大方程組的病態(tài)程度就愈嚴重。 矩陣的條件數(shù)具有如下性質(zhì)：；；；。容易驗證單位矩陣的條件數(shù)為1，且無論是2-范數(shù)、1-范數(shù)、范數(shù)中的任一范數(shù)，單位矩陣的條件數(shù)均達到最小值。對于任取的非奇異矩陣A，有 特別需要指出的，2-范數(shù)下正交矩陣Q其條件數(shù)也為1，即也達到條件數(shù)的最小值。線性方程組的誤差分析如前所述，可以把求解過程的誤差看作是原始數(shù)據(jù)的擾動問題，首先考慮系數(shù)矩陣的擾動對方程組解的影響。設是方程組的精確解，是方程組的精確解，即假設是系數(shù)矩陣有一個擾動之后的方程組的精確解：展開并化簡得到則，從而有，得最后寫成如下形式 (3.40)該式給出了一個系數(shù)矩陣的擾動對最終所求解的誤差界，顯然影響解的精度的主要因素是系數(shù)矩陣的條件數(shù)，從該式可以看出，系數(shù)矩陣的一個微小的擾動，在計算解中被放大了倍，其實條件數(shù)就是誤差的放大因子。 再看右端項的擾動，設是方程組右端項有一個擾動后方程組的精確解，即則有，從而有，再由，以下兩不等式相除可得 (3.41)對于右端項的擾動，其相對誤差也被放大了倍。 最后看系數(shù)矩陣與右端項均有擾動的情況。設是系數(shù)矩陣有一個擾動同時右端項也有一個擾動之后的方程組的解，即則有，從而，當擾動較小時，可以在分析中忽略擾動的二次項，即忽略上式中的項，兩邊取范數(shù)有從而有 (3.42)根據(jù)向后誤差分析的觀點，在實際計算過程中產(chǎn)生的誤差及誤差的傳播、積累可以看作是原始數(shù)據(jù)的兩個擾動導致的誤差，從(3.42)式可以看出這些誤差被放大了cond(A)倍。因此我們可以用條件數(shù)cond(A)來度量求解線性方程組的解的病態(tài)程度，條件數(shù)愈大愈病態(tài)，條件數(shù)愈小愈良性。 現(xiàn)在我們看一個嚴重病態(tài)的矩陣—Hilbert矩陣 (3.43)其逆矩陣為 (3.44)此矩陣的條件數(shù)見下表：Hilbert矩陣的條件數(shù)表n581012此表說明，以Hilbert矩陣為系數(shù)矩陣的線性方程組，將是嚴重病態(tài)的，當方程組的階數(shù)大于15時，不論是用列主元Gauss消去法還是用正交分解方法求解所得解已經(jīng)面目全非，讀者可以編程一試。病態(tài)線性方程組的求解線性方程組是否為病態(tài)是由其系數(shù)矩陣的條件數(shù)確定的，但是求矩陣的條件數(shù)將是十分困難的，通常并不會去計算矩陣的條件數(shù)，而是根據(jù)計算過程中的某些現(xiàn)象來判斷方程組是否為病態(tài)，比如消元過程出現(xiàn)小的主元，則有可能是病態(tài)的；接近奇異的系數(shù)矩陣基本肯定是病態(tài)的等等。關于病態(tài)方程組的求解方法，有以下幾種可選擇的方法：最理想的方法是提高計算機的字長，增加基本運算精度，用雙精度、四倍精度甚至更高的精度進行計算，這些已經(jīng)超出本課程討論的范圍;改善系數(shù)矩陣的條件數(shù)、即尋找預處理矩陣P，對矩陣做預處理： (3.45)這里，使得方程組變成良態(tài)的，但預處理矩陣的求解是非常困難的；迭代改善、迭代改善是最為經(jīng)典的求解病態(tài)方程組的算法，其思想是用單精度求得方程組的解，之后用雙精度計算其殘量，然后再依此殘量為右端項用單精度求解新構成的方程