2022電力系統(tǒng)分析中的計算方法

電力系統(tǒng)分析中的計算方法（原書第2版）目 錄第1章 引

……………………1第2章 線性方程組的解

…………………32.1Gauss消去法 42.2LU分解法 82.2.1采用部分選主元的LU分解法 152.2.2采用完全選主元的LU分解法 192.3條件數(shù)與誤差傳播 202.4松弛法 212.5共軛梯度法 266廣義最小殘差法 302.7問題 35第3章 非線性方程組的解法 391不動點迭代法 393.2Newton-Raphson迭代法 453.2.1收斂性 473.2.2用于求解非線性方程組的Newton-Raphson法 483.2.3Newton-Raphson法的改進 513.33.4

……………………52……………………553.5數(shù)值微分法 573.6在電力系統(tǒng)中的應(yīng)用 603.6.1潮流計算 613.6.2調(diào)節(jié)變壓器 683.6.3解耦潮流算法 723.6.4快速解耦潮流算法 743.6.5PV曲線與連續(xù)潮流計算 77目 錄 Ⅶ6.6三相潮流計算 823.7問題 83第4章 稀疏矩陣求解技術(shù) 861存儲方法 864.2稀疏矩陣的表示方法 904.3排序方案4.3.1方案0

…………………92………………984.3.2方案Ⅰ 994.3.3方案Ⅱ 1044.3.4其他方案 1064在電力系統(tǒng)中的應(yīng)用 1074.5問題 110第5章 數(shù)值積分 1145.1單步法 1151.1基于Taylor級數(shù)的算法 1155.1.2向前Euler法 1155.1.3Runge-Kutta法 1165.2多步法 1165.2.1Adams算法 1225.2.2Gear法 1255.3精度與誤差分析 1265.4數(shù)值穩(wěn)定性分析 1295.5剛性系統(tǒng) 1335.6步長選擇 1375.7微分代數(shù)方程 1405.8在電力系統(tǒng)中的應(yīng)用 1415.8.1暫態(tài)穩(wěn)定性分析 1428.2中期穩(wěn)定性分析 1495.9問題 152第6章 優(yōu)化方法 1571最小二乘狀態(tài)估計 1576.1.1加權(quán)最小二乘狀態(tài)估計 1616.1.2壞數(shù)據(jù)的檢測 1636.1.3非線性最小二乘狀態(tài)估計 1666.2線性規(guī)劃 167Ⅷ 電力系統(tǒng)分析中的計算方法（2版）6.2.1單純形法 1686.2.2內(nèi)點法 1716.3非線性規(guī)劃 1766.3.1二次規(guī)劃 1766.3.2最速下降法 1786.3.3序貫二次規(guī)劃法 1824在電力系統(tǒng)中的應(yīng)用 1846.4.1最優(yōu)潮流 1846.4.2狀態(tài)估計 1946.5問題 198第7章 特征值問題 2027.1冪法 2027.2QR法 2057.2.1移位QR法 2117.2.2縮減法 2123Arnoldi法 2127.4奇異值分解 2197.5模態(tài)辨識 2227.5.1Prony法 2237.5.2矩陣束法 2257.5.3Levenberg-Marquardt法 2267.5.4Hilbert變換 2287.5.5舉例 2296在電力系統(tǒng)中的應(yīng)用 2357.7問題 236第1章 引 論在今天放松管制的大環(huán)境下，本國電網(wǎng)已不能按照當(dāng)初設(shè)計時設(shè)想好的方式運行。因此，系統(tǒng)分析和計算對于預(yù)測并不斷地更新電網(wǎng)的運行狀態(tài)就顯得非常重要。系統(tǒng)分析和計算的內(nèi)容包括估算當(dāng)前的潮流和電壓（潮流分析和狀態(tài)估計），確定系統(tǒng)的穩(wěn)定極限（連續(xù)潮流、用于暫態(tài)穩(wěn)定計算的數(shù)值積分和特征值分析），以及費用最小化計算（最優(yōu)潮流）等。本書講述各種計算方法的入門知識，這些計算方法構(gòu)成了電力系統(tǒng)及其他工程和科學(xué)領(lǐng)域很多分析計算的基礎(chǔ)，也是無數(shù)商業(yè)軟件所使用算法的基礎(chǔ)。通過了解支撐這些算法的理論，讀者可以更好地使用那些商業(yè)軟件并做出更明智的決策（即在仿真軟件中選擇更好的數(shù)值積分方法和確定更合適的積分步長等）。由于電網(wǎng)的規(guī)模巨大，電力系統(tǒng)分析計算依靠手工幾乎是不可能的，計算機是唯一真正可行的手段。電力工業(yè)界是計算機技術(shù)的最大用戶之一，也是在大型計算機出現(xiàn)后最早接受計算機分析的工業(yè)部門之一。雖然電力系統(tǒng)分析的第一批算法是在20世紀(jì)40年代開發(fā)出來的，但直到60年代計算機才在電力工業(yè)界得到廣泛使用。今天用于大系統(tǒng)仿真和分析的很多技術(shù)和算法起初是為電力系統(tǒng)應(yīng)用而開發(fā)的。隨著電力系統(tǒng)更多地運行在重載狀態(tài)下，計算機仿真將會在電力系統(tǒng)的控制和安全性評估方面發(fā)揮更大的作用。將商業(yè)化軟件包用于仿真重載的電力系統(tǒng)時，往往不能運行或者給出錯誤的結(jié)果。為了正確理解和解釋商業(yè)化軟件包的仿真結(jié)果，必須了解其內(nèi)在的數(shù)值計算方法。例如，是系統(tǒng)確實呈現(xiàn)出了仿真得到的行為特征，還是僅僅因為數(shù)值計算不精確而導(dǎo)致的偽結(jié)果？如何彌補商業(yè)化軟件包存在的算法缺陷？是通過選擇更好的仿真參數(shù)，還是通過將問題描述成更適合于數(shù)值計算的形式？對諸如此類的問題，受過數(shù)值計算訓(xùn)練的用戶可以做出更好的判斷。本書將給出多種廣泛使用的數(shù)值計算方法的背景知識，這些計算方法構(gòu)成了眾多電力系統(tǒng)分析和設(shè)計商業(yè)化軟件包的基礎(chǔ)。本書是為研究生計算方法課程而寫的教材，課程長度是一個學(xué)期。盡管本書中的大多數(shù)例子是基于電力系統(tǒng)應(yīng)用的，但數(shù)值計算的理論是按照一般性的方式闡述的，因而可應(yīng)用于大范圍的工程系統(tǒng)。雖然完全領(lǐng)會本書中某些內(nèi)容的微妙之處需要具有電力系統(tǒng)工程的一些知識，但這些知識對理解算法本身并不是前提條件。相對于內(nèi)容廣泛的計算方法領(lǐng)域，本書的敘述和例子僅僅是一個入門的導(dǎo)引，并不能包羅萬象。本書闡述的很多算法一直存在著眾多的改進算法，并且仍2 2 電力系統(tǒng)分析中的計算方法（2版）然是當(dāng)前研究的目標(biāo)；由于本書的意圖是打基礎(chǔ)，因此很多新的進展并沒有明確地包含在敘述中，而是作為參考文獻列出，以方便感興趣的讀者進一步學(xué)習(xí)。為了讓讀者能夠容易地將書中的算例復(fù)現(xiàn)出來，本書特意挑選了簡單而足夠完整的問題作為算例。很多“真實世界”的問題其規(guī)模和范圍要大得多，然而本書給出的處理問題的方法，應(yīng)該可以為讀者提供足夠的工具，以應(yīng)對其所遇到的困難。本書中的大部分算例是通過Matlab編程實現(xiàn)的，雖然這是本書作者所采用的平臺，但實際上，任何計算機語言都可用來實現(xiàn)本書中的算法，沒有理由偏好某個特定的平臺或語言。第2章 線性方程組的解法在很多工程和科學(xué)領(lǐng)域，希望根據(jù)一組物理學(xué)上的關(guān)系式，能夠用數(shù)學(xué)的方法來確定系統(tǒng)的狀態(tài)。這些物理學(xué)上的關(guān)系式可以是描述諸如電路拓撲、質(zhì)量、重量、力等物理特性的方程，例如，電路系統(tǒng)中針對每個節(jié)點的由注入電流、網(wǎng)絡(luò)拓撲和支路阻抗構(gòu)成的約束方程。在很多情況下，描述已知量（即輸入）與未知狀態(tài)（即輸出）之間的關(guān)系式是線性的。因此，一般性地可以將線性系統(tǒng)描述為Ax=b （2.1）式中，b是一個n1向量，為已知量；x是一個n1的未知狀態(tài)向量；而A是將x與b聯(lián)系起來的nn階矩陣。暫時我們假定矩陣A是可逆的，即非奇異的，這樣，每個向量b將對應(yīng)一個唯一的向量x?。也就是矩陣A-1存在，且x?=A-1b （2.2）是式（21）的唯一解。對式（2.1）進行求解的自然方法是直接計算矩陣A的逆，并將它與向量b相乘。計算A-1的方法之一是采用Cramer法則：（A1（i，j）det1A）（Aij）T 對于i=1，…，n，j=1，…，n（2.3）式中，A1（i，j）A1的第i行、j列元素，Aij是矩陣A的元素aij的余子式。這種方法需要計算（n1）個行列式，從而需要2（n1）！次乘法運算。當(dāng)n值較大時，計算量爆炸式增長，不可控制。因此，必須開發(fā)替代的算法。（一般地說，求解式（21）有兩種基本的方法：直接法，也稱為消去法。該方法通過有限次算術(shù)運算找到上述方程的精確解（指在計算機的精度內(nèi)）。如果不計計算機的舍入誤差，采用直接法得到的解x將是完全精確的。迭代法。這種方法的每次計算都基于相同的計算過程，試圖產(chǎn)生一系列不斷逼近真解的近似解。當(dāng)所獲得的近似解達到預(yù)先設(shè)定的精度指標(biāo)或者已確定迭代過程不再對精度有改進時，終止迭代。求解方法的選擇通常取決于所研究系統(tǒng)的結(jié)構(gòu)，某些系統(tǒng)更適合于采用特定的方法進行求解。一般地，直接法更適合于滿矩陣，而迭代法更適合于稀疏矩陣。但是，正像對于大多數(shù)一般性結(jié)論一樣，上述經(jīng)驗法則存在很多的例外。PAGEPAGE4 電力系統(tǒng)分析中的計算方法（2版）第第2章 線性方程組的解法 PAGE1121 Gauss消去法Gauss消去法是一種不需要顯式計算A-1而對式（2.1）進行求解的方法。因為x是直接求出的，因此Gauss消去法是一種求解線性方程組的直接法，并且是一種常用的直接法。Gauss消去法的基本思路是利用第1個方程消去其余方程中的第1個未知量，依次重復(fù)這個過程就能消去第2個未知量、第3個未知量……直到整個消去過程完成，此時，第n個未知量可以直接從輸入向量b中求出。然后，對上述方程進行遞歸回代就能求出所有的未知量。Gauss消去法是這樣一個過程，它通過一系列初等行變換運算，將n（n+1）階的增廣矩陣[A|b]變換為n（n+1）階的矩陣 其中，

