線性方程組的解法

1、第2章 線性方程組的解法 直接法與迭代法 在自然科學(xué)與社會科學(xué)的研究中，常常需要求解線性代在自然科學(xué)與社會科學(xué)的研究中，常常需要求解線性代數(shù)方程組，這些方程組的系數(shù)矩陣大致分為兩種：一種是低數(shù)方程組，這些方程組的系數(shù)矩陣大致分為兩種：一種是低階稠密矩陣（例如：階數(shù)大約為小于等于階稠密矩陣（例如：階數(shù)大約為小于等于150），另一種是），另一種是大型稀疏矩陣（即矩陣階數(shù)高且零元素較多）。大型稀疏矩陣（即矩陣階數(shù)高且零元素較多）。 在計(jì)算機(jī)上求解線性代數(shù)方程組在計(jì)算機(jī)上求解線性代數(shù)方程組 Ax=b 的常用的數(shù)的常用的數(shù)值解法：值解法：n 1、直接法：、直接法：就是經(jīng)過有限次算術(shù)運(yùn)算，可求得方程就是經(jīng)

2、過有限次算術(shù)運(yùn)算，可求得方程組精確解的方法（若計(jì)算過程中沒有舍入誤差）。但實(shí)際組精確解的方法（若計(jì)算過程中沒有舍入誤差）。但實(shí)際計(jì)算中由于舍入誤差的存在和影響，這種方法也只能求得計(jì)算中由于舍入誤差的存在和影響，這種方法也只能求得線性方程組的近似解。線性方程組的近似解。 這類方法是解低階稠密矩陣及大型帶狀方程組的有效這類方法是解低階稠密矩陣及大型帶狀方程組的有效方法。方法。n 2、迭代法：、迭代法：就是用某種極限過程去逐步逼近線性方程組就是用某種極限過程去逐步逼近線性方程組精確解的方法，迭代法具有需要計(jì)算機(jī)的存貯單元較少，程精確解的方法，迭代法具有需要計(jì)算機(jī)的存貯單元較少，程序設(shè)計(jì)簡單，原始系數(shù)

3、矩陣在計(jì)算機(jī)中始終不變等優(yōu)點(diǎn)，但序設(shè)計(jì)簡單，原始系數(shù)矩陣在計(jì)算機(jī)中始終不變等優(yōu)點(diǎn)，但存在收斂性和收斂速度等問題。存在收斂性和收斂速度等問題。 迭代法是解大型稀疏矩陣方程組的重要方法。迭代法是解大型稀疏矩陣方程組的重要方法。 3/ )3(23nnn2.1 Gauss消去法消去法2.2 直接三角分解法直接三角分解法2.3 矩陣的條件數(shù)與病態(tài)方程組矩陣的條件數(shù)與病態(tài)方程組求解線性方程組的直接方法求解線性方程組的直接方法 2.4 迭代法迭代法求解線性方程組的迭代方法求解線性方程組的迭代方法n元線性方程組11 11221121 1222221 122nnnnnnnnnna xa xa xba xa xa

4、 xba xa xa xb(2.1)方程組的矩陣形式(2.2)1111121212222212nnnnnnnnxbaaaaaaxbaaaxb簡記為.Ax b設(shè)設(shè)A非奇異，則方程組有唯一解。非奇異，則方程組有唯一解。2.1 Gauss消去法消去法n高斯消去法步驟高斯消去法步驟：(1) 首先將增廣陣 A, b 化為上三角陣；(2) 用三角方程組，回代求解 。 Gauss消去法是一個(gè)古老的求解線性代數(shù)方消去法是一個(gè)古老的求解線性代數(shù)方程組的方法（早在公元前程組的方法（早在公元前250年我國就掌握了解年我國就掌握了解三元一次聯(lián)立方程組的方法）。但由于它改進(jìn)、三元一次聯(lián)立方程組的方法）。但由于它改進(jìn)、變

5、形得到的主元素消去法，三角分解法仍然是目變形得到的主元素消去法，三角分解法仍然是目前計(jì)算機(jī)上常用的有效方法。前計(jì)算機(jī)上常用的有效方法。例例1 1 用消去法解方程組12323123 645221xxxxxxxx(1)(2)(3)解解 (1) 化上三角方程組12323123 645221xxxxxxxx 1232323 645411xxxxxxx +(-2) + 123233 64526xxxxxx 123233 64526xxxxxx (2)回代過程. 得到下同解方程組后,如下處理3( 6)/( 2)3x 從下向上逐步求解從下向上逐步求解把x3的值代入求x2用x3, x2的值求x1235 / 4

6、2xx1236()1xxx對應(yīng)的增廣矩陣的變化1116(| )04152211A b1116041504111111604150026(-2)r1 + r3r3r2 + r3 r3高斯消去法的基本思想高斯消去法的基本思想n將原方程組逐次消去未知元， 變?yōu)榕c之同解的上三角方程組， 再回代求解 。n用矩陣語言敘述，上述過程是使用初等行變換把增廣陣約化為上三角陣，使用上三角方程組，回代求解 。 高斯消去法的一般性敘述高斯消去法的一般性敘述11121121222212nnnnnnnaaabaaabaaab用行變換用行變換根據(jù)下面的上三角方程組,逐次回代求解 xk(1)(1)(1)(1)11111221

