數(shù)值分析 第5章 解線性方程組的直接方法

1、數(shù)理學(xué)院數(shù)理學(xué)院SCHOOL OF MATHEMATICS AND PHYSICS5 .1 引言與預(yù)備知識引言與預(yù)備知識5 .2 高斯消去法高斯消去法5 .3 高斯主元素消去法高斯主元素消去法5 .4 矩陣三角分解法矩陣三角分解法5 .5 向量和矩陣的范數(shù)向量和矩陣的范數(shù)5 .6 誤差分析誤差分析Ch5 解線性方程組的直接方法解線性方程組的直接方法 在自然科學(xué)和工程技術(shù)中有很多問題的解決常常歸結(jié)為解在自然科學(xué)和工程技術(shù)中有很多問題的解決常常歸結(jié)為解線性代數(shù)方程組如三次樣條函數(shù)問題，用最小二乘法求實(shí)驗(yàn)線性代數(shù)方程組如三次樣條函數(shù)問題，用最小二乘法求實(shí)驗(yàn)數(shù)據(jù)的曲線擬合問題，解非線性方程組問題，用差

2、分法或者有數(shù)據(jù)的曲線擬合問題，解非線性方程組問題，用差分法或者有限元方法解常微分方程、偏微分方程的邊值問題等都導(dǎo)致求解限元方法解常微分方程、偏微分方程的邊值問題等都導(dǎo)致求解線性代數(shù)方程組，而這些方程組的系數(shù)矩陣大致分為兩種，一線性代數(shù)方程組，而這些方程組的系數(shù)矩陣大致分為兩種，一種是低階稠密矩陣，另一種是大型稀疏矩陣。種是低階稠密矩陣，另一種是大型稀疏矩陣。 關(guān)于線性方程組的數(shù)值解法一般有兩類：關(guān)于線性方程組的數(shù)值解法一般有兩類：1.直接法直接法2.迭代法迭代法5.1引言與預(yù)備知識引言與預(yù)備知識5.1.1 引言引言 本章討論本章討論n元線性方程組元線性方程組nnnnnnnnnnbxaxaxab

3、xaxaxabxaxaxa.22112222212111212111 （5.1） 的直接解法。方程組的直接解法。方程組(5.1)的矩陣形式為的矩陣形式為 A Ax=b=b其中其中 nn2n1nn22221n11211a.aa.a.aaa.aaAn21x.xx,xn21b.bb,b若矩陣若矩陣A A非奇異非奇異,即即det(A A)0,則方程組則方程組(2.1)有唯一解。有唯一解。 所謂直接解法是指，若不考慮計(jì)算過程中的舍入誤差，所謂直接解法是指，若不考慮計(jì)算過程中的舍入誤差，經(jīng)過有限次算術(shù)運(yùn)算就能求出線性方程組的精確解的方法。經(jīng)過有限次算術(shù)運(yùn)算就能求出線性方程組的精確解的方法。 但由于實(shí)際計(jì)算

4、中舍入誤差的存在，用直接解法一般但由于實(shí)際計(jì)算中舍入誤差的存在，用直接解法一般也只能求出方程組的近似解。也只能求出方程組的近似解。 Cramer法則是一種不實(shí)用的直接法，本章將介紹幾種法則是一種不實(shí)用的直接法，本章將介紹幾種實(shí)用的直接法。實(shí)用的直接法。5.1.2 預(yù)備知識預(yù)備知識mnmmnijnmaaaaaaaaaaARA2124222111211)(M行行n列矩陣列矩陣. nnxxxxRx21n維列向量維列向量.列。同理的第為其中iAa，aaaAin),(21行行的第的第為為iAbbbATiTmT,1 矩陣的基本運(yùn)算矩陣的基本運(yùn)算:(1)矩陣的加法矩陣的加法是非奇異矩陣是非奇異矩陣AAdRA

5、RcAccAcRAAAbnnnnnT 0)det()(.,),det()det()(),det()det()(7)矩陣的行列式矩陣的行列式 行列式性質(zhì)行列式性質(zhì): (a) det(AB)=det(A)det(B)(6)非奇異矩陣非奇異矩陣(5)單位矩陣單位矩陣(4)轉(zhuǎn)置矩陣轉(zhuǎn)置矩陣(3)矩陣與矩陣的乘法矩陣與矩陣的乘法(2)矩陣與標(biāo)量的乘法矩陣與標(biāo)量的乘法5.1.3矩陣特征值與譜半徑矩陣特征值與譜半徑定義定義1設(shè)設(shè),nnijRaA若存在一個數(shù)若存在一個數(shù)(實(shí)數(shù)或復(fù)數(shù)實(shí)數(shù)或復(fù)數(shù))和非零向量和非零向量,),(21nTnRxxxx使使, xAx(1.1)則稱則稱為為A的特征值的特征值,x為為A對應(yīng)對

6、應(yīng)的特征向量的特征向量,A的全體特征值稱為的全體特征值稱為A的譜，記作的譜，記作記即.,)(),(21nAA稱為稱為A的譜半徑的譜半徑.|,|max)(1iniA(1.2)由式由式(1.1)知知, 可使齊次方程組可使齊次方程組0)(xAI有非零解有非零解,故系數(shù)行列式故系數(shù)行列式,0)det( AI記記0)det()(111212222111211nnnnnnnnnncccaaaaaaaaaAIp)(p稱為矩陣的特征多項(xiàng)式，方程稱為矩陣的特征多項(xiàng)式，方程(1.3)稱為的特征方程稱為的特征方程(1.3)在復(fù)數(shù)域中有在復(fù)數(shù)域中有n個根個根,21n故故).()()(21np由行列式由行列式(1.3)

7、展開可知：展開可知：ActrAacnnnnniiindet) 1() 1(,211211nnijRaA的的n個特征值個特征值 故故n,21是它的特征是它的特征方程方程(1.3)的幾個根的幾個根,并有并有.,det1121niiniiinatrAA(1.4)(1.5)A的跡的跡.A的特征值的特征值和特征向量和特征向量x還有以下性質(zhì)還有以下性質(zhì): (1) AT與與A有相同的特征值有相同的特征值 及相同的特征向量及相同的特征向量x . (2) 若若A非奇異，則非奇異，則A-1的特征值為的特征值為-1,特征向量為特征向量為x. (3) 相似矩陣相似矩陣 B=S-1AS有相同的特征多項(xiàng)式有相同的特征多項(xiàng)

8、式.例例1 求求242422221A的特征值及譜半徑的特征值及譜半徑.解解: A的特征方程為的特征方程為, 0)7()2(28243242422221)det(223AI故故A的特征值為的特征值為7, 2321A的譜半徑為的譜半徑為. 7)(A5.1.4 特殊矩陣特殊矩陣如果正交矩陣）對任意非零向量，（）如果（對稱正定矩陣如果設(shè)埃爾米特矩陣如果對稱矩陣時，如果當(dāng)上海森伯格陣時，如果當(dāng)上三角矩陣時，如果當(dāng)三對角矩陣時，如果當(dāng)對角矩陣)8(. 0),( ,)7()(,)6(.)5(. 01)4(. 0)3(. 01|)2(. 0) 1 (.)(AxxxAxRxbAAaAAAAAAajiajiaji

9、ajiRaATnTTHHnnTijijijijnnij到的矩陣。由初等置換陣的乘積得置換陣，列），得到的矩陣記為列與第第行（或交換行與第由單位矩陣交換第初等置換陣。如果設(shè)酉矩陣)11()10(,)9(1BAIAAjijiAACAijijijHnn定理定理1.，則下述命題等價：設(shè)nnRA.)()5(.)4(. 0)det()3(. 00)2(.) 1 (1nArankAAAxAxbAxRbn的秩存在只有唯一解齊次方程組有唯一解，方程組對任何定理定理2.為對稱正定陣，則：設(shè)nnRA)., 2 , 1(0)det(,)4(), 2 , 1(0)3(), 2 , 1(,), 2 , 1(,)2(.)