I|b?]Ax=bA-1Ax=A-1bx=A-1b=b?x?=b?這樣，如果存在一系列初等行變換運算可以將矩陣A變換成一個單位矩陣I，那么，對向量b進行相同系列的初等行變換運算就能將向量b變換為方程組的解x?。初等行變換運算包括如下3種矩陣運算：1）交換矩陣中的任意兩行。用一個常數(shù)乘任意一行。將任意行的線性組合加到另一行上。選擇一系列初等行變換運算將矩陣A變換成一個上三角矩陣，該上三角矩陣的對角元為1，而所有下三角元素為零。這個過程被稱為“向前消去過程”。向前消去過程的每一步可以通過對矩陣A乘以一個初等矩陣ξ來獲得，這里的初等矩陣ξ指的是通過對單位矩陣進行初等行變換運算就能得到的矩陣。例2.1尋找一系列初等矩陣將如下矩陣變換成上三角矩陣。〓1 3 4 8〓A〓2 1 2 3A〓4 3 5 8〓〓〓〓9 2 7 4〓〓〓解2.1 為了使上述矩陣上三角化，所采用的行變換運算應(yīng)能使此矩陣對角線以下每一列的元素為零。這可以通過如下的行變換運算來實現(xiàn)：將對角線以下的每一行用其本身減去一個常數(shù)乘對角線行來代替，該常數(shù)的選擇應(yīng)使此列中對角線以下的元素為零。因此，A的第2行應(yīng)該用（行2-2行1）來代替，而對應(yīng)此初等行變換運算的初等矩陣是〓 1 0 0 0〓ξ1=〓-2 1 0 0〓〓 0 0 1 0〓〓 0 0 0 1〓和〓1 3 4 8〓ξA1〓0 -5 -ξA1〓4 3 5 8〓〓〓〓9 2 7 4〓〓〓注意除了第2行其他的行沒有變化，而第2行在第1個對角元下的元素變?yōu)榱肆?。類似地，完成消去?列任務(wù)的另外2個初等矩陣是〓ξ2〓

1 0 0 0〓0 1 0 0〓〓-4 0 1 0〓〓ξ3〓〓

0 0 0 1〓1 0 0 0〓0 1 0 0〓0 0 1 0〓〓-9 0 0 1〓〓 〓并且有〓1 3 4 8ξξξA3ξξξA3〓0 -9 -11 -24〓〓〓0 -25 -29 -68〓〓

（2.4）現(xiàn)在向前消去過程將針對第2列展開，首先是將第2個對角元下面的元素全部消去，然后將第2個對角元變換為1。具體過程是〓1 0 0 0〓〓0 1 0 0〓〓ξ4=〓0 -9 1 0〓〓 5 〓〓〓〓0 0 0 1〓〓〓〓1 0 0 0〓〓0 1 0 0〓ξ5=〓0 0 1 0〓〓 25 〓〓0 -5 0 1〓〓〓1 〓0 -

〓10 0〓10 0〓ξ6=〓 5 〓〓0 0 1 0〓〓0 0 0 1〓并且有

〓 〓〓1 3 4 8〓〓0 1 6 13〓〓ξ6ξ5ξ4ξ3ξ2ξ1A〓〓

5 51 3

（2.5）〓0 0 1 -3〓0 0 -5 -〓0 0 1 -3繼續(xù)向前消去過程，有

〓 〓〓1 0 0 0〓0 1 0 0ξ7=〓 〓從而有

〓0 0 1 0〓0 0 5 1〓〓1 0 0 0〓ξ8ξ8〓0 0 -5 0〓0 0 0 1〓〓1 3 4 8〓〓0 1 6 13〓ξ8ξ7ξ6ξ5ξ4ξ3ξ2ξ1A=〓 5 5

（2.6）〓0 0 1 3〓〓0 0 0 -6〓最后有

〓 〓〓1 0 0 0〓〓0 1 0 0〓〓〓ξ9=〓0 0 1 0〓〓〓1〓〓〓0 0 0 -6〓〓〓和〓1 3 4 8〓〓0 1 6 13〓ξ9ξ8ξ7ξ6ξ5ξ4ξ3ξ2ξ1A=〓 5 5

（2.7）〓0 0 1 3〓〓0 0 0 1〓〓 〓從而完成了向前消去過程，將矩陣變換成了上三角矩陣。 ?一旦將矩陣變換成了上三角矩陣，式（2.1）的解向量x的逐次代換（也稱為回代”）而求得。例2.2 利用例2.1的上三角矩陣，求如下方程的解。

就可以通過狀態(tài)量〓1 3 4 8〓〓

〓1〓2 1 2 3