7、(2)(2)(2)22122(1)(1)(1)1,11,1( )( ) nnnnnnnnnnnnnnnnnnnnba xa xa xbaxaxaxaxbaxb(2.5)11nnxxx)()()2(2)2(2)2(22)1(1)1(1)1(12)1(11000nnnnnnnbabaabaaa22n-1,n-12.1.1 順序順序Gauss消去法消去法n在使用高斯消去法的過程中，僅對方程組做在使用高斯消去法的過程中，僅對方程組做倍倍加加變換，就形成了變換，就形成了順序高斯消去法。順序高斯消去法。1. 消元過程消元過程對于對于k =1,2,1,2,n -1-1執(zhí)行執(zhí)行(1)(1)如果如果 , ,則算

8、法失效則算法失效, ,停止計(jì)算停止計(jì)算; ;否則轉(zhuǎn)否則轉(zhuǎn)(2)(2)(2)(2)對于對于 i = k + +1, ,k + +2, , ,n 計(jì)算計(jì)算 ( )0kkka(1)(1)( ,1,2,., );(1,2,., ).ijijiiaai jnbbinn算法算法如下：記如下：記( )( )(1)( )( )(1)( )( )/(1,2,., )kkikikkkkkkijijikkjkkkiiikkmaaaam ajkknbbm b2. 回代過程回代過程( )( )( )( )( )1/(1,2,.,1)nnnnnnnkkkkkkjjkkj kxbaxbaxaknn 順序高斯消去法求解順序高

9、斯消去法求解n 元線性方程組的乘除運(yùn)算總次數(shù)為：元線性方程組的乘除運(yùn)算總次數(shù)為：32(3) /3nnn順序高斯消去法計(jì)算過程中出現(xiàn)的順序高斯消去法計(jì)算過程中出現(xiàn)的 稱為稱為主元素主元素. 出現(xiàn)出現(xiàn) ，消元過程就進(jìn)行不下去了。，消元過程就進(jìn)行不下去了。( )kkka( )0kkka即使即使 det A 0，也可能對某個(gè)也可能對某個(gè) k a=0.0120 0.010 0.1670 ;1 0.8334 5.910;3200 1200 4.2a = 1.0e+003 * 0.0000 0.0000 0.0002 0.0010 0.0008 0.0059 3.2000 1.2000 0.0042 b=0

10、.6781 12.1 981b = 0.6781 12.1000 981.0000 abans = 17.4593 -45.7600 5.5460 解：列選主元高斯消去法解：列選主元高斯消去法 1To MATLAB MV=0.0120 0.010 0.1670 0.6781;1.000 0.8334 5.910 12.10;3200 1200 4.200 981.0MV = 1.0e+003 * 0.0000 0.0000 0.0002 0.0007 0.0010 0.0008 0.0059 0.0121 3.2000 1.2000 0.0042 0.9810 y=LieZhuYuan(MV)

11、 JiaoHuan 1 hang Yu 3 hang. hang cheng shu mik(2,1 )=a(2,1)/a(1,1)=0.000313 hang cheng shu mik(3,1 )=a(3,1)/a(1,1)=0.000004 解：列選主元高斯消去法解：列選主元高斯消去法 2MV = 1.0e+003 * 3.2000 1.2000 0.0042 0.9810 0 0.0005 0.0059 0.0118 0 0.0000 0.0002 0.0007 JiaoHuan 2 hang Yu 2 hang. hang cheng shu mik(3,2 )=a(3,2)/a(2

12、,2)=0.011998 MV = 1.0e+003 * 3.2000 1.2000 0.0042 0.9810 0 0.0005 0.0059 0.0118 0 0 0.0001 0.0005y = 17.4593 -45.7600 5.5460 不選主元高斯消去法不選主元高斯消去法 Ab=0.0120 0.010 0.1670 0.6781;1.000 0.8334 5.910 12.10;3200 1200 4.200 981.0Ab = 1.0e+003 * 0.0000 0.0000 0.0002 0.0007 0.0010 0.0008 0.0059 0.0121 3.2000 1

13、.2000 0.0042 0.9810 y=ShunXuGauss(Ab) mik=aik/akk=83.333333 mik=aik/akk=266666.666667 Ab = 1.0e+005 * 0.0000 0.0000 0.0000 0.0000 0 0.0000 -0.0001 -0.0004 -0.0000 -0.0147 -0.4453 -1.7985 mik=aik/akk=-14670000 Ab = 1.0e+008 * 0.0000 0.0000 0.0000 0.0000 0 0.0000 -0.0000 -0.0000 -0.0000 0 -1.7619 -6.5

14、17y = -104.0 100.0 5.5460 2.2.1 Doolittle分解法與分解法與Crout分解法分解法111ijlL=單位下三角陣單位下三角陣記記 (1)(1)(1)11121(2)(2)222( ).0.00nnnnnaaaaaUaLUA 上三角陣上三角陣三角分解三角分解三角分解與解方程組三角分解與解方程組2.2 直接三角分解法直接三角分解法定義定義 A = LU 叫做叫做 A 的的三角分解三角分解。 L，U 是是下、上三角陣下、上三角陣.若若 L 是單位下三角矩陣，是單位下三角矩陣， U 是上三角矩陣，則是上三角矩陣，則 A =LU 叫叫 A 的的Doolittle 分解

15、分解；若若 L 是下三角矩陣，是下三角矩陣，U 是單位上三角矩陣，是單位上三角矩陣， A =LU 叫叫 A 的的 Crout 分解分解。三角分解用處三角分解用處 如果方程組如果方程組 Ax =b 的系數(shù)陣的系數(shù)陣 A 能分解為能分解為A =LU, 其中，其中，L 是下三角矩陣，是下三角矩陣，U 是上三角矩陣是上三角矩陣. 這時(shí)解方程組這時(shí)解方程組 Ax=b， 可化為求解兩個(gè)三角方程組可化為求解兩個(gè)三角方程組Ly =b, Ux =y . 先由先由 Ly =b 解出向量解出向量 y，再由再由 Ux =y 解出向量解出向量 x， 則則 x 是是原方程組原方程組 Ax =b 的解向量。的解向量。111

16、212122212100100100nnnnnnuuuluuAllu 矩陣矩陣 ARnn 的的 Doolittle 分解分解下面，考察第下面，考察第 k 次順序消元時(shí)，次順序消元時(shí)，Ax =b 增廣陣的變化情況增廣陣的變化情況. .單位下三角單位下三角陣陣上三角陣上三角陣22,2,kkkkkrmrr,nn kknrmrr11,1,kkkkkrmrr第第 k 次消元使用的行變換次消元使用的行變換:(1)(1)(1)(1)11121(1)(1)(1)1,11,1( )( )( )( )( ),( )( )( ),00,0000nkkkkkknkkkkkkk kk nkkkkn kn nnaaaba

17、abAbaabaab 第第k個(gè)增廣陣個(gè)增廣陣,22,11,(),(),()()()()nn kkkkkkkkkkkrmrrmrrmrEEEP做這做這 n- -k 次行變換相當(dāng)與給次行變換相當(dāng)與給 A(k),b(k) 乘一些乘一些倍加倍加初等方陣初等方陣,( )( )/.kkikikkkmaa其中其中, ,Pk定義定義 下面的 nn 矩陣叫做初等下三角矩陣初等下三角矩陣。 1,1000010(1,2,.,1)0010000001kkkn kPknpp 由上面的討論，可引出下面的定義由上面的討論，可引出下面的定義。 其中，pik 是行乘數(shù) -mik .1kP1.11, 1knkkmm 111121

18、.nPPP111jim,記為記為L初等下三角陣初等下三角陣的逆陣之積為的逆陣之積為單位下三角陣單位下三角陣初等下三角陣初等下三角陣的逆陣為的逆陣為初等下三角陣初等下三角陣定理定理2.3 矩陣矩陣 Ann (n2) 有唯一的有唯一的 Doolittle 分解充要分解充要條件是條件是 A 的前的前n-1個(gè)順序主子式個(gè)順序主子式 Dk0 ( k = 1,2, n-1)。 證明證明 充分性充分性 因?yàn)橐驗(yàn)?Dk0 ( k = 1,2, n-1) ，依定理，依定理2.1，可對矩陣進(jìn)行順序高斯消去法，把可對矩陣進(jìn)行順序高斯消去法，把 A 化為上三角矩陣化為上三角矩陣 U，其變換過程相當(dāng)于做矩陣乘法其變換過