10、1 (11111nkAAniAnkaaaaAnkAAAAAkikkkkkkk即的順序主子式都大于零的特征值其中正定陣也是對稱則的順序主子陣為記也是對稱正定陣為非奇異矩陣，且定理定理3.), 2 , 1(0), 2 , 1(0det為對稱正定矩陣則的特征值或）（為對稱矩陣，如果設(shè)AniAnkARAiknn定理定理4 (Jordan標(biāo)準(zhǔn)型標(biāo)準(zhǔn)型) 設(shè)設(shè)A為為n階矩陣階矩陣,則存在一個非奇異矩陣則存在一個非奇異矩陣P使得使得)()()(22111rrJJJApP其中其中:）塊。為若當(dāng)（且JordannnrinJriiinniiiiiii.), 2 , 1( 1111(1)當(dāng)當(dāng)A的若當(dāng)標(biāo)準(zhǔn)型中所有若當(dāng)

11、塊的若當(dāng)標(biāo)準(zhǔn)型中所有若當(dāng)塊Ji均為一階時均為一階時,此標(biāo)準(zhǔn)型變成此標(biāo)準(zhǔn)型變成對角矩陣對角矩陣.）。，（陣其若當(dāng)標(biāo)準(zhǔn)型必為對角的特征值各不相同，則如果ndiagA21)2(求解求解bxA 高斯消去法高斯消去法( (逐次消去法逐次消去法) )及消去法和矩陣三角分解之及消去法和矩陣三角分解之間的關(guān)系：間的關(guān)系：5.2 高斯消去法高斯消去法,22112222212111212111nnnnnnnnbxaxaxabxaxaxabxaxaxann,2121212222111211nnnnnnnnbbbxxxaaaaaaaaa或例例2 用消去法解方程組用消去法解方程組 12254632132321xxxxx

12、xxx解解 第第1步步.將方程將方程(2.2)乘上乘上-2加到方程加到方程(2.4)上上,消去未知數(shù)消去未知數(shù)x1,得到得到(2.2)(2.3)(2.4).11432 xx第第2步步.將方程將方程(2.3)加到方程加到方程(2.5)上去上去,消去方程消去方程(2.5)中的中的x2,得到與原方程組等價的三角形方程組得到與原方程組等價的三角形方程組62546332321xxxxxx解為解為:Tx)3 , 2 , 1(* 首先舉一個簡單的例子來說明消去法的基本思想首先舉一個簡單的例子來說明消去法的基本思想.上述過程相當(dāng)于上述過程相當(dāng)于 65620014011111561401401111561221

13、40111)|(bA332331*)2(rrrrrr 思思路路首先將首先將A化為上三角陣化為上三角陣 /* upper-triangular matrix */，再回代求解再回代求解 /* backward substitution */。=消元消元記記,)()1()1(nnijaAA )1()1(1)1(.nbbbbStep 1：設(shè)設(shè) ，計(jì)算因子，計(jì)算因子0)1(11 a)., 2(/)1(11)1(11niaamii將增廣矩陣將增廣矩陣/* augmented matrix */ 第第 i 行行 mi1 第第1 1行行，得到，得到)1(1)1(1)1(12)1(11.baaan)2(A)

14、2(b).,2,()1(11)1()2()1(11)1()2(njibmbbamaaiiijiijij 其中其中Step k：設(shè)設(shè) ，計(jì)算因子，計(jì)算因子且計(jì)算且計(jì)算0)( kkka)., 1(/)()(nkiaamkkkkikik )., 1,()()()1()()()1(nkjibmbbamaakkikkikikkjikkijkij回代回代)()(/nnnnnnabx ) 1., 1()(1)()( niaxabxiiinijjiijiii只要只要 A 非奇異，即非奇異，即 A 1 存在，則可通過逐次消元存在，則可通過逐次消元及行交換，將方程組化為三角形方程組，求出唯一解。及行交換，將方程組

15、化為三角形方程組，求出唯一解。共進(jìn)行共進(jìn)行 ？ 步步 )()2(2)1(121)()2(2)2(22)1(1)1(12)1(11.nnnnnnnnbbbxxxaaaaaan 1注意注意2:設(shè)設(shè)Ax=b,其中其中A為非奇異矩陣為非奇異矩陣,如果如果, 011 a由于由于A為非奇異矩陣為非奇異矩陣,所以所以A的第一列一定有元素不等于零的第一列一定有元素不等于零.例如例如計(jì)算。斯消去法照樣可以進(jìn)行矩陣。繼續(xù)這過程，高階非奇異右下角矩陣為時然后進(jìn)行消元計(jì)算，這）位置，調(diào)到（將于是交換兩行元素（）（111), 02111nAarraiii定理定理5 設(shè)設(shè)Ax=b,其中其中nnRA(1)如果如果), 2