〓x2〓 1〓 〓〓 〓=〓〓〓4 3 5 8〓〓x3〓 〓1〓〓9 2 7 4〓〓x4〓 〓1〓解2.2 注意一系列下三角矩陣的乘積仍然是下三角矩陣。因此，矩陣乘積W=ξ9ξ8ξ7ξ6ξ5ξ4ξ3ξ2ξ1 （2.8）是一個下三角矩陣。將W與矩陣A相乘將得到一個上三角矩陣，WA=U （2.9）即U是一個通過向前消去過程而得到的上三角矩陣。在式（2.1）的兩邊同時左乘W，可以得到 WAx=Wb （2.10）Ux=Wb （2.11）=b′ （2.12）式中Wbb′。21

1 0 0 0〓2 1 〓〓W〓〓5 -5 〓W〓2 9 5 0〓1--11--1〓且

146〓〓〓1

56〓1〓1

〓6〓b′=W〓

〓=〓5〓〓〓 〓6〓〓1〓 〓 〓〓〓 〓3〓〓1〓

〓2〓這樣，

〓1 3 4 8〓〓x1〓〓

〓1〓〓52〓52

13〓〓x5

〓1〓〓5〓〓 〓

〓〓6

（2.13）〓0 0 1 3〓〓3

〓〓〓 〓〓〓

〓3〓〓0 0 0 1通過觀察，x4=3/2。由第3行可得

〓2〓x3=6-3x4 （2.14）將x4的值代入到式（214）中，得到x33/2。類似地，x2=1-6x3-13x4 （2.15）5 5 5將x3和x4的值代入到式（215）中，得到x2=-11/2。以類似的方法求解x12x1=1-3x2-4x3-8x4 （2.16）2從而得到

=-1〓〓〓x1〓 〓-1〓〓

（2.17）〓x2〓=1〓-11〓〓x3〓 2〓 3〓〓 〓 〓 〓〓x4〓 〓 3〓逐次將已求得的x值代回到方程中進行求解的步驟，在Gauss消去法的求解過程中被稱為“回代”。因此，Gauss消去法包括兩個主要步驟：“向前消去”過程和“回代”過程。向前消去過程是將矩陣A變換成其三角因子，而回代過程則根據(jù)輸入向量b和A的矩陣因子求解未知向量x。Gauss消去法同時也為LU因子分解過程提供了框架。2 LU分解法Gauss消去法的向前消去過程產(chǎn)生了一系列與矩陣A相關(guān)的上三角矩陣和下三角矩陣，如式（2.9）所示。其中，矩陣W是一個下三角矩陣，而矩陣U是一個對角元為1的上三角矩陣。考慮到下三角矩陣的逆仍然是一個下三角矩陣，因此如果定義L=?W-1那么A=LU正是基于矩陣L和矩陣U，將矩陣的因子分解或消去算法命名為矩陣的LU分解”。事實上，對于任何非奇異的矩陣A，存在一些置換矩陣P（有可能PI），使得LU=PA （2.18）式中，U是對角元為1的上三角矩陣，L是對角元為非零的下三角矩陣，而P是由單位矩陣I通過行交換或列交換而形成的元素僅為0或1的矩陣。一旦選定一個合適的矩陣P，那么上述的LU分解就是唯一的[6]。如果P、L和U已經(jīng)確定，那么方程組Ax=b （2.19）就可以迅速得到求解。在式（219）兩邊同時左乘矩陣P得PAx=Pb=b′ （2.20）LUx=b′ （2.21）式中，b′僅僅是對向量b的一個重新排列。如果引入一個“啞”向量y，使得Ux=y （2.22）那么考慮到式（2.23）的結(jié)構(gòu)為

Ly=b′ （2.23）〓l11 0 0 … 0〓〓〓〓l21 l22 0 … 0〓〓y2

〓b′1〓〓b′2〓〓〓〓〓〓〓l31 l32 l33 … 0〓〓y3〓=〓b′3〓〓〓〓〓〓? ? ? ? ?〓〓?〓 〓?〓〓ln1 ln2 ln3 … lnn〓

n

〓

n〓向量y的元素可以通過直接代換得到：y1=b′1ly1=y2=y3=

l1（b′2-l21y1）2233l1（b′3-l31y1-l32y2）2233）?）lyn=1lnn

（b′n

n-1- njy- j=1求出y以后，x就可以很容易根據(jù)如下方程求得：PAGEPAGE10 電力系統(tǒng)分析中的計算方法（2版）第第2章 線性方程組的解法 PAGE11〓1 u12 u13 … u1n〓〓〓〓〓〓0 1 u23 … u2n〓〓x2〓〓

〓y1〓〓y2〓〓n〓〓〓〓n〓〓

〓〓x3〓=〓y3〓〓? ? ? ? ?〓〓

〓?〓〓0 0 0 ? 1類似地，解向量x可以通過回代得到：xn=ynxn-1=yn-1-un-1，nxn

〓

n

〓

n〓xn-2=yn-2-un-2，nxn-un-2，n-1xn-1n?nx1=y1-j=2u1jxjLU分解的價值在于，一旦矩陣A被分解成上三角矩陣和下三角矩陣，那么尋找解向量x就是一件直截了當(dāng)?shù)氖铝?。注意上述求解過程并沒有顯式地對矩陣A求逆。存在多種方法對矩陣進行LU因子分解，每種方法具有各自的優(yōu)缺點。一種常用的LU分解方法被稱為Croutb算法[6]。定義矩陣Q為〓l11 u12 u13 … u1n〓〓l21 l22 u23 … u2n〓〓〓Q?L+U-I=〓l1 l2 l3 … u3n〓〓

（2.24）〓? ? ? ? ?〓〓ln1 ln2 ln3 … lnn〓Crout算法逐列逐行地計算Q的元素，首先計算Q的列元素，然后計算Q的行元素，如圖2.1所示。Q中的每個元素qij只依賴于A的元素aij和前面已經(jīng)計算得到的Q的元素。對矩陣A進行LU分解的Crout算法對矩陣Q進行初始化，即清零，并令j1；根據(jù)下式計算Q的第j列（即矩陣L的第j列）：

圖2.1矩陣Q列和行元素的計算次序qkj

=akj

j-1- - i=1

對k=j，…，n （2.25）如果jn，則停止；假定qjj≠0，根據(jù)下式計算Q的第j行（即矩陣U的第j行）：jkqjjikq =1jkqjjik

-j-1q

） 對k=j+1，…，n （2.26）jki=ji令jj1jki=ji一旦完成LU分解，通過“前代”就能得到啞向量y：q- yk=q- kk

（bk

k-1kjj=1

） 對k=1，…，n （2.27）n類似地，通過“回代”就能得到解向量x：nxk=yk-=k+1qkjxj 對k=n，n-1，…1 （.28）衡量LU分解過程計算量的一種指標(biāo)是，獲得解向量所需要的乘法和除法運算次數(shù)，因為兩者都是浮點數(shù)運算。計算Q的第j列（L的第j列）需要的乘法和除法次數(shù)為n nj=1k=j（j-1）類似地，計算Q的第j行（U的第j行）需要的乘法和除法次數(shù)為∑∑

n-1 njj=1k=j+1nj=1jnnj=1（n-j）n這樣，可以算出LU分解所需要的乘法和除法次數(shù)為31（n3-n）3而前代和回代所需要的乘法和除法次數(shù)為n2。因此，求解式（2.1）線性方程組所需要的總的乘法和除法次數(shù)為31（n3-n）+n2 （2.29）3將這個乘除法次數(shù)與采用Cramer法則所需要的2（n+1）！次乘除法相比，對于任何較大階數(shù)的矩陣，顯然采用LU分解和前代/回代方法的計算效率要高得多。例2.3 采用LU分解和前代/回代算法，求解例2.2中方程的解。解2.3 第1步是求矩陣A的LU因子：〓1 3 4 8〓A〓2 1 2 3A〓4 3 5 8〓〓〓〓9 2 7 4〓〓〓從j=1開始，式（2.25）表示Q的第1列與A的第1列完全相同。類似地，根據(jù)式（2.26），Q的第1行變?yōu)閝12=a12=3=3q11 1q13=a13=4=4q11 1q14=a14=8=8q11 1這樣，對于j=1，矩陣Q變?yōu)?/p>