19、程相當(dāng)于做矩陣乘法Pn- -1P2 P1 A = U (2.9)其中其中, Pk 為初等下三角陣為初等下三角陣, Pk中中pik 的是行乘數(shù)的是行乘數(shù)- -mik ; U 的主對的主對角線元素角線元素 . ( )0, (1,2,1)kkkkkuakn若若 A 非奇異，則非奇異，則 unn 0 ; 若若 A 奇異，則奇異，則unn =0 。 由式由式(2.9) 得得 .111121nAP PP U記記 ，則，則 L 是一個(gè)單位下三角陣，于是有是一個(gè)單位下三角陣，于是有A=LU.111121nLP PP 唯一性唯一性 設(shè)有分解設(shè)有分解.ALULU當(dāng)當(dāng) A 非奇異時(shí)，非奇異時(shí)，L , L* ; ;U

20、, U* 均非奇異，從上式可得均非奇異，從上式可得11().U UL L 由于由于 L-1L* 是單位下三角矩陣，是單位下三角矩陣，UU*-1 是單位上三角陣是單位上三角陣，且二者相等且二者相等, ,它們只能為單位陣。即它們只能為單位陣。即, ,L-1L* = I， UU*-1 = I . 由此得到由此得到 L = L* ， U=U* 。 當(dāng)當(dāng) A 奇異時(shí)，唯一性也成立。證明見教材奇異時(shí)，唯一性也成立。證明見教材 P20。 必要性必要性 設(shè)設(shè) A 有唯一的有唯一的 D-分解分解 A = LU ，則必有，則必有 uii 0 ( i= 1, 2, n- -1) ; 否則存在否則存在 ukk = 0