16、, 1(0)(nkakkk則可通過高斯消去法將則可通過高斯消去法將Ax=b約化為等價的三角形線性方程組約化為等價的三角形線性方程組(2.10),且計(jì)算公式為且計(jì)算公式為:)., 1(/)()(nkiaamkkkkikik消元計(jì)算消元計(jì)算(k=1,2,n-1)nkibmbbnkjiamaakkikkikikkjikkijkij., 1,., 1,)()()1()()()1(回代計(jì)算回代計(jì)算)()(/nnnnnnabx ) 1., 1()(1)()( niaxabxiiinijjiijiii(2)如果如果A為非奇異矩陣為非奇異矩陣,則可通過高斯消去法則可通過高斯消去法(及交換兩行的初等及交換兩行的

17、初等變換變換)將方程組將方程組Ax=b約化為方程組約化為方程組(2.10).定理定理6 約化的主元素約化的主元素aii(i) 0(i=1,2,k)的充要條件是矩陣的充要條件是矩陣A的所有的所有順序主子式順序主子式 /* determinant of leading principal submatrices */ Di0(i=1,2,k) .即即., 2 , 1, 0, 01111111kiaaaaDaDiiiii(2.12)推論推論 如果如果A的順序主子式的順序主子式Dk0(k=1,2, ,n-1),則則., 3 , 2,/1)(1)1(11nkDDaDakkkkk5.2.2 三角分解法三角

18、分解法 /* Matrix Factorization */ 高斯消元法的矩陣形式高斯消元法的矩陣形式 /* Matrix Form of G.E. */：Step 1:)0(/111111 aaamii記記 L1 =1.11121nmm ，則，則 )1()1(1bAL)1(1)1(1)1(11.baan)2(A) 2 (bStep n 1: )()2(2) 1(1)()2(2)2(22) 1(1) 1(12) 1(11121.nnnnnnnnnbbbaaaaaabALLL其中其中 Lk =1.11, 1knkkmm 1kL1.11, 1knkkmm 111211.nLLL111jiM,記為記

19、為L單位下三角陣單位下三角陣/* unitary lower-triangular matrix */記記 U =)()2(2)2(22)1(1)1(12)1(11.nnnnnaaaaaaLUA 定理定理7 若若A的所有順序主子式的所有順序主子式 /* determinant of leading principal submatrices */ 均不為均不為0，則，則 A 的的 LU 分解唯一（其分解唯一（其中中 L 為為單位單位下三角陣）。下三角陣）。證明：證明：由由1中定理可知，中定理可知，LU 分解存在。下面證明唯一性。分解存在。下面證明唯一性。若不唯一，則可設(shè)若不唯一，則可設(shè) A =

20、 L1U1 = L2U2 ，推出，推出121111UULL211122211LLUULLUpper-triangularLower-triangularWith diagonal entries 1I L 為一般下三角陣而為一般下三角陣而 U 為為單位單位上三角陣的分解稱為上三角陣的分解稱為Crout 分解分解。 實(shí)際上只要考慮實(shí)際上只要考慮 A* 的的 LU 分解，即分解，即 ，則，則 即是即是 A 的的 Crout 分解。分解。ULA* *LUA 121UU例例3 對于例對于例2,系數(shù)矩陣系數(shù)矩陣 122140111A由高斯消去法由高斯消去法,LUAmmm 2001401111120100

21、01, 1, 2, 0323121故故5.3.1 列主元消去法列主元消去法:在順序消元過程中在順序消元過程中,只要只要即可進(jìn)行計(jì)算即可進(jìn)行計(jì)算,) 1, 2 , 1(0)(nkakkk但如果但如果|)(kkka很小很小,則將導(dǎo)致舍入誤差增長則將導(dǎo)致舍入誤差增長,使解的誤差很大使解的誤差很大.例例4 用用Gauss 消去法求解方程組消去法求解方程組232120001. 02121xxxx5 .3 高斯主元素消去法高斯主元素消去法解解:因因099997. 3, 00001. 021故方程有唯一解故方程有唯一解,且精確解為且精確解為Tx)49998750. 0 ,25001875. 0(*若用若用G

22、auss消去法取四位有效數(shù)字計(jì)算消去法取四位有效數(shù)字計(jì)算,可得解可得解Tx)5000. 0 ,0000. 0(*xx與比較比較,誤差很大誤差很大,若將兩個方程互換為若將兩個方程互換為120001. 02322121xxxx用用Gauss消去法取四位有效數(shù)字計(jì)算消去法取四位有效數(shù)字計(jì)算,可得解可得解Tx)5000. 0 ,2500. 0( 本例表明通過行交換可避免舍入誤差增長本例表明通過行交換可避免舍入誤差增長,這就是列主元消這就是列主元消去法的基本思想去法的基本思想. 其計(jì)算步驟如下其計(jì)算步驟如下:第第1步步，在，在 )()1()1(ijaAA中的第中的第1列選主元，即列選主元，即 1111|

23、,|max|1iaainii行為主元行為主元,若若, 11i將 |)1()1(bA的第的第i1行與第行與第1行互換，再按消元行互換，再按消元公式計(jì)算得到公式計(jì)算得到|)2()2(bA假定上述過程已進(jìn)行假定上述過程已進(jìn)行(k-1)步，得到步，得到 |)()(kkbA第第k步，在步，在 )(kA中的第中的第k列選主元，即列選主元，即 |,|max|)()(kiknikkkiaak若若, kik則在則在 |)()(kkbA中將中將ik行與第行與第k行互換，再按消元行互換，再按消元公式計(jì)算得到公式計(jì)算得到|)1()1(kkbA對對k=1,2, ,n-1,重復(fù)以上過程則得重復(fù)以上過程則得|)()(nnb

24、A如果某個如果某個k出現(xiàn)主元出現(xiàn)主元 ），或0(0|)(kkka方程沒有唯一解或嚴(yán)重病態(tài)，否則可由方程沒有唯一解或嚴(yán)重病態(tài)，否則可由(3.2.4)求得解求得解.則表明則表明detA=0，它也表明當(dāng)非奇異時，存在排列矩陣它也表明當(dāng)非奇異時，存在排列矩陣(若干初等排列矩陣若干初等排列矩陣的乘積的乘積),使使PA=LU,其中其中L為單位下三角矩陣為單位下三角矩陣,其元素其元素|lij|=1,U為為上三角矩陣上三角矩陣.ULILILIAUAAILILILnniiiniininnn112211)(111212111121121上述每步行交換后再消元相當(dāng)于上述每步行交換后再消元相當(dāng)于1, 2 , 1,|)

25、1()1()()(1nkbAbAILkkkkkikk其中其中1kL是指標(biāo)為是指標(biāo)為k的初等下三角矩陣的初等下三角矩陣,kikI為初等排列矩陣為初等排列矩陣kik時時,IIkik表示不換行表示不換行,經(jīng)過經(jīng)過(n-1)步換行與消元步換行與消元,A化為上三角矩陣化為上三角矩陣.UAn)(即即:解解:例例5用列主元消去法解用列主元消去法解x=b,其中其中4178. 745625. 5996. 33816. 1078125. 014 . 022002. 0|bA4 . 022002. 03816. 1078125. 014178. 745625. 5996. 3|bA消元消元47471. 00010.

26、 161077. 0040371. 00020. 20028. 204178. 745625. 5996. 3|)2()2(bA消元消元35159. 039047. 00040371. 00020. 20028. 204178. 745625. 5996. 3|)3()3(bA消元結(jié)束由回代公式求得解消元結(jié)束由回代公式求得解Tx)90043. 0 ,69850. 0,92730. 1 (此例的精確解為此例的精確解為Tx)900423. 0 ,698496. 0,9300. 1 (*可見結(jié)果精度較高若不選列主元可見結(jié)果精度較高若不選列主元auss消去法，求得解消去法，求得解Tx)88888. 0

27、 ,68695. 0,930. 1 (，誤差較大，誤差較大除列主元消去法外，還有一種消去法，是在的所有元素除列主元消去法外，還有一種消去法，是在的所有元素aij中中選主元，稱為全主元消去法因計(jì)算量較大且應(yīng)用列主元已能選主元，稱為全主元消去法因計(jì)算量較大且應(yīng)用列主元已能滿足實(shí)際要求，故不再討論目前很多數(shù)學(xué)軟件庫都有列主元滿足實(shí)際要求，故不再討論目前很多數(shù)學(xué)軟件庫都有列主元消去法，可直接調(diào)用消去法，可直接調(diào)用注注:為了減少計(jì)算的舍入誤差為了減少計(jì)算的舍入誤差,使用消去法通常都要選主元使用消去法通常都要選主元.目前最目前最常用的是列主元消去法常用的是列主元消去法,也就是每步消元之前選主元也就是每步消

28、元之前選主元,當(dāng)當(dāng)A=(aij)第一步選第一步選A中第中第1列的主元列的主元,即即max|ai1|=ai1.然后將然后將i1行與第行行與第行互換，再進(jìn)行消元，以后每步消元做法類似，先選主元，再消互換，再進(jìn)行消元，以后每步消元做法類似，先選主元，再消元元5.3.2 高斯若當(dāng)消去法高斯若當(dāng)消去法消去對角線上方和下方的元素消去對角線上方和下方的元素.假設(shè)已經(jīng)完成假設(shè)已經(jīng)完成k-1步步,得到與方程得到與方程Ax=b等價的方程組等價的方程組1111, 1, 11111)()(11),(knnnnkkkknkkkknkkkknkkkbaabaabaabaabA 高斯消去法有很多變形高斯消去法有很多變形,有

29、的是高斯消去法的改進(jìn)、改寫，有的是高斯消去法的改進(jìn)、改寫，有的是用于某一類特殊性質(zhì)矩陣的高斯消去法的簡化。有的是用于某一類特殊性質(zhì)矩陣的高斯消去法的簡化。5.3.1 直接三角分解法直接三角分解法. 將高斯消去法改寫為緊湊形式將高斯消去法改寫為緊湊形式,可以直接從可以直接從A的元素得到計(jì)算的元素得到計(jì)算L,U元素的遞推公式元素的遞推公式,而不需要任何中間步驟而不需要任何中間步驟,這就是所謂直接三角這就是所謂直接三角分解法分解法. 一旦實(shí)現(xiàn)一旦實(shí)現(xiàn)A的的LU分解分解,那么求解那么求解Ax=b的問題就等價于求解兩的問題就等價于求解兩個三角形方程組個三角形方程組 (1) Ly=b,求求y; (2) U

30、x=y,求求x. yUxbLybLUx.4 矩陣矩陣三角分解法三角分解法 /* Matrix Factorization */.不選主元的三角分解法不選主元的三角分解法 設(shè)設(shè)A為非奇異矩陣為非奇異矩陣,且有分解式且有分解式A=LU,其中其中L為單位下三為單位下三角角,U為上三角即為上三角即 nnnnnnnnuuullaaaa.1.11.1111211111 L,U元素可以由元素可以由n步直接計(jì)算定出步直接計(jì)算定出,其中第其中第r步定出步定出U的第的第r行行和和L的第的第r列元素列元素.由上式有由上式有:.1), 2(/1), 2 , 1(1111111111列元素列元素的第的第得得行元素；行元

31、素；的第的第得得LniualulaUniuaiiiiii 故故), 1,(11nrriulaukirkrkriri 同樣有同樣有:rrirrkkriknkkrikirululula 111 設(shè)已經(jīng)定出設(shè)已經(jīng)定出U的第的第1行到第行到第r-1行元素與行元素與L的第的第1列到第列到第r-1列列元素元素.利用矩陣乘法利用矩陣乘法(注意當(dāng)注意當(dāng)rk時時,lrk=0),有有rikirkrkkinkrkriuulula 111得得irl通過比較法直接導(dǎo)出通過比較法直接導(dǎo)出L 和和 U 的計(jì)算公式。的計(jì)算公式。思思路路 nnnnnnnnuuullaaaa.1.11.1111211111 ),min(1jik

32、jkkiul jia固定固定 r : :對對 i = r, r+1, , n 有有rikirkrkriuula11lii = 1), 1,(11nrriulaukirkrkriria將將 r ，i 對換，對對換，對 r = i, i+1, , n 有有iijikiikjkjiulula 11);, 1(/ )(11nrnriuulalrrikkrikirir且b結(jié)論結(jié)論:用直接三角分解法解用直接三角分解法解Ax=b(要求要求A的所有順序主子式都不等的所有順序主子式都不等于于0)的計(jì)算公式如下的計(jì)算公式如下.Step 1: u1i = a1i; li1 = ai1 / u11; ( i = 2,

33、 , n )計(jì)算計(jì)算U的第的第r行行,L的第的第r列元素列元素(r=2,3,，n）Step 2:), 1,(11nrriulaukirkrkriri);, 1(/ )(11nrnriuulalrrrkkrikirir且求解求解Ly=b, Ux=y 的計(jì)算公式：的計(jì)算公式：Step 3:Step 4:1111);, 3 , 2(ikkikiiniylbybyStep 5:nikiikikiinnnnnniuxuyxuyx1).1 , 2, 1(/ )(;/例例5 用直接三角分解法解用直接三角分解法解201814513252321321xxx解解 用分解公式計(jì)算得用分解公式計(jì)算得.24004103

34、21153012001LUA求解求解.)3 , 2 , 1 (,)72,10,14(,)72,10,14(,)20,18,14(TTTTxUxyLy得得2.選主元的三角分解法選主元的三角分解法 當(dāng)當(dāng)urr=0時計(jì)算中斷，或者當(dāng)時計(jì)算中斷，或者當(dāng)urr絕對值很小時，按分解公式計(jì)絕對值很小時，按分解公式計(jì)算可能引起舍入誤差的積累。算可能引起舍入誤差的積累。 但如果但如果A非奇異，可以通過交換非奇異，可以通過交換A的行實(shí)現(xiàn)矩陣的行實(shí)現(xiàn)矩陣A的的LU分解，分解，因此可采用與列主元消去法類似的方法，將直接三角分解法修改因此可采用與列主元消去法類似的方法，將直接三角分解法修改為（部分）選主元的三角分解法。

35、為（部分）選主元的三角分解法。設(shè)第設(shè)第r-1步分解已完成，這時有步分解已完成，這時有nnnrrnnnrnrrrrrrnrrrrrrrnrrnrraalllaallluuulluuuuluuuuuA1,2,11,21, 1, 11, 12, 11 , 1221, 22221111, 11211第第r步分解需用到（步分解需用到（3.2）及（）及（3.3）式，為避免用小的數(shù)）式，為避免用小的數(shù)urr做除數(shù)，引進(jìn)量做除數(shù)，引進(jìn)量11)., 1,(rkkrikirinrriulas于是有：于是有：)., 1(/,nrisslsuriirrrr取取|maxirinirss 交換交換A的的r行與行與ir行元

36、素，將行元素，將irs調(diào)到調(diào)到(r,r)位置位置(將（將（i，j）位置的新元素仍記為）位置的新元素仍記為lij 及及 aii),于是有于是有| lir |1時時,aij=0,且：且：. 0| )();1, 3 , 2(0,|,| )(; 0| )(11 nniiiiiabcnicacabbcba nnnnnnnfffxxxbacbacbacb212111122211Step 1: 對對 A 作作Crout 分解分解 111121nnnA 直接比較等式兩邊直接比較等式兩邊的元素，可得到計(jì)的元素，可得到計(jì)算公式。算公式。).1, 3 , 2(), 2(,111111 nicnirbracbiiii

37、iiiii 為下三角，為單位上三角為下三角，為單位上三角注意當(dāng)注意當(dāng)j=1時有時有對對j=2，3，n求得求得L的元素的元素，這就是，這就是A的的Cholesky分解，然后再解兩個三角方程組，得分解，然后再解兩個三角方程組，得這就是對稱正定方程組的平方根法這就是對稱正定方程組的平方根法 另外，由于另外，由于 故有故有這表明分解過程中矩陣這表明分解過程中矩陣L中元素中元素因此平方根法計(jì)算是數(shù)值穩(wěn)定的因此平方根法計(jì)算是數(shù)值穩(wěn)定的. 11111111,lalalii ijlnilylbylbyiiikkiki, 3 , 2,/ )(,/111111 1 , 1,/ )(,/1 nilxlyxlyxjj

38、knikkiiinnnn的數(shù)量級不增長，的數(shù)量級不增長，), 2 , 1(0), 2 , 1(12njlnjaljjjjjkjk 而而jkljjnjjknknjal 111max|maxStep 2: 追追即解即解fyL ,111 fy ),.,2()()(111nibyafyrfyiiiiiiiiiii Step 3: 趕趕即解即解 ：yxU )1,., 1(,1 nixyxyxiiiinn xyUxyfLyfLUx，求，求求求,求解求解xf等價于等價于Step :計(jì)算計(jì)算 i 的遞推公式的遞推公式)/(,/1111 iiiiiabcbc 將計(jì)算系數(shù)將計(jì)算系數(shù)nnyyy 21121及及 為追

39、的過程為追的過程將計(jì)算方程組的解將計(jì)算方程組的解11xxxnn 為趕的過程為趕的過程定理：定理：設(shè)有三對角線方程組設(shè)有三對角線方程組xf，其中滿足對角占優(yōu)的條件，其中滿足對角占優(yōu)的條件，則為非奇異矩陣且追趕法計(jì)算公式中的則為非奇異矩陣且追趕法計(jì)算公式中的 ii ,滿足：|);1, 2( |02);1, 2 , 1( 1|01nnnnniiiiiiiababniababcni 。定理定理 若若 A 為為對角占優(yōu)對角占優(yōu) /* diagonally dominant */ 的三對角的三對角陣，且滿足陣，且滿足 ，則追趕，則追趕法可解以法可解以 A 為系數(shù)矩陣的方程組。為系數(shù)矩陣的方程組。0,0,

40、0|, 0|11 iinncaabcb如果如果 A 是是嚴(yán)格嚴(yán)格對角占優(yōu)陣，則不要求三對角線上的對角占優(yōu)陣，則不要求三對角線上的所有元素非零。所有元素非零。根據(jù)不等式根據(jù)不等式 可知：分解過程中，矩陣元素不會過分增大，算法可知：分解過程中，矩陣元素不會過分增大，算法保證穩(wěn)定。保證穩(wěn)定。運(yùn)算量為運(yùn)算量為 O(6n)。|,1|1iiiiiiiiabbab 5.5.1 內(nèi)積與向量范數(shù)內(nèi)積與向量范數(shù)為了研究方程組為了研究方程組Ax=b解的誤差和迭代法收斂性，需對向量解的誤差和迭代法收斂性，需對向量 nRx 及矩陣及矩陣 nnRA 的的大小大小引進(jìn)一種度量，就要定義范數(shù)，引進(jìn)一種度量，就要定義范數(shù)， 它

41、是向量它是向量長度長度概念的直接推廣，通常用概念的直接推廣，通常用 nR表示表示n維實(shí)向量維實(shí)向量空間，空間， nC表示表示n維復(fù)向量空間維復(fù)向量空間.定義定義2 設(shè)設(shè) ).(),(,),(2121nnTnTnCRyyyyxxxx或或 將實(shí)數(shù)將實(shí)數(shù)）或或復(fù)復(fù)數(shù)數(shù)iniiTiniiHyxxyyxyxxyyx 11),(),(稱為向量稱為向量x,y的數(shù)量積的數(shù)量積.非負(fù)實(shí)數(shù)非負(fù)實(shí)數(shù) ，或或21122122112212)|(),(|()(),(| niiniixxxxxxxx稱為向量稱為向量x的歐氏范數(shù)或的歐氏范數(shù)或2-范數(shù)范數(shù). 5.5 向量和矩陣范數(shù)向量和矩陣范數(shù)定理定理12設(shè)設(shè) ),(,nnC

42、Ryx或或 則內(nèi)積有以下性質(zhì)：則內(nèi)積有以下性質(zhì)： 時成立；時成立；當(dāng)且僅當(dāng)當(dāng)且僅當(dāng)0, 0),( xxx(1)(2)為為復(fù)復(fù)數(shù)數(shù)，為為實(shí)實(shí)數(shù)數(shù)（或或 ),(),().,(),(yxyxyxyx (3),(),()(,(),( xyyxxyyx或或（(4);,(),(),(2121yxyxyxx (5) (柯西柯西-施瓦茨不等式施瓦茨不等式),| ),( |22yxyx 等式當(dāng)且僅當(dāng)?shù)仁疆?dāng)且僅當(dāng)x與與y線性相關(guān)時成立線性相關(guān)時成立;(6) 三角不等式三角不等式 ,|222yxyx 定義定義3(向量范數(shù)向量范數(shù))如果向量如果向量 )C(nnRx或或 的某個實(shí)值函數(shù)的某個實(shí)值函數(shù)|,|)(xxN 滿

43、足條件滿足條件:三角不等式），三角不等式），）（），），或或）（），），時等號成立（正定條件時等號成立（正定條件當(dāng)且僅當(dāng)當(dāng)且僅當(dāng)(|3(, |20, 0|)1(yxyxCRxxxx 則稱則稱 .|)(上上的的一一個個向向量量范范數(shù)數(shù)是是nRxxN 對于對于 ,)(|2/1122 niixx由內(nèi)積性質(zhì)可知它滿足定義由內(nèi)積性質(zhì)可知它滿足定義2的三個條件，的三個條件， 故它是一種向量范數(shù)故它是一種向量范數(shù).此外還有以下幾種常用的向量范數(shù)此外還有以下幾種常用的向量范數(shù). 范數(shù)）范數(shù)）稱為稱為 ( |max|1inixx范數(shù)）范數(shù)）稱為稱為 1( |11niixx容易驗(yàn)證容易驗(yàn)證 1|xx及及 均滿足定

44、義均滿足定義2的三個條件的三個條件.更一般的還可定義更一般的還可定義pnipipxx/11)(| 但只有但只有p=1,2,時的三種范數(shù)是常用的向量范數(shù)時的三種范數(shù)是常用的向量范數(shù).例如給定例如給定 ,)3, 2 , 1(Tx 則可求出則可求出 3| ,14| , 6|21 xxx定理定理14設(shè)設(shè) |)(xxN 是是nR則則N(x)是向量是向量x的分量的分量 上任一種向量范數(shù)，上任一種向量范數(shù)，nxxx,21的連續(xù)函數(shù)的連續(xù)函數(shù). 定理定理15(向量范數(shù)的等價性向量范數(shù)的等價性)設(shè)設(shè) tS| 與與是是nR上任意兩種向量范數(shù)，則存在常數(shù)上任意兩種向量范數(shù)，則存在常數(shù) , 021 cc使使.|21S

45、tSxcxxc 5.5.2 矩陣范數(shù)矩陣范數(shù) 矩陣矩陣nnijRaA )(可看成可看成nn維向量，如果直接將向量維向量，如果直接將向量的的2-范數(shù)用于矩陣范數(shù)用于矩陣A，則可定義，則可定義2/11,2)(|)( njiijFaAAF稱為矩陣稱為矩陣A的的Frobenius范數(shù)，簡稱范數(shù)，簡稱F-范數(shù)范數(shù).它顯然滿足向它顯然滿足向量范數(shù)的三條性質(zhì)，但由于矩陣還有乘法運(yùn)算，因此矩陣量范數(shù)的三條性質(zhì)，但由于矩陣還有乘法運(yùn)算，因此矩陣范數(shù)的定義中應(yīng)增加新條件范數(shù)的定義中應(yīng)增加新條件.定義定義4 如果如果 nnijRaA )(的某個非負(fù)實(shí)函數(shù)的某個非負(fù)實(shí)函數(shù)N(A)，記作，記作A，滿足條件：滿足條件：;

46、,|,|)4(;,(|3;, |20, 0|)1(nnnnRBABAABRBABABARAAAA 三角不等式），三角不等式），）（）（），），時等號成立（正定條件時等號成立（正定條件當(dāng)且僅當(dāng)當(dāng)且僅當(dāng) 則稱則稱 .|)(上上的的一一個個向向量量范范數(shù)數(shù)是是nnRAAN 顯然顯然 FA|滿足定義中的四個條件滿足定義中的四個條件, (3)，(4)兩條均可由兩條均可由Cauchy-Schwarz不等式證明，故不等式證明，故 FA|是一種矩陣范數(shù)是一種矩陣范數(shù). 除矩陣自身的運(yùn)算外，在解方程中矩陣乘向量的運(yùn)算即除矩陣自身的運(yùn)算外，在解方程中矩陣乘向量的運(yùn)算即Ax，也是必不可少的也是必不可少的.因此要求所

47、引進(jìn)的范數(shù)應(yīng)滿足條件：因此要求所引進(jìn)的范數(shù)應(yīng)滿足條件：|xAAx 上式稱為相容性條件上式稱為相容性條件.為使引進(jìn)的矩陣范數(shù)滿足條件為使引進(jìn)的矩陣范數(shù)滿足條件(4.5)，我們給出以下定義我們給出以下定義.(4.5)定義定義6(矩陣的算子范數(shù)矩陣的算子范數(shù)) 設(shè)設(shè) ,nnnRARxnvnnnRxRARx是| ,當(dāng)給定向量范數(shù)當(dāng)給定向量范數(shù) v|時可定義時可定義 vxvvxvAxxAxAv|max|max|1|0稱為矩陣的算子范數(shù)或從屬范數(shù)稱為矩陣的算子范數(shù)或從屬范數(shù). (4.6)定理定理17 設(shè)設(shè) 上的一種向量范數(shù)，上的一種向量范數(shù)，則由則由(4.6)定義的定義的vA|是一種矩陣范數(shù)，且滿足相容性

48、條件是一種矩陣范數(shù)，且滿足相容性條件|xAAx 證明證明 因因 nvRA是|中有界閉集中有界閉集1| ,),(|1vTnxxxxxD上的連續(xù)函數(shù)，故上的連續(xù)函數(shù)，故 vA|在在D上有最大值，即上有最大值，即 Dx 0使使,|max|1|0vvxvAAxAxv而對而對 vnxxxxRx|, 0,令故故DxxAxxAxAxvvvvv,|所以所以 vvvvxvxAxxAxA|max|0從而當(dāng)從而當(dāng) |0 xAAxx 時成立，而成立，而x=0時時顯然也成立顯然也成立.|xAAx 定理定理17 設(shè)設(shè) ,nnnRARx則則)( |max|11的行范數(shù)AaAnjijni)( |max|111的列范數(shù)AaAn

49、iijnj)2( )(|2范數(shù)的AAAAT這里這里)(為矩陣的譜半徑為矩陣的譜半徑.例例7 已知已知.| ,| ,| ,|,432112AAAAAF求解解. 6| ,477. 530| , 7|1AAAF043020141410)det(,201414102AAIAATT465. 5866.29| ,22115212A從定理可以看出，計(jì)算從定理可以看出，計(jì)算 1|AA及較容易，而計(jì)算 2|A時因?yàn)橐髸r因?yàn)橐?AAT的特征值，所以較為困難的特征值，所以較為困難.但當(dāng)?shù)?dāng)A對稱時，有對稱時，有 ).()()()()(|max2maxmax22AAAAAAAATT定理定理19.|2）（為對稱矩陣

50、，則如果AARAnn定理定理18 對任何對任何| , nnRA為任一種從屬范數(shù)則為任一種從屬范數(shù)則|,|)(AA 反之，對任意反之，對任意0，至少存在一種從屬范數(shù)，至少存在一種從屬范數(shù) ,|使使.)(|AA證明證明:設(shè)設(shè) 為為A的特征值，則的特征值，則 ,Ax, 0 xx使由由|xAxxAAx及得得. |max)(|AAA即BIBRBnn則, 1| ,非奇異，且非奇異，且|11|)(|1BBI證明證明用反證法用反證法.假定假定(I+B)奇異，則齊次方程奇異，則齊次方程 0)(xBI有非零解有非零解 使即0,0 xRxn1|0000 xBxBxx故而而 1|max|000 xBxxBxBx與B1

51、的假設(shè)矛盾，故(I+B)非奇異. 又又 ，即IBIBBIIBIBI111)()(,)(得得，11)()(BIBIBI取范數(shù)取范數(shù) ，|)(|1|)(|11BIBBI得得|11|)(|1BBI定理定理20設(shè)設(shè) 5.6.1矩陣條件數(shù)與擾動方程組誤差界矩陣條件數(shù)與擾動方程組誤差界在解方程組在解方程組Ax=b時，由于各種原因，時，由于各種原因，A或或b往往有往往有誤差，從而使得解也產(chǎn)生誤差誤差，從而使得解也產(chǎn)生誤差.例例8 方程組方程組 0001. 880001. 626221xx的準(zhǔn)確解為的準(zhǔn)確解為 Tx)1 , 1(* ,當(dāng)當(dāng)A與與b有微小變化時，如變?yōu)榉匠逃形⑿∽兓瘯r，如變?yōu)榉匠?0002. 8

52、89999. 526221xx則準(zhǔn)確解為則準(zhǔn)確解為 ,)2,10(Tx 它表明它表明A,b的微小擾動引起方程解的微小擾動引起方程解x的很大的很大變化，這就是病態(tài)方程變化，這就是病態(tài)方程.5.6 誤差分析與病態(tài)方程組誤差分析與病態(tài)方程組定義定義7 求解線性方程組求解線性方程組Ax=b時，若時，若A或或b有微小擾動有微小擾動 時，時，或或|b| A解解x的誤差的誤差 | x 很大很大, ，則稱此方程組為，則稱此方程組為病態(tài)方程組，相應(yīng)的系數(shù)矩陣病態(tài)方程組，相應(yīng)的系數(shù)矩陣A稱為病態(tài)矩陣稱為病態(tài)矩陣. 反之，若此時反之，若此時 | x 很小很小, ，則稱此方程組為，則稱此方程組為良態(tài)方程組，相應(yīng)的系數(shù)

53、矩陣良態(tài)方程組，相應(yīng)的系數(shù)矩陣A稱為良態(tài)矩陣稱為良態(tài)矩陣. 注意方程組是否病態(tài)與用什么數(shù)值方法無關(guān)，它是注意方程組是否病態(tài)與用什么數(shù)值方法無關(guān)，它是由方程自身性質(zhì)決定的由方程自身性質(zhì)決定的. 在例在例8中因?yàn)樾辛惺街幸驗(yàn)樾辛惺?, 0102det4 A因此出現(xiàn)病態(tài)因此出現(xiàn)病態(tài).但有時但有時A從表面上看性質(zhì)很好，也可能是病態(tài)的從表面上看性質(zhì)很好，也可能是病態(tài)的.那么如何判斷那么如何判斷A是否病態(tài)是否病態(tài)?先給出如下定義先給出如下定義.例例9 方程組方程組Ax=b表示為表示為 3133233210957910685657787104321xxxx它的準(zhǔn)確解它的準(zhǔn)確解 ,)1 , 1 , 1 , 1

54、(*Tx A對稱正定且對稱正定且 . 1det A表面看性質(zhì)表面看性質(zhì)較好較好，但若對右端，但若對右端b作微小變化，如方程改為作微小變化，如方程改為 ,)9 .30, 1 .33, 9 .22, 1 .32()(TxxA 則解變?yōu)閯t解變?yōu)?.)1 . 1, 5 . 4 , 6 .12, 2 . 9(Txx 這里這里b的相對誤差大約只有的相對誤差大約只有 ,2001但解的相對誤差卻很大，但解的相對誤差卻很大，故故A也是病態(tài)矩陣也是病態(tài)矩陣.定義定義8 設(shè)設(shè) nnRA 非奇異，非奇異，v為矩陣的任一種從屬范數(shù)，則為矩陣的任一種從屬范數(shù)，則 ）或或 2 , 1(|)(1vAAACondvvv稱為矩陣

55、稱為矩陣A的條件數(shù)的條件數(shù). 從定義看到矩陣條件數(shù)依賴于范數(shù)的選取，如范數(shù)為從定義看到矩陣條件數(shù)依賴于范數(shù)的選取，如范數(shù)為2-范數(shù)，范數(shù)， 則記為則記為 .)()(|)(minmax2122AAAAAAACondTT 同理有同理有 1)(,)(ACondACond 等等等等.條件數(shù)有以下性質(zhì)：條件數(shù)有以下性質(zhì)： );()(, 1)(1 ACondACondACond(1)(2); 0,)()(1 RACondACond(3) U為正交矩陣，則為正交矩陣，則 ;)()()(, 1)(2222AUCondUACondACondUCond (4) 若若 1 與與n 為為A的按模最大與最小特征值，則的

56、按模最大與最小特征值，則.|)(1nACond 若若A對稱，則對稱，則 .|)(12nACond 下面給出擾動方程組解的誤差分析下面給出擾動方程組解的誤差分析.先考察先考察b有擾動有擾動 , b 則擾動方程為則擾動方程為0,)( bbbxxA 由于由于Ax=b,故得故得 ,1bAxbxA 于是于是 |,|1bAx 再由再由Ax=b,有有 |,|bxA 即即,|1bAx 故得故得.|)(|1bbACondbbAAxx 下面再研究方程下面再研究方程Ax=b，當(dāng)，當(dāng)A有擾動有擾動 A 時，其解時，其解 xx 的誤差分析的誤差分析. 此時擾動方程為此時擾動方程為 bxxAA )( 因因Ax=b,故有故

57、有 xAxAA)()( 因因1 A存在存在,)(1AAIAAA 若假定若假定 1|11 AAAA 則由定理則由定理20可知可知 )(1AAI 非奇異非奇異,并有并有:|11|)(|111AAAAI (5.6)由由(5.6)可得可得 ,)()(111xAAAAIx 因此因此|)(|1|11AAxAAx .|)(1|)(|)(|1|11AAACondAAACondAAAAAAAAxx (5.7)定理定理22 設(shè)設(shè)A為非奇異矩陣為非奇異矩陣,Ax=b0,且且bxxAA )( 如果如果, 1|1 AA 則則(5.7)式成立式成立.從從(5.7)看到，當(dāng)看到，當(dāng)A的條件數(shù)的條件數(shù)Cond(A)很大時，解

58、的相對誤差很大時，解的相對誤差 |xx 也很大，故方程組為病態(tài)也很大，故方程組為病態(tài).在例在例9中中 ,0001. 8| A而而5 .60000| ,1000010000300005 .3000011 AA于是于是480010|)(1 AAACond條件數(shù)很大，故方程是嚴(yán)重病態(tài)的條件數(shù)很大，故方程是嚴(yán)重病態(tài)的. 例例10 Hibert矩陣是一個著名的病態(tài)矩陣，記作矩陣是一個著名的病態(tài)矩陣，記作 1211111131211211nnnnnH它是一個對稱正定矩陣，當(dāng)它是一個對稱正定矩陣，當(dāng)n3時它是病態(tài)矩陣時它是病態(tài)矩陣.例如例如 ,611| , 33 Hn408| ,18018030180192

59、36303691313 HH故故.748)(3 HCond另外還有另外還有 764410*9 . 2)(,10*84. 2)( HCondHCond等等等等.因此因此Hn是嚴(yán)重的病態(tài)矩陣，且是嚴(yán)重的病態(tài)矩陣，且n越大越大 Cond(Hn)越大越大.例例11 在例10的方程組中可算出A的特征值 ,01015. 0,8431. 0,858. 3,2887.304321 故故2984)(412 ACond例中例中,)1 . 2, 5 . 3 , 6 .13, 2 . 8(,)31,33,23,32(,)1 . 0, 1 . 0 , 1 . 0, 1 . 0(TTTxbb 實(shí)際相對誤差是實(shí)際相對誤差是

60、 1985. 8239695.16|22 xx 而根據(jù)而根據(jù)(5.6)的誤差估計(jì)為的誤差估計(jì)為 943. 9|)(|22222 bbACondxx 這與實(shí)際相差不大，即相對誤差放大了將近這與實(shí)際相差不大，即相對誤差放大了將近3 000倍倍.故方程故方程為病態(tài)方程組為病態(tài)方程組. 定理定理23(事后誤差估計(jì)事后誤差估計(jì))設(shè)方程組設(shè)方程組 , 0, bbAx ，則，則若實(shí)際求得解為若實(shí)際求得解為x.|)(|brACondxxx 證明記剩余證明記剩余 , 0 xAbr則則,)(1rAxxrxxA .|)(|,|1|,|11bxAbACondrbAAxxxbAxrAxx 故故有有又又它表明如果方程組病