〓1 3 4 8〓Q〓2 Q〓4 〓〓9 〓對于j=2，將分別計算矩陣Q主對角元下面的列元素和右面的行元素。對于Q的第2列有q22=q22-q21q12=1-（2）（3）=-5q32=a32-q31q12=3-（4）（3）=-9q42=a42-q41q12=2-（9）（3）=-25Q中的每個元素通過A中的對應(yīng)元素和Q中已經(jīng)求得的元素得到。注意乘積的第2個下標(biāo)總是相同的，而第1個下標(biāo)與正要計算的元素下標(biāo)相同。這個規(guī)則對于列計算和行計算都是成立的。Q的第2行計算如下：q23=q1（a23-q21q13）=1（2-（2）（4））=622 -5 5q24=（a24-q21q14）=1（3-（2）（8））=1322 -5 5對j=2處理完以后，矩陣Q變?yōu)椤? 3 4 8〓5〓2 -5 6 13〓5Q=〓 5 〓4 -9 〓〓9 -25 〓繼續(xù)計算j=3，Q的第3列計算如下：〓

6〓 1〓〓q33=a33-（q31q13+q32q23）=5-〓（4）（4）+（-9）5〓=-5〓〓q43a43（q41q13q42q23）7〓（9）（4）（25）〓1而Q的第3行變?yōu)?

〓 5〓q34=q33（a34-（q31q14+q32q24））（5）〓8〓（4）（8）（9）〓13〓〓〓3〓 〓 〓5〓〓〓得到〓1 3 4 8〓5〓2 -5 6 13〓5〓〓Q=〓 51 〓〓〓〓4 -9 -5 3〓〓9 -25 1 〓最后，計算j=4，最后一個對角元為q44=a44-（q41q14+q42q24+q43q34）4〓（9）（8）（25）〓13〓（3）（1）〓6這樣，

〓 〓5〓 〓1 3 4 8〓Q〓 6 13〓Q〓2 -〓

51 5〓5〓〓4 -9 - 3〓5〓〓 〓〓9 -25 1 -6〓〓1 0 0 0〓〓2 -5 0 0〓〓L=〓4 -9 -1 0〓〓9 -25 1 -〓9 -25 1 -6〓 〓〓1 3 4 8〓〓0 1 6 13〓U=〓 5 5〓〓0 0 1 3〓〓〓〓0 0 0 1〓〓〓檢查上述計算是否正確的一種方法是驗證LU是否等于A，對于本例，結(jié)果是成立的。1一旦求得了矩陣LU，下一步就是基于矩陣L和向量b通過前代來計算啞向量y，采用前代解方程Ly=b得到y(tǒng)為1y1=b1y1=

1=1y2=（b2LL1y1）=（1-（2）（1））=122 -5 5y3（b3（L31y1L32y2））（〓1〓（4）（1）（9）1〓〓6L33 〓 44y4=（b4-（L41y1y2+L43y3））44

5〓〓=〓1-〓（9）（1）+（-25）〓1〓+（1）（6）〓〓 ==〓 這樣

-6〓5〓〓1〓〓1〓y〓5y〓〓〓6〓〓〓

〓〓 322〓3〓2類似地，采用回代解方程Ux=y，可以得到解向量x為x4=y4=32 〓3〓 3〓〓x3y3U34x46（3）〓2〓=2〓〓x2=y2-（U24x4+U23x3）=1

〓〓〓〓〓6〓〓=-115x1=y1-（U14x4+U13x3+U12x2）

〓〓5〓〓

〓5〓〓2〓〓 21〓（8）〓3〓（4）〓3〓（3）〓11〓〓=-1〓最終得到解向量

〓2

〓2〓〓

〓1

2〓〓 2x=1〓11〓2〓 3〓〓 3〓此解與例2.2中采用Gauss消去法和回代方法所得到的解是一樣的。檢查所得到的解是否正確的一種快速方法是將解向量x代回到線性方程組Ax=b中，看是否成立。2.2.1 采用部分選主元的LU分解法〓10-10 1〓〓x〓10-10 1〓〓x1〓=〓〓2〓〓2 1〓〓x2〓

〓5〓 （2.30）通過觀察，很容易得出此線性方程組的解為x12x21而A的LU分解因子為

L=〓10-10 0 〓〓2 （1-2×1010）〓U=〓1 1010〓y得到：

〓0 1〓10y1=10y= 2（5-2y= 2（1-2×1010）x采用回代方法求解Uy=x得到x2

=y21x1=100-100x20這里，x2的解是正確的，而x1的解是完全不對的。為什么會發(fā)生這種情況呢？問題是式（230）1010對大多數(shù)計算機來說是太接近于零了。但是，如果將式（2.30）進行整理，變?yōu)?2 1 x 5〓2〓10-0 1〓〓x1〓=1 （〓2然后，再進行LU分解，則得到2 0L=〓 1 〓〓10-10

〓1-2×10-10〓〓〓 〓 〓〓〓1 1〓U=〓 2〓〓0 1〓啞向量y變?yōu)?/p>

y1=52〓 52

-10〓=〓1-2×10 〓1y2 〓 1 〓〓1-2×10-10〓通過回代，得到x為

〓x212x1≈52

〓2-1（1）=22這與通過觀察得到的解相同。因此，即使對角元不是精確等于零，將方程組的次序進行重新排列以使最大數(shù)值的元素位于對角元上，仍然是一種好的做法。這種做法被稱為選主元”，并導(dǎo)致了式（218）中的置換矩陣P。由于Crout算法從第1列和第1行開始逐列和逐行計算矩陣Q，因此只能實施部分選主元”，也就是只有Q（對應(yīng)地A）的行可以相互交換，而列必須保持不變。為了選擇最好的主元，需要在第j個對角元（在LU分解的第j步）下面的列元素中選擇絕對值最大的元素，然后將該元素所在行與第j行相互交換。上述選主元策略可以簡潔地表述為部分選主元策略（1）在LU分解的第j步，選擇第k行作為交換行，k的選擇滿足如下條件：（2）

|qjj|=max|qkj| 對k=j，…，n （2.32）將第j行與第k行進行交換，同時更新A、P和Q。。置換矩陣P是由0和1構(gòu)成的矩陣，它等于一系列初等置換矩陣Pj，k的乘積，這里Pj，k表示第j行與第k行相交換的初等置換矩陣。初等置換矩陣Pj，k如圖2.2所示，它由單位矩陣通過交換第j行和第k行得到。通過左乘合適的Pj，k就能實現(xiàn)選主元的目的。因為這僅僅是行的交換，未知向量中的元素次序沒有變化。例2.4采用部分選主元方法重做例2.3。解2.4 為方便起見將A重

圖2.2初等置換矩陣Pj，k寫如下：

〓1 3 4 8〓2 1 2 3A=〓 〓〓4 3 5 8〓9 2 7 4〓對于j=1，Q的第1列就是A的第1列，利用式（2.32）的選主元策略，第1列中q41的絕對值最大，因此將第4行和第1行交換。對應(yīng)的初等置換矩陣P1，4為〓0 0 0 1〓P1，4〓0 1 0 0P1，4〓0 0 1 0〓1 0 0 0〓對應(yīng)的矩陣A變?yōu)?/p>

〓9 2 7 4〓A〓2 1 2 3A〓4 3 5 8〓1 3 4 8〓而j=1時的矩陣Q為

〓 2 79〓 9 99Q〓2〓4〓1

4〓9〓〓〓〓當(dāng)j=2時，計算Q第2列得到的結(jié)果是〓9 2 7