21、 , 而而 u11 , u22 , uk-1, k-1 均均不為零。這時(shí)，前不為零。這時(shí)，前 k + 1 階主子矩陣階主子矩陣11,11,1001kkkkTTkkkkAyUsLAxaur其中其中Ak , Lk , Uk 分別是分別是 A, L, U 的的 k 階順序主子矩陣階順序主子矩陣. 可知，可知，,TTTkkxr UxU r 因因 Uk 奇異，故奇異，故 r 不存在或不唯一。這與不存在或不唯一。這與A 有唯一有唯一D-分解矛盾分解矛盾.由于由于 uii 0 ( i= 1, 2, n- -1) 及及 Ak = LkUk ，可知可知Dk = detAk = u11u22 ukk 0, ( k

22、 = 1, 2, n- -1) .推論推論2.3 矩陣矩陣 Ann (n2) 有唯一的有唯一的 Crout 分解的充要條分解的充要條件是件是 A 的前的前 n -1 個(gè)順序主子式個(gè)順序主子式 Dk0 ( k = 1,2, n-1)。 證證 僅須證明僅須證明 A 有唯一的有唯一的 Crout 分解等價(jià)于分解等價(jià)于A 有唯一有唯一的的 Doolittle 分解分解. . 其中其中, D = diag(u11,u22,unn), 且且 D 與與 都是唯一的都是唯一的. 于是有于是有其中其中 是下三角陣是下三角陣, 也是唯一的也是唯一的. UALULDULU LLDL A 有唯一的有唯一的 D-分解分

23、解 A=LU ，則，則 U 的前的前 n- -1 個(gè)主對角元個(gè)主對角元 uii 0 ( i = 1, 2, n- -1) ; 這時(shí)這時(shí), U 可分解為對角陣可分解為對角陣 D 與單位上三角陣與單位上三角陣 之積之積: .UDUU反之反之, 同理可證同理可證.求解線性方程組的求解線性方程組的Doolittle分解法分解法111212122212,nnnnnnaaaaaaAaaa設(shè)有設(shè)有 Doolittle 分解分解 A = LU ,11121212221210010.0100nnnnnnuuuluuLUllu下三角陣下三角陣 L 與上三角陣與上三角陣 U 的求解的求解由矩陣乘法由矩陣乘法 及及

24、矩陣相等定義矩陣相等定義，知111 111(1,2,., )(2,3,., )jjiiuajnl uain第1步，解出U 的第1行與L 的第1 列元！(2.11)當(dāng)當(dāng) k k = 2,3,n 2,3,n 時(shí)，有時(shí)，有111111(,1,., )(1,2,., ;)nkk jktt jktt jk jttnkikittkittkikkkttal ul uujk knal ul ul uikkn knlkt = 0, t k utk = 0, t k解之，有解之，有1111 (,1,., )()/ (1,2,., ;)kk jk jktt jtkikikittkkktual ujk knlal u

25、uikkn kn第2步，解出U 與L 的其它元！(2.12)公式公式 (2.11) 與與 (2.12)，求解過程圖示如下：，求解過程圖示如下：11, (1,2,., )jjuajn1111/, (2,3,., )iilauin 11kk jk jktt jtual u, jk 11()/,kikikittkkktlal uuik 11kk jk jktt jtual u, jk121,1,1,jjkjkjkk kk kuuuullu?上三角陣上三角陣 U 中元素的求法中元素的求法與與ukj同行同列的數(shù)被同行同列的數(shù)被用到用到! !與與ukj同行同列的數(shù)被同行同列的數(shù)被用到用到! !11,1,1

26、212,1,kkkk kk kkkkjiii ki kuulullullll? 11()/,kikikittkkktlal uuik下三角陣下三角陣 L 中元素的求法中元素的求法下、上三角方程組下、上三角方程組 Ly = b 與與 Ux = y 的計(jì)算公式的計(jì)算公式用同樣的方法可得用同樣的方法可得 Crout 分解的計(jì)算公式：分解的計(jì)算公式：1/()/ (1,2,.,1)nnnnniiittiit ixyuxyu xuinn 1111 (2,3,., )iiiitttybybl yin(2.13)左行乘右列求u12左行乘右列求u1n Crout 分解分解111212122212nnnnnnaa

27、aaaaAaaa11121212221200101.0001nnnnnnluullulll下三角陣下三角陣 L 與單位上三角陣與單位上三角陣 U 的求解的求解對對 k =1,2,n，均有，均有1111 (,1,., )()/ (1,2,., ;)kikikittktkkjkjkttjkktlal uik knual uljkkn kn(2.14)左行乘右列求l11下、上三角方程組下、上三角方程組 Ly=b 與與 Ux=y 的計(jì)算公式的計(jì)算公式1111111/()/ (2,3,., ) (1,2,.,1)iiiittiitnnniiittt iyblybl ylinxyxyu xinn (2.1

28、5)例例 用矩陣的杜利脫爾（用矩陣的杜利脫爾（Doolittle）分解解方程組）分解解方程組解：解：設(shè)設(shè) .201814513252321321xxxLUuuuuuulll332322131211323121111513252321比較兩邊系數(shù)得比較兩邊系數(shù)得:2454132321333223223121131211uluulluuu2441321153121UL于是201814153121321yyyLy解方程組,72,10,14321yyy解得7210142441321321xxxUx再解方程組. 1, 2, 3123xxx解得練習(xí)：練習(xí)： 用矩陣的杜利脫爾（用矩陣的杜利脫爾（Doolit

29、tle）分解：）分解：A=LU 解方程組解方程組答案：答案：731395222211212032114321xxxx0, 1, 1, 34321yyyy1, 0, 1, 01234xxxx例例2 在四位十進(jìn)制的限制下，用在四位十進(jìn)制的限制下，用 Doolittle 分解法分解法 求解下列方程組求解下列方程組1231231238.1 2.31.56.10.56.230.872.32.51.510.21.8xxxxxxxxx解解 此方程組的系數(shù)陣此方程組的系數(shù)陣 A 的各階順序主子式不為零，有的各階順序主子式不為零，有 D-分解分解 A =LU，且，且 U 非奇異，用公式非奇異，用公式(2.12)

30、，得分解結(jié)果，得分解結(jié)果 1112132122233132338.12.31.50.061736.3720.96260.30860.12410.78uuuluullu再用公式(2.13), 可得結(jié)果 1231236.1, 1.923, 0.1560.01447,0.2996,0.8408yyyxxx 此方程組精確解舍入到四位有效數(shù)字是1230.01445,0.2997,0.8409xxx n-1-1個(gè)順序主子式不為零？好苛刻的個(gè)順序主子式不為零？好苛刻的條件吆條件吆! ! 判別判別, ,還要浪費(fèi)計(jì)算時(shí)間！還要浪費(fèi)計(jì)算時(shí)間！2.2.2 選主元的選主元的 Doolittle 分解法分解法矩陣矩陣A

31、有有 D-分解的條件是什么分解的條件是什么?還有一種選主元的還有一種選主元的D- -分解！分解！所要求的條件弱一些！所要求的條件弱一些！前前n-1 個(gè)順序主子式不為零！個(gè)順序主子式不為零！定義定義 對 n 階單位陣 In 做一次交換 2 行的操作,得到的方陣叫初等置換陣初等置換陣； 對 In 做多次交換 2 行的操作，得到的方陣 Q 叫 置換陣置換陣. 置換陣 Q 可以表為一些初等置換陣的乘積. 定理定理2.4 若矩陣若矩陣 A 非奇異，則非奇異，則 存在置換陣存在置換陣 Q，使，使 QA 可可做做 Doolittle 分解分解 QA = LU 。 證證 若若 A 非奇異非奇異, 依定理依定理

32、 2.2，可對，可對 A 進(jìn)行列主元消去，進(jìn)行列主元消去，把把 A變?yōu)樯先顷囎優(yōu)樯先顷?U。變換過程相當(dāng)于矩陣運(yùn)算。變換過程相當(dāng)于矩陣運(yùn)算Pn-1 Qn-1 P2 Q2 P1 Q1 A = U. 其中，其中，Pk , Qk 分別是分別是初等下三角陣初等下三角陣 和和 初等置換陣初等置換陣。kkknkkikQikQQ1101111011.,11其中置換陣置換陣: 初等置換陣初等置換陣:因?yàn)?，因?yàn)?，QkQk = I， 故，故，Pn-1 Qn-1P2 Q2 P1(Q2 Qn-1)(Qn-1 Q2 ) Q1 A = U.其中，其中，L-1 = Pn-1 Qn-1 P2 Q2 P1 (Q2 Qn-1

33、) 是單位下三角陣。是單位下三角陣。以以 n = 4 說明，說明，L-1 = P3 Q3 P2 Q2 P1(Q2 Q3) =P3(Q3 P2 Q3 ) (Q3 Q2 P1 Q2 Q3)是單位下三角陣，是單位下三角陣， Q = Qn-1 Q2 Q1 是置換陣。故是置換陣。故QA = LU 。選主元選主元 D-分解的實(shí)施方法分解的實(shí)施方法 采取與列主元類似的方法，在 D -分解中也選主元素。設(shè)使用 D-分解公式 (2.12) 已進(jìn)行了 k -1 步，第 k -1 次分解矩陣為1111212222(1)1,11,21,11,1,1,11,1knknkkkkkkkkknkk kkkknnn knknn

34、uuuluuulluuuAllaallaa進(jìn)行第 k 步分解，(1) 先尋找主元素。計(jì)算中間量11(,1,., )kiikittktsal uik kn滿足 的 就是第 k 步的主元素，以 的值作為 。為此，交換矩陣 A(k-1) 的第 k 行與第 ik 行。 得到 QA 。| max |kiik i nss kiskiskku(2) 再按公式，對 QA 進(jìn)行第 k 步的D-分解運(yùn)算。這時(shí)，原方程組變?yōu)?LUx = Qb, 它可歸結(jié)為求解下面兩個(gè)方程組：Ly = Qb, Ux = y選主元選主元 D-分解的算法分解的算法(1) 做分解做分解 QA = LU . 定義整型數(shù)組 M(n), 其第