4〓〓 9 9 9〓〓2 5 〓〓=Q 〓 9 〓〓=4〓〓 19 〓4〓〓 9 〓〓1 25 〓〓 9 〓對第j列主對角元下面的元素進行最大值搜索，第j列（這里為第2列）的第4行元素絕對值最大，因此第2行與第4行交換，產(chǎn)生的初等置換矩陣P2，4為〓1 0 0 0〓P2，4〓0 0 0 1P2，4〓0 0 1 0〓〓0 1 0 0〓類似地，更新后的A為

〓9 2 7 4〓〓1 3 4 8〓〓4 3 5 8〓〓 〓而產(chǎn)生的矩陣Q為

〓2 1 2 3〓〓9 2 7 4〓〓 9 9 9〓〓1 25 29 68〓〓 9 25 25〓〓〓Q=〓4 19 〓〓〓〓 9 〓〓〓〓2 5 〓〓〓〓 9 〓當(dāng)j=3時，計算Q第3列得到的結(jié)果是〓9 2 7

4〓〓 9 9 9〓〓1 25 29 68〓〓〓〓 9 25 25〓〓〓Q=〓4 19

14 〓〓 9 25 〓〓〓〓2 5 -1 〓〓〓〓 9 5 〓這種情況下，主對角元具有最大的絕對值，因此不用再選主元。繼續(xù)計算Q的第3行得到〓9 2 7 4〓〓 9 9 9〓〓1 25 29 68〓〓〓〓 9 25 25〓〓〓Q=〓4 19

-

-12〓〓 9 25 14〓〓〓〓2 5 -1 〓〓〓〓 9 5 〓最后，計算q44得到最終的矩陣Q為〓9 2 7 4〓〓 9 9 9〓〓1 25 29 68〓Q=〓

9

25〓4〓 19 14 12〓4〓〓〓 9 -25 -14〓〓〓〓2 5 -1 3〓〓 9 5 7〓置換矩陣P等于上述2個初等置換矩陣的乘積：P=P2，4P1，4I〓0 0 0 1〓=〓1 0 0 0〓=〓0 0 1 0〓〓〓〓0 1 0 0〓〓〓上述結(jié)果可以通過檢查PA=LU而得到驗證。而前代和回代計算將對修正過的向量b′=Pb進行。2.2.2 采用完全選主元的LU分解法采用完全選主元的另一種LU分解方法是Gauss方法。這種方法將產(chǎn)生2個置換矩陣：一個用于行交換，如在部分選主元中那樣；另一個用于列交換。采用這種方法，LU因子滿足P1AP2=LU （2.33）因此，為了對線性方程組Ax=b進行求解，需要采用略微不同的方法。與部分選主元情況相同，置換矩陣P1左乘線性方程組得P1Ax=P1b=b′ （2.34）現(xiàn)在，定義一個新的向量z為

x=P2z （2.35）然后，將式（235）代入到式（234）中得P1AP2z=P1b=b′LUz=b′

（2.36）其中，式（2.36）可以通過前代和回代求解出z。一旦求出z，解向量x就可以根據(jù)式（235）求出。在全主元方法中，LU分解的每一步，都會對行和列進行交換以使絕對值最大的元素交換到對角元位置。主元從余下的矩陣元素中搜索，包括對角元下面的元素和對角元右面的元素。全主元策略（1）在LU分解的第j步，按如下條件選擇主元：（2）

|qjj|=max|qkl| 對k=j，…，n和l=j，…，n （2.37）交換對應(yīng)的行并相應(yīng)地更新A、P和Q。對A進行LU分解的Gauss算法將矩陣Q初始化為零矩陣，令j1；設(shè)令矩陣Q的第j列（矩陣L的第j列），使之等于殘余矩陣A（j）的第j列，其中A（1）=A，即kjqkj=a（j）kj

對k=j，…，n （2.38）如果jn，則停止。假定qjj≠0，設(shè)令Q的第j行（U的第j行）為a（j）qjk=

jkqjj

對k=j+1，…，n （2.39）A（j）A（j1）即ikikijjka（j+1）=a（j）-qqikikijjk