35、k 個(gè)元素 Mk 用于記錄第 k 個(gè)主元所在的行號。 對于對于 k =1,2,n 執(zhí)行執(zhí)行 . (1.1) 計(jì)算中間量計(jì)算中間量11(,1,., )kiikittktsal uik kn(1.2) 選行號選行號 i k 使使 ，令，令 Mk = i k .| max |kiik i nss (1.3) 若若i k=k,則轉(zhuǎn)則轉(zhuǎn)(1.4)；否則；否則, 交換交換 與與 與與 及及 與與 中的數(shù)據(jù)，中的數(shù)據(jù)，轉(zhuǎn)轉(zhuǎn) (1.4) 。,(1,2,.,1)ki tltk,k tl,k ta,(1,2,.,1)ki tatkkskis (1.4) 計(jì)算計(jì)算11(1,2,., ;)/(1,2,., ;)kkk

36、kkjkjkttjtikikkusual ujkkn knlsuikkn kn(2) 求求 Qb.對于對于 k = 1,2,n-1 執(zhí)行執(zhí)行 .(2.1) 令令 t = Mk.(2.2) 交換交換 bk 與與 bt 所含數(shù)據(jù)。所含數(shù)據(jù)。 矩陣矩陣 L，U 中的元素中的元素1(1,2,.,2,1)/()/nnnnniiittiit iinnxyuxyu xu 1111(2,3,., )iiiktttybybl yin(3) 求解方程組求解方程組 Ly = Qb 和和 Ux = y .例：例：用選主元的杜利脫爾（用選主元的杜利脫爾（Doolittle）分解解方程組）分解解方程組.432121111