對i=j+1，…，n和k=j+1，…，n （2.40）設(shè)令jj1，轉(zhuǎn)到步驟（2）。Gauss矩陣LU分解算法的乘除法次數(shù)與Crout的LU分解算法相同。Crout算法對矩陣A的每個元素只計算1次，而Gauss算法的每一步都要對矩陣A的元素進行更新。Crout算法相比于Gauss算法的一個優(yōu)點是矩陣A的每個元素只使用一次。由于每個qjk是ajk的函數(shù)，而ajk以后再也不會被使用，因此元素qjk蓋掉元素ajk。這樣，就不需要在內(nèi)存中存儲2nn階矩陣（AQ），只要存儲一個矩陣就夠了。CroutGaussLU分解算法中的兩種其他算法還包括Doolittle算法和雙因子算法[20，26，49等。這些算法中的大多數(shù)具有相同的乘除法次數(shù)當(dāng)在傳統(tǒng)的串行計算機上實現(xiàn)時其性能只有微小的差別但是當(dāng)考慮開放內(nèi)存存儲量和并行計算時這些算法會有巨大差別因此應(yīng)當(dāng)明智地選擇LU分解的算法，使其適合所應(yīng)用的場合和所采用的計算機結(jié)構(gòu)。23 條件數(shù)與誤差傳播GaussLU分解法被認為是直接法因為它們通過有限步就能夠計xA1b，并不需要迭代過程。在具有無限精度的計算機上，直接x?但是由于計算機是有限精度的因而所獲得的解也只有有限精度。矩陣的條件數(shù)”是一種用來衡量解的精確程度的有用指標(biāo)。A的條件數(shù)一般定義為λmaxλminκ（A）= （2.λmaxλmin式中，λmax和λmin表示矩陣ATA的最大和最小特征值。不管矩陣A本身的特征值是實數(shù)還是復(fù)數(shù)，ATA的特征值是實數(shù)并且非負。矩陣的條件數(shù)是該矩陣特征向量線性獨立性的一種度量。奇異矩陣至少有一個零特征值，并且包含至少一個退化行（即該行可以通過其他行的線性組合來表示）。單位矩陣的所有特征值為1，其特征向量是最線性獨立的，它的條件數(shù)為1。如果一個矩陣的條件數(shù)大大超過1，那么就稱該矩陣是“病態(tài)的”。條件數(shù)越大，求解過程相對于A的元素的微小變化就越敏感，解向量包含數(shù)值誤差的可能性就越大。由于在求解過程中引入的數(shù)值誤差，計算得到的式（2.1）的解x將不同于其精確解x?，設(shè)兩者存在的誤差為Δx。其他誤差，如近似誤差、測量誤差和舍入誤差等，也會引入到矩陣A和向量b中。Gauss消去法產(chǎn)生的解大致具有下式描述的十進制正確位數(shù)：tlog10βlog10κ（A） （2.42）式中，t是尾數(shù)的位長（對于典型的32位二進制字，t24），β是基（對二進制運算，β2），而κ是矩陣A的條件數(shù)。式（242）的一種解釋是，在Gauss消去過程（因此也是LU分解過程）中，方程的解將會損失大約log10κ的精確位數(shù)。基于矩陣元素的已知精度、條件數(shù)和機器精度，可以對數(shù)值解x的精度做出預(yù)測[35。24 松弛法松弛法本質(zhì)上是迭代法，它產(chǎn)生一系列向量，理想條件下會收斂于解x1b松弛法可以以多種方式應(yīng)用于求解式（21）。在所有情況下，采用松弛法的一個主要優(yōu)勢來自于不需要對一個大型線性方程組進行直接求解，并且可以有效地利用該方程組的一些潛在特性（這些特性目前是相對不變的）此外隨著并行處理技術(shù)的出現(xiàn)，松弛法比直接法更容易實現(xiàn)并行計算。最常用的兩種松弛法是Jacobi法和Gauss-Seidel法[56]。這些松弛法可以應(yīng)用于求解線性方程組：Ax=b （2.43）松弛法的一種通用做法是定義一個分裂矩陣”M，使得式（2.43）可以重寫成如下的等價形式：Mx=（M-A）x+b （2.44）這種分裂方法可以導(dǎo)出如下的迭代過程：Mxk1（MA）xkbk1，…，∞ （2.45）式中，k是迭代指標(biāo)。對于一個初始的猜想值x0，此迭代過程產(chǎn)生一系列向量x1x2通過選擇不同的矩陣M，可以導(dǎo)出不同的迭代方法。松弛法的目標(biāo)是選擇一個分裂矩陣M使得上述向量序列容易計算，并且能夠快速收斂到解。設(shè)A被分裂成L+D+U，其中L是嚴格下三角矩陣，D是對角矩陣，而U是嚴格上三角矩陣。注意，這些矩陣是不同于LU分解時所得到的L和U的。那么采用Jacobi松弛法求解向量x的迭代式為xk+1=-D-1（（L+U）xk-b） （2.46）n或等價地寫成標(biāo)量形式：nijaiixk+1ijaii

=-

〓aij〓xk+bi

1≤i≤n，k≥0 （2.47）j≠i〓aii在Jacobi松弛法中，向量xk+1所有元素的更新只用向量xj≠i〓aiiGauss-Seidel松弛法是類似的：xk+1=-（L+D）-1（Uxk-b） （2.48）或用標(biāo)量形式：ijj1〓aiiji1〓aiixk+1=-i-1ijj1〓aiiji1〓aii

〓aij〓xk

1≤i≤n，k≥0 （2.49）ijaiiGauss-Seidel法的優(yōu)勢是每個新的更新xk+1僅依賴于前面迭代已計算得到的值ijaii12i-xk+1、xk+1、…、xk+1。12i-值所占的位置上，從而減小了存儲的需求。由于松弛法是迭代的，確定在什么條件下能保證收斂到如下的精確解是必須的：x?=A-1b （2.50）眾所周知，對于任意的初始猜想值x0，Jacobi松弛法收斂的充分必要條件是如下矩陣的所有特征值落在復(fù)平面的單位圓內(nèi)[56]。MJ?-D-1（L+U） （2.51）類似地，對于Gauss-Seidel松弛法，相應(yīng)的收斂條件是如下矩陣的所有特征值落在復(fù)平面的單位圓內(nèi)。MGS?-（L+D）-1U （2.52）實際上，這些條件是難以確認的。但存在幾個容易確認的更一般性的條件，可以保證收斂。特別地，如果A是嚴格對角占優(yōu)的，那么不管是Jacobi松弛法還是Gauss-Seidel松弛法，都能保證收斂到精確解。初始向量x0可以是任意的，但是，如果存在解的更好的猜想值，就應(yīng)該使用之，這樣能夠更快速地收斂到指定精度下的精確解。一般地對大多數(shù)問題，Gauss-Seidel法比Jacobi法收斂速度快。如果A是一個下三角矩陣，Gauss-Seidel法可以通過1次迭代收斂到精確解，而Jacobi法則需要n次迭代才能收斂到精確解。但Jacobi法也有其優(yōu)點，即每次迭代計算時，ijixk1是獨立于所有其他的xk1（ji）的；這樣，所有xk1的計算可以并行處理。iji因此，這種方法非常適合于并行處理[36]。不管是Jacobi法還是Gauss-Seidel法，都可以被一般化為分塊Jacobi法和分塊Gauss-Seidel法。即A被分裂成分塊矩陣L+D+U，其中，D是分塊對角矩陣，L和U分別是分塊下三角矩陣和分塊上三角矩陣。分塊情況與標(biāo)量情況相同，也存在收斂的充分必要條件，即MJ和MGS的特征值必須落在復(fù)平面的單位圓內(nèi)。例2.5 采用（1）Gauss-Seidel法和（2）Jacobi法求解如下方程?！?10 2 3 6

〓1〓〓 0 -9 1 4〓x=〓2

（2.53）〓 2 6 -12 2〓〓〓 3 1 0 -8〓〓

〓3〓〓〓〓4〓〓〓解2.5 式的-，當(dāng)初始向量x=[0000]T時得到如下的迭代序列：kx1x2x3x41-0.0000-0.0000-0.0000-0.00002-0.1000-0.2222-0.3778-0.56533-0.5969-0.5154-0.7014-0.78834-0.8865-0.6505-0.8544-0.91375-1.0347-0.7233-0.9364-0.97846-1.1126-0.7611-0.9791-1.01247-1.1534-0.7809-1.0014-1.03018-1.1747-0.7913-1.0131-1.03949-1.1859-0.7968-1.0193-1.044310-1.1917-0.7996-1.0225-1.046811-1.1948-0.8011-1.0241-1.048212-1.1964-0.8019-1.0250-1.048913-1.1972-0.8023-1.0225-1.049214-1.1976-0.8025-1.0257-1.049415-1.1979-0.8026-1.0259-1.0495161.19800.80271.02591.0496Gauss-Seidel法已收斂到解x=[-1.1980 -0.8027 -1.0259 -1.0496]T采用由式（2.47）給出的Jacobi法，當(dāng)初始向量x=[0000]T時得到如下的迭代序列：kx1x2x3x41-0.0000-0.0000-0.0000-0.00002-0.1000-0.2222-0.2500-0.50003-0.5194-0.4722-0.4611-0.56534-0.6719-0.5217-0.6669-0.75385-0.8573-0.6314-0.7500-0.81766-0.9418-0.6689-0.8448-0.90047-1.0275-0.7163-0.8915-0.93688-1.0728-0.7376-0.9355-0.97489-1.1131-0.7594-0.9601-0.994510-1.1366-0.7709-0.9810-1.012311-1.1559-0.7811-0.9936-1.022612-1.1679-0.7871-1.0037-1.031113-1.1772-0.7920-1.0100-1.036314-1.1832-0.7950-1.0149-1.040415-1.1877-0.7974-1.0181-1.043116-1.1908-0.7989-1.0205-1.045117-1.1930-0.8001-1.0221-1.046418-1.1945-0.8009-1.0233-1.047419-1.1956-0.8014-1.0241-1.048020-1.1963-0.8018-1.0247-1.048521-1.1969-0.8021-1.0250-1.048922-1.1972-0.8023-1.0253-1.049123-1.1975-0.8024-1.0255-1.049224-1.1977-0.8025-1.0257-1.0494251.19780.80261.02581.0494Jacobi法所收斂到的解與GaussSeidel法相同。迭代過程中的誤差如圖23ii所示，該圖采用半對數(shù)刻度，而誤差定義為（xkx?）（i1，…，4）的最大值。iiGauss-Seidel法和Jacobi法都呈現(xiàn)出線性收斂特性，但Gauss-Seidel法的斜率更陡，因而在相同的初始值下能更快地收斂到誤差限內(nèi)。23GaussSeidelJacobi法的收斂速度例2.6 采用Jacobi迭代法重做例2.2。解2.6 再次采用例2.5中的迭代過程，得到如下的Jacobi迭代序列：kx1x2x3x41000021.00001.00000.20000.25003-4.8000-2.1500-1.6000-2.8500436.650022.35009.890014.92505-225.0100-136.8550-66.4100-110.6950顯然，上述迭代序列是不收斂的。為了理解為什么它們是發(fā)散的，考察Ja-迭代矩陣：MJ=-D-1（L+U）〓 0.00 -3.00 -4.00 -8.00〓=〓-2.00 0.00 -2.00 -3.00〓=〓-2.25 -0.50 -1.75 0.00〓〓-2.25 -0.50 -1.75 0.00MJ的特征值為

〓 〓〓66212〓〓 4.3574〓〓 1.2072〓〓 1.0566〓都大于1并落在單位圓外。因此，不管取什么初始值，Jacobi法都不能收斂，因而不能用來求解例2.2的方程組。如果迭代矩陣MJ和MGS的最大特征值小于1但幾乎等于1，那么收斂過程將非常慢。這種情況下，希望引入一個加權(quán)因子ω來改善收斂的速度。根據(jù)xk+1=-（L+D）-1（Uxk-b） （2.54）得到xk+1=xk-D-1（Lxk+1+（D+U）xk-b） （2.55）帶有加權(quán)因子ω的一種新的迭代方法為xk1xkωD1（Lxk1（DU）xkb） （2.56）這種方法被稱為逐次超松弛（SOR）”法，其中松弛系數(shù)ω0。注意，如果松弛迭代收斂，它們會收斂到解xA1b。SOR法收斂的一個必要條件是0ω<2[27]。除了很少的簡單情況外，計算ω的最優(yōu)值是困難的。通常通過試錯的方法確定最優(yōu)值，但是分析表明，對于n30的線性方程組，最優(yōu)的SOR方法可以比Jacobi方法快40倍[27]。隨著n的增大，對收斂速度的改進將更加明顯。25 共軛梯度法另一種求解方程組Ax=b的常用迭代法是“共軛梯度”法。這種方法可以認為是對如下函數(shù)逐次沿射線方向進行最小化的方法。E（x）=‖Ax-b‖2 （2.57）這種方法的一個有吸引力的特點是，如果矩陣A是正定的，那么可以保證至多在n步內(nèi)收斂（忽略舍入誤差）。如果矩陣A非常大且是稀疏的，共軛梯度法比Gauss消去法要常用得多，這種情況下可以在不到n步內(nèi)得到解，特別是當(dāng)矩陣A的性態(tài)很好時。如果矩陣A是病態(tài)的，那么舍入誤差可能會導(dǎo)致n步迭代后仍然得不到足夠精確的解。在共軛梯度法中，每次迭代需要確定搜ρk和參數(shù)αk，以使函數(shù)f（xkαkρk沿著ρk方向最小化。一旦設(shè)定xk1等于xk-αkρk，新的搜索方向就確定了。隨著共軛梯度迭代的進展，每個誤差函數(shù)都與一個特定的射線即正交擴展相關(guān)。因此，共軛梯度法可簡化為產(chǎn)生正交向量并尋找合適系數(shù)來表示所期望解的過程。共軛梯度法可以用24來說明，設(shè)x表示精確解（但未知），xk是一個近似解，而Δxk=xk-x?，給

圖2.4共軛梯度法定任何搜索方向ρk，從此直線到x的最小距離可以通過構(gòu)建垂直于此直線而連接x?的Δxk1來實現(xiàn)。由于精確解是未知的，因而構(gòu)造的殘差須與ρk垂直。不管新的搜索方向怎樣選擇，此殘差的范數(shù)將不會增加。所有求解Ax=b的迭代法都定義一個如下的迭代過程：xk1xkαk1ρk1 （2.58）式中，xk1是需要更新的值，αk是步長，ρk定義一個在Rn空間的方向，算法沿著此方向更新估計值。設(shè)第k步的殘差向量為而誤差函數(shù)由下式給出：

rk=Axk-b （2.59）Ek（xk）=‖Axk-b‖2 （2.60）那么在k+1步使誤差函數(shù)最小化的系數(shù)為=ak+1=

‖ATrk‖2‖Aρk+1

（2.61）上式的幾何解釋是沿著ρk+1定義的射線方向使Ek+1最小化。更進一步，一種改進的算法是在由兩個向量張成的平面內(nèi)尋找Ek+1的最小值，即xk+1=xk+αk+1（ρk+1+βk+1σk+1） （2.62）式中，射線ρk+1和σk+1在Rn空間張成一個平面。當(dāng)所選擇的向量正交，即滿足Aρk1Aσk10 （2.63）時，為了使誤差函數(shù)Ek1最小化而選擇的方向向量和系數(shù)是最優(yōu)的。式（2.63）中的表示內(nèi)積。滿足式（263）正交條件的向量被稱為關(guān)于算子ATA相互共軛，其中AT是A的共軛轉(zhuǎn)置。選擇合適向量的一種方法是選擇σk1與ρk正交，這樣就不需要在每一步中確定兩個正交向量。雖然這簡化了迭代過程，但在產(chǎn)生向量ρ時存在隱性的遞歸依賴性。求解Ax=b的共軛梯度算法初始化：令k=0，和如果‖rk‖≥ε，則

r0=Ax0-b （2.64）ρ0=-ATr0 （2.65）αk+1

‖ATrk‖2‖Aρk‖2

（2.66）=xk+1=xk+αk+1ρk （2.67）== +rk+1=Axk+1-b = +Bk+1

‖ATrk‖2‖ATrk‖2

（2.69）ρk+1=-ATrk+1+Bk+1ρk （2.70）k=k+1 （2.71）對于任意的非奇異正定矩陣A，共軛梯度法至多在n步內(nèi)得到解（忽略舍入誤差）。這是由n個方向向量ρ0、ρ1、…必定會張成一個解空間而得出的直接結(jié)果。有限步終止是共軛梯度法相對于諸如松弛法等其他迭代法的一個巨大優(yōu)勢。例2.7用共軛梯度法重新計算例2.5。解2.7為方便起見將例2.5的問題重寫如下：〓-10 2 3 6

〓1〓〓 0 -9 1 4〓x=〓2

（2.72）〓 2 6 -12 2〓 3 1 0 -8

〓3〓〓4〓〓并且x0=[0000]T初始化：令k=1，并且

〓 〓〓〓〓〓3r0Ax0b〓〓3〓4

（2.73）〓ρ0ATr0〓

〓86

（2.74）〓12〓〓12〓 〓初始的誤差為第1次迭代

E0=‖r0‖2=（5.4772）2=30 （2.75）‖ 01α=‖ATr0‖2=0.‖ 01〓x1x0α1ρ0〓

0.038900292

（2.77）〓0〓00583〓 〓〓21328〓r= xr= x-〓10553〓〓〓〓33874〓〓〓

（2.78）‖T r01B=‖ATr‖T r01〓7〓〓86195ρ1ATr1B1ρ0〓〓86195〓35414

（2.80）〓 〓對應(yīng)的誤差為E1=‖r1‖2=（4.9134）2=24.1416 （2.81）類似地，第2～4次迭代如下：kαkxkrkBkρk‖rk‖20.0464〓03211〓03820〓05504〓-0.2225〓15391〓00030〓0225435649〓2.5562〓249948〓174017〓147079〓-28.77653.8895（續(xù)）kαkxkrkBkρk‖rk‖30.02370.0085〓〓09133〓07942〓08988〓-0.904311981〓08027〓10260〓-1.0496〓15779〓06320〓06146〓-0.2998〓00000〓〓0.0000〓0.0000〓〓0.0000〓0.79490.0000〓335195〓10021〓149648〓-17.1040〓00000〓〓0.0000〓〓0.0000000001.83220.00004正如該算法所保證的，迭代在第4步收斂。不幸的是，對于一般性的線性方程組，共軛梯度法所需要的乘除法次數(shù)大大超過LU分解法。共軛梯度法在如下兩種情況下更有競爭力：第1種情況是矩陣非常大并且很稀疏；第2種情況是矩陣具有特殊結(jié)構(gòu)，LU分解不容易處理。在有些情況下，共軛梯度法的收斂速度可以通過“預(yù)處理”而得到改進。正如在Gauss-Seidel和Jacobi迭代中所看到的，迭代算法的收斂速度與迭代矩陣的特征值頻譜緊密相關(guān)。因此，采用重新定標(biāo)或矩陣變換的方法將原始方程組變換成具有更好矩陣特征值頻譜的方程組，將能夠大大提高收斂的速度。這個過程被稱為預(yù)處理。已開發(fā)出了針對稀疏矩陣的多個系統(tǒng)性的預(yù)處理方法，其中最基本的做法是將方程組變換成等價的方程組

Ax=bM-1Ax=M-1b式中，M-1用來近似A-1。例如，此種方法可用于非完整的LU分解，即使用LU分解但所有的填入被忽略。如果矩陣A是對角占優(yōu)的，M-1的一個簡單近似是〓〓〓M1〓〓〓〓

1A（1，1）A 1 A（2，2）

〓〓〓〓〓〓? 1 A（n，n）〓

（2.82）這種預(yù)處理策略將線性方程組重新定標(biāo)，使得對角線上的元素全部相等。這種處理可以補償元素大小上的數(shù)量級的差別。注意，重新定標(biāo)對矩陣固有的病態(tài)特性并沒有任何作用。例2.8 采用預(yù)處理的方法重做例2.7。解2.8 設(shè)M-1采用式（2.82）的形式，因此有〓-1 〓1〓 10 1-〓 -M-1〓 9 M-1-〓〓〓 1 -〓〓〓〓〓 12 〓〓- - 8

（2.83）〓AM1A

10000 02000 03000 0600000000 10000 01111 04444

（2.84）〓-0.3750 -0.1250 0.0000 1.0000〓〓-0.3750 -0.1250 0.0000 1.0000〓〓01000〓b′= bb′= b〓02500〓〓〓〓05000〓〓〓

〓（2.85）求解A′x=b′，采用共軛梯度法得到如下的誤差序列：k Ek0 0.60981 0.55002 0.35593 0.1131 4 0.0000雖然也要用4步迭代才收斂，但注意在k<4時，其誤差與例2.7的情況相比要小很多。對于大型線性方程組，可以想象誤差將會足夠快速地下降，以至于在第n步迭代之前就可以終止。如果只需要x的一個近似解，這種方法也是很有用的。6 廣義最小殘差法如果矩陣A既不是對稱的，也不是正定的，那么就不能保證下式<Aρk+1，Aσk+1>為零，從而搜索向量不是相互正交的。為了構(gòu)成解空間的基，相互正交是必須的。因此，這個基必須采用顯式方法構(gòu)造。共軛梯度方法的擴展，被稱為廣義最小殘差法（GMRES），在如下向量張成的子空間中使殘差的范數(shù)最小化：r0，Ar0，A2r0，…，Ak-1r0式中，向量r0是初始殘差r0=‖b-Ax0‖，而解的第k次近似就是從這個空間中選擇的。這個子空間是一個“Krylov子空間”，通過著名的Gram-Schmidt方法將其正交化；當(dāng)應(yīng)用于Krylov子空間時，這個正交化的過程被稱為Arnoldi算法[37]。在第k步，GMRES方法將Arnoldi算法應(yīng)用于構(gòu)成第k個Krylov子空間的k個正交基向量上，以產(chǎn)生下一個基向量。Arnoldi算法將在7.3節(jié)更詳細地描述。每步計算時，用矩陣A乘前面已求出的Arnoldi向量vj，然后，將所得到的向量wj對所有前面已求出的向量vi正交化。列向量集Vv1v2，…，vk]構(gòu)成了Krylov子空間的正交基，而H是矩陣A在此子空間上的正交投影。像Arnoldi算法那樣的正交矩陣三角化需要確定一個n×n正交矩陣Q，使得0QT〓0

（2.86）式中，R是一個m×m上三角矩陣。這樣，求解過程就簡化為求解三角矩陣方程Rx=Py，這里的P由Q的前m行組成。為了將矩陣化為上三角矩陣，可以應(yīng)用Givens旋轉(zhuǎn)變換一次清除一個元素。Givens旋轉(zhuǎn)變換是基于如下矩陣的變換：〓1 … 0 … 0 … 0〓〓? ? ? ?〓〓〓〓0 … cs … sn … 0〓〓〓〓? ? ? ?〓〓〓0 … -sn … cs … 0〓〓

（2.87）〓? ? ? ?〓0 … 0 … 0 … 1〓式中，cs=cos（?），sn=sin（?），?為經(jīng)過適當(dāng)選擇的旋轉(zhuǎn)角。采用Givens旋轉(zhuǎn)變換可以將元素Aki變?yōu)榱恪MRES方法的困難之一是隨著k的增加，需要存儲的向量個數(shù)增加，而乘法次數(shù)等于k2n/2（對于n×n矩陣）。為了克服這個困難，此算法可以按照迭代的方式應(yīng)用，即可以每隔m步重啟動一次，其中m是某個固定的整數(shù)。這種方法通常被稱為GMRES（m）算法。用于求解Ax=b的GMRES（m）算法初始化：令k=0，并且r0=Ax0-be1=[1 0 0 … 0]T當(dāng)‖rk‖≥ε和k≤kmax時，設(shè)置

j=1當(dāng)j<m時Arnoldi過程

v1=r/‖r‖s=‖r‖e1cs=[000…0]Tsn=[000…0]T）構(gòu)成矩陣H，使得H（i，j）（Avj）Tvi，i1，…，ji=）令w=Avj-∑j H（i，j）i=（c）令H（j+1，j）=‖w‖（d）令vj+1=w/‖w‖Givens旋轉(zhuǎn)）計算〓H（

〓

cs（i） sn（i）〓

H（i，j）

i=1，…，j-1〓H（i1，j））令）

〓-sn（i） cs（i）〓〓H（i+1，j）〓cs（）=j H（j，j） cs（）=H（j+1，j）2+H（j，j）2sn（）=j H（j+1，j） sn（）=H（j+1，j）2+H（j，j）2求殘差范數(shù)的近似值

α=cs（j）s（j））令

s（j+1）=-sn（j）s（j）s（j）=α誤差=|s（j+1）|H（j，j）=cs（j）H（j，j）+sn（j）H（j+1，j）H（j+1，j）=0（a）（（c）

如果誤差≤ε，則求解方程Hy=s，求出y作近似值更新此過程已收斂，返回。

x=x-Vy令j=j1如果jm并且誤差ε（重新啟動GMRES）（a）令s=[100…0]）解Hy=sy）計算x=xVy）計算rAxb）令k=k1例2.9 采用GMRES方法重做例2.5。解2.9 為了方便起見將例2.5重寫如下：求解方程〓-10 2 3 6