37、011321xxx解：解：設(shè)設(shè) 1111Q1010110011P)1(11110100011AAQP11101010000122PQ設(shè)UAAQP:100110011)2() 1 (22于是AQQQPQPAQQQPQPAQPQPU)(1221221221221122即21221221212)(QPQQPQPLQQQQ令1010110012121QPQLQbLy 解方程組, 1, 2, 2321yyy解得yUx 再解方程組. 1, 1, 1123xxx解得n定義定義 設(shè)矩陣設(shè)矩陣 A = aij，如果存在兩個(gè)正整數(shù)，如果存在兩個(gè)正整數(shù) r 和和 s，使得，使得當(dāng)當(dāng) i - -j r 以及以及 j

38、- -i s 時(shí)，有時(shí)，有 aij=0 ， 則則 稱稱 A 是上半帶寬為是上半帶寬為 s，下半帶寬為，下半帶寬為 r 的的帶狀矩陣帶狀矩陣，線性，線性方程組方程組 Ax =b 稱為稱為帶狀方程組帶狀方程組。即，帶狀矩陣。即，帶狀矩陣 A 為為2.2.3 三角分解法解帶狀線性方程組三角分解法解帶狀線性方程組111,11,1,srn s nn n rnnaaaAaaa定理定理2.5 (保帶狀三角分解保帶狀三角分解) 設(shè)矩陣設(shè)矩陣 A 滿足下面兩個(gè)條件：滿足下面兩個(gè)條件：(1) A = aij 是上半帶寬為是上半帶寬為 s，下半帶寬為，下半帶寬為 r 的帶狀矩陣；的帶狀矩陣；(2) A 的的前前 n

39、 -1個(gè)順序主子式個(gè)順序主子式均不為零。均不為零。 則則 A 有唯一的有唯一的 D-分解分解 A=LU，其中，其中，L 是下半帶寬為是下半帶寬為 r 的的單位三角陣單位三角陣，U 是上半帶寬為是上半帶寬為 s 的上三角陣，即，的上三角陣，即，11121,121,1,11,1100100000100sn s nrnnn n rn nnnuuuluAlullu 證證 條件 (2) 說明 A 有 D-分解： A = LU.當(dāng)當(dāng) i - -j r 時(shí)時(shí), 由條件 (1) 知 aik = 0 ( k = 1, 2, j ), 再由分解公式11()/,kikikittkkktlal uu故，L 是下半帶寬

40、為 r 的單位下三角陣. 1111221 12223/,()/,0,0.iiiiiiijlaulal uull可推知11222130,0,0,0.jjjjjijuaualuu可推知 故，U 是上半帶寬為 s 的下三角陣. 當(dāng)當(dāng) j - - i s 時(shí)時(shí), 由條件(1)知 akj = 0 ( k = 1, 2, i ), 再由分解公式11.kk jk jktt jtual u推論推論2.5 設(shè)矩陣設(shè)矩陣 A = aij 滿足定理滿足定理 2.5 的條件，則的條件，則 A 有唯有唯一的一的 Crout 分解分解 A =LU，其中，其中，L 是下半帶寬為是下半帶寬為 r 的三角陣的三角陣，U 是上半

41、帶寬為是上半帶寬為 s 的單位上三角陣。的單位上三角陣。帶狀分解算法：帶狀分解算法： (1) 做分解做分解 A =LU。 對對 k = 1, 2, n 計(jì)算計(jì)算 1max1,1max1,(,1,.,min, )(1,2,.,min, ;), ()/,kk jk jktt jtk r j skikikittkkkti r k sjk kks nikkkr nknual ulal uu(2) 求解求解 Ly = b, Ux = y 。max, 1(1,2,.,2,1)/()/nnnnti s niiittiit iinnxyuxyu xu 111max1,(2,3,., )iiikttti rin

42、ybybl y1122231nnnacdacdAcda 2.2.4 追趕法求解三對角線性方程組追趕法求解三對角線性方程組設(shè) n 元線性方程組 Ax = b 的系數(shù)矩陣為非奇異的 三對角陣三對角陣此方程組叫做三對角方程組三對角方程組。 A 是上、下半帶寬為1的帶狀矩陣。依推論依推論2.5，有唯一的，有唯一的 Crout 分解分解.由矩陣乘法及矩陣相等定義由矩陣乘法及矩陣相等定義，有下面的等式：1111 1,apcp q對于對于 i = 1,2,n-1 有有11111,iiiiiiiiidrar qpcp q1122231001000100000001nnnpqrpqrAqrp (2.18)算法算

43、法(1) Crout 分解計(jì)算分解計(jì)算11,pa對于對于 i = 1,2,n-1 有有(2) 求解方程組求解方程組 L y = b 與與 Ux = y 追追求yi123nyyyy111/iiiiiiiqcppadq(2.19)1111(2,3,., )/()/iiiiiinybpybd yp(2.20)追追趕趕求xi1(1,2,.,1)nniiiixyxyq xinn(2.20)注解注解 追趕法的存貯與計(jì)算量都很小，乘除追趕法的存貯與計(jì)算量都很小，乘除運(yùn)算總次數(shù)為運(yùn)算總次數(shù)為 5n-4。例例3 在四位精度下，在四位精度下， 求解下列方程組。求解下列方程組。123454111410.514111

44、413142xxxxx 解解 用公式用公式 (2.19) 可有系數(shù)陣的分解式可有系數(shù)陣的分解式 A =LU ，410.2513.7510.266713.73310.267913.73210.267813.7321用公式用公式 (2.20) 可求出下面的解可求出下面的解123450.25,0.06667,0.2858,0.8805,0.3001yyyyy 此方程組的精確解為此方程組的精確解為543210.3301,0.8001,0.5001,0.2001,0.2000 xxxxx 543210.3,0.8,0.5,0.2,0.2xxxxx 例：例：試用試用“追趕法追趕法”求解線性代數(shù)方程組：求解

45、線性代數(shù)方程組： 022112111131124321xxxx解：解：設(shè)設(shè) ULqqqprprprp111112001110013100123214433221比較兩邊系數(shù)得比較兩邊系數(shù)得:37235531522512124433322211prqprqprqp13515212113725312512UL于是022137253125124321yyyy解方程組2,37,53,214321yyyy解得237532113510521002114321xxxx再解方程組. 0, 11, 21234xxxx解得.200031002310023100224321xxxx練習(xí)練習(xí) 試用試用“追趕法追趕法”

46、求解線性代數(shù)方程組：求解線性代數(shù)方程組： A0=2 3 3 3;Ad=0 1 1 1 ;Ac=2 2 2 0;n=length(A0);P =zeros(1,n); R=Ad; Q=P;P(1)=A0(1); for i=1:(n-1) Q(i)=Ac(i)/P(i); P(i+1)=A0(i+1)-Ad(i+1)*Q(i); fprintf(nQ%d = C%d / P%d = %f n,i,i,i, Q(i) fprintf(nP(%d+1) = a(%d+1)-d(%d+1)Q(%d) = %fn ,i,i,i,i, P(i+1) endRPQ答案答案:Q1 = C1 / P1 = 1.

47、000000 P(1+1) = a(1+1)-d(1+1)Q(1) = 2.000000 Q2 = C2 / P2 = 1.000000 P(2+1) = a(2+1)-d(2+1)Q(2) = 2.000000 Q3 = C3 / P3 = 1.000000 P(3+1) = a(3+1)-d(3+1)Q(3) = 2.000000 R = 0 1 1 1P = 2 2 2 2Q = 1 1 1 0 1, 0, 0, 04321yyyy1, 1, 1, 11234xxxx2.3 矩陣的條件數(shù)與病態(tài)線性方程組矩陣的條件數(shù)與病態(tài)線性方程組n2.3.1 矩陣的條件數(shù)與線性方程組的性態(tài)矩陣的條件數(shù)與

48、線性方程組的性態(tài) 給定線性方程組給定線性方程組 Ax =b，現(xiàn)在考察，系數(shù)矩陣，現(xiàn)在考察，系數(shù)矩陣 A 和常數(shù)列和常數(shù)列 b 有了微小變化有了微小變化 A，b ，它如何影，它如何影響解向量響解向量 x，即，解向量，即，解向量 x 的變化量的變化量 x 何樣？何樣？ 由于由于A （或（或 b）的元素是測量得到的，或者是）的元素是測量得到的，或者是計(jì)算的結(jié)果，在前種情況下，計(jì)算的結(jié)果，在前種情況下， A （或（或 b）常常帶有）常常帶有某些觀測誤差，在后種情況下，某些觀測誤差，在后種情況下， A （或（或 b）包含舍）包含舍入誤差，因此我們處理的實(shí)際矩陣是入誤差，因此我們處理的實(shí)際矩陣是A + A

49、 （或（或 b+ b ）。）。1122112112, 11.0001211.00012.0001xxxx 精確解精確解 (2, 0)精確解精確解 (1, 1) 考察方程組考察方程組 Ax = b, 當(dāng)當(dāng) A 或或 b 有微小擾動時(shí)有微小擾動時(shí), 對解的影響，對解的影響，首先看一個(gè)例子：首先看一個(gè)例子： n定義定義 如果矩陣如果矩陣 A 或常數(shù)項(xiàng)或常數(shù)項(xiàng) b 的微小的微小變化，引起方程組變化，引起方程組 Ax =b 的解的巨大變的解的巨大變化，則稱方程組為化，則稱方程組為“病態(tài)病態(tài)”方程組，矩方程組，矩陣陣 A 稱為稱為“病態(tài)病態(tài)”矩陣；否則稱方程組矩陣；否則稱方程組為為“良態(tài)良態(tài)”方程組，方程

50、組， A 稱為稱為“良態(tài)矩良態(tài)矩陣陣”。定理定理2.6 設(shè)設(shè) 非奇異，非奇異，b0, x 是方是方程組程組 Ax =0 的解向量。的解向量。 若若 則有則有,n nnAARbbRA 11AA (1) 方程組方程組（2.21） 有唯一的解向量有唯一的解向量 ()()AA xxbb .xx(2) 下列估計(jì)式成立：下列估計(jì)式成立：（2.22）111AAxAbxAbAA 證證 (1) 因 ，又111AAAA 1()AAA IAA 由定理1.5知， 非奇異，及 A 非奇異，知 非奇異，方程組 (2.21) 有唯一解 .1IAAAA AAxxbbAxA xAxb 依方程組 (2.21) 及 b =Ax ,

51、 可得 又,bAx11AA 在上面最后一式兩端，同時(shí)左乘 A-1 ，得1xAbA xAx 再估計(jì)范數(shù)，有1xAbAxAx 11xAAxAbAx 111AAxAbAx 1111AAxAbAAx ,bxAxAbbbxxb 1111bAAxAAxAAxb 11.1AAxbAxbAAA 111AAxAbAx 用不等式放用不等式放! !11.1AAxbAxbAAAAA結(jié)論 (2.22) 可寫成: 由此可以看出，量由此可以看出，量 |A| |A-1| 越小，由越小，由A （或（或 b）的相對誤差引起的解的相對誤差的相對誤差引起的解的相對誤差| x | / | x |就越就越小小 ；量；量 |A| |A-1

52、| 越大，解的相對誤差就可能愈大，越大，解的相對誤差就可能愈大，所以量所以量 |A| |A-1| 實(shí)際上刻劃了解對原始數(shù)據(jù)變化實(shí)際上刻劃了解對原始數(shù)據(jù)變化的靈敏程度，即刻劃了方程組的的靈敏程度，即刻劃了方程組的“病態(tài)病態(tài)”程度。程度。誤差分析誤差分析定義定義 對于可逆矩陣對于可逆矩陣 A，稱，稱 |A| |A-1| 為為 A 的條件數(shù)，記作的條件數(shù)，記作cond(A) = | A | | A-1 | .矩陣的條件數(shù)與所取的范數(shù)有關(guān)，常用的條件數(shù)是矩陣的條件數(shù)與所取的范數(shù)有關(guān)，常用的條件數(shù)是cond(A) = | A | | A-1 | ， cond(A)1 = | A | 1 | A-1 |

53、1， cond(A)2 = | A | 2 | A-1 | 2 .定義定義 若可逆矩陣若可逆矩陣 A 的條件數(shù)的條件數(shù) cond(A) 相對很大，則稱相對很大，則稱A是是病態(tài)矩陣，相應(yīng)的方程組病態(tài)矩陣，相應(yīng)的方程組Ax =b叫病態(tài)方程組；反之，稱二叫病態(tài)方程組；反之，稱二者是良態(tài)的。者是良態(tài)的。條件數(shù)的性質(zhì)條件數(shù)的性質(zhì) cond(A) = | A | | A-1 | n(1) 對任何可逆陣對任何可逆陣 A，cond(A) 1.n(2) A 是可逆陣是可逆陣, k 0，則，則 cond(kA) = cond(A).n(3) A 可逆實(shí)對稱陣，則可逆實(shí)對稱陣，則 cond(A)2 = | l l1 | | l ln |.其中，其中， l l1 和和l ln 分別是分別是 A 的模最大和最小的特征值。的模最大和最小的特征值。n(4) A 正交陣，則有正交陣，則有 cond(A)2 = 1。n(5)設(shè)設(shè)A為可逆陣為可逆陣 ，P為正交陣，則有為正交陣，則有 cond(PA)2 = cond(AP)2 = cond(A)2 例例4 設(shè)線性方程組 Ax = b 為121 1.00012112xx (1) 求 cond(A) 和 Ax = b 的解 x 。(2) 設(shè) b +b =(2.0001, 2)T, 求 A( x +x) = b +b 的解 x +x .(3) 計(jì)算bxbx