〓1〓〓 0 -9 1 4〓 2 6 -12 2x〓 3 1 0 -8x

〓2〓3〓4

（2.88）〓并且x0=[0000]T，ε=10-3

〓 〓〓j=1 采用Arnoldi算法得到-4.〓H〓〓〓

62369 0 0 0〓0 0 0 0〓0 0 0 0〓〓-0.1826 -0.9084 0 0〓V〓-0.3651 -0.2654 0 0V〓-0.5477 0.0556 0 0〓〓-0.7303 0.3181 0 0〓〓 〓應(yīng)用Givens旋轉(zhuǎn)變換得到cs=[-0.5430 0 0 0]sn=[0.8397 0 0 0]〓7.4274 0 0 0〓〓H〓〓〓

0 0 0 0〓0 0 0 0〓0 0 0 0〓s=[-2.9743 -4.5993 0 0]T由于誤差（=|s（2）|=4.5993）大于ε，j=j+1再重復(fù)。j2 采用Arnoldi算法得到〓7.4274 2.6293 0 0〓〓H〓〓〓

0 125947 0 0〓0 19321 0 0〓0 0 0 0〓〓-0.1826 -0.9084 -0.1721 0〓V〓-0.3651 -0.2654 0.6905 0V〓-0.5477 0.0556 -0.6728 0〓〓-0.7303 0.3181 0.2024 0〓〓 〓應(yīng)用Givens旋轉(zhuǎn)變換得到cs=[-0.5430 0.9229 0 0]sn=[0.8397 0.3850 0 0]〓7.4274 -12.0037 0 0〓〓H〓〓〓

0 50183 0 0〓0 0 0 0〓0 0 0 0〓s=[-2.9743 -4.2447 1.7708 0]T由于誤差（=|s（3）|=1.7708）大于ε，j=j+1再重復(fù)。j3 采用Arnoldi算法得到〓7.4274 -12.0037 -3.8697 0〓〓H〓〓〓

0 50183 02507 0〓0 0 131444 0〓0 0 26872 0〓〓-0.1826 -0.9084 -0.1721 -0.3343〓V〓-0.3651