1、第二章第二章 解線性方程組的直接法解線性方程組的直接法第二章第二章 解線性方程組的直接法解線性方程組的直接法2.1 引言引言n小行星軌道問題：小行星軌道問題： 天文學家要確定一小行星的軌道，在軌道平面建天文學家要確定一小行星的軌道，在軌道平面建立以太陽為原點的直角坐標系。在坐標軸上取天文丈立以太陽為原點的直角坐標系。在坐標軸上取天文丈量單位量單位(一天文單位為地球到太陽的平均間隔：一天文單位為地球到太陽的平均間隔：9300萬萬英里，約英里，約1.5億千米億千米)，對小行星作，對小行星作5次察看，測得軌道次察看，測得軌道上上5個點的坐標數據如下：個點的坐標數據如下： x 5.7640 6.286

2、0 6.7590 7.1680 7.4800 y 0.6480 1.2020 1.8230 2.5260 3.3600橢圓的普通方程橢圓的普通方程: a1x2 + a2xy + a3y2 + a4x + a5y + 1 = 0將數據逐個代入，可得五個方程的方程組，求解該線性將數據逐個代入，可得五個方程的方程組，求解該線性方程組即可得行星軌道方程。方程組即可得行星軌道方程。nnnnnnaaaaaaaaaA212222111211nbbbb2112nxxxx由以前所學內容知，當且僅當矩陣由以前所學內容知，當且僅當矩陣A行列式不為行列式不為0時，即時，即A非奇特時，方程組存在獨一解，可根據非奇特時，

3、方程組存在獨一解，可根據Cramer法那么求解。法那么求解。其算法設計如下：其算法設計如下：1 輸入系數矩陣輸入系數矩陣A和右端向量和右端向量b；2計算系數矩陣計算系數矩陣A的行列式值的行列式值D，假設，假設D=0，那么輸，那么輸出錯誤信息，終了，否那么進展第出錯誤信息，終了，否那么進展第(3)步；步；3 對對k=1,2,n，用，用b交換交換A的第的第k列數據，并計算列數據，并計算交換后矩陣的行列式值交換后矩陣的行列式值Dk；4 計算并輸出計算并輸出x1 = D1 / D，x2 = D2 / D，xn=Dn/D， 終了。終了。n但但Cramer法那么只適用于低階方程組，高階方程組任法那么只適用

4、于低階方程組，高階方程組任務量太大，故普通用數值方法求解。數值方法分兩類：務量太大，故普通用數值方法求解。數值方法分兩類：n1. 直接法直接法n2. 迭代法迭代法2.2 Gauss消元法消元法u第一步第一步: 消元過程，將方程組消元化為等價的消元過程，將方程組消元化為等價的上三角形方程組上三角形方程組; )1()1()1(2)1(22)1(2211212111nnnnnnnnnnbxabxaxabxaxaxaGauss消元的目的：消元的目的：原始方程組原始方程組約化方程組約化方程組11 11221121 1222221 122nnnnnnnnnna xa xa xba xa xa xba xa

5、 xa xb消元過程消元過程(化普通方程組為上三角方程組化普通方程組為上三角方程組)4444343242141343433323213124243232221211414313212111bxaxaxaxabxaxaxaxabxaxaxaxabxaxaxaxa444434241334333231224232221114131211baaaabaaaabaaaabaaaaA)()()()()()()()()()()()()(1414414314213134133132121241231221141312111baaabaaabaaabaaaaA)0(11 a第一輪消元：第一輪消元：u計算計算3個

6、數個數: m21 m31 m41T = a21 a31 a41T / a11u用用-m21乘矩陣第一行后加到矩陣第二行乘矩陣第一行后加到矩陣第二行;u用用-m31乘矩陣第一行后加到矩陣第三行乘矩陣第一行后加到矩陣第三行;u用用-m41乘矩陣第一行后加到矩陣第四行乘矩陣第一行后加到矩陣第四行;其系數增廣矩陣變為：其系數增廣矩陣變為：)()()()()()()()()()()(2424424323234233121241231221141312112baabaabaaabaaaaA)0()1(22 a計算計算2個數：個數：m32 m42T = a32(1) a42(1)T / a22(1)u 用用

7、-m32乘矩陣第二行后加到矩陣第三行乘矩陣第二行后加到矩陣第三行;u用用-m42乘矩陣第二行后加到矩陣第四行乘矩陣第二行后加到矩陣第四行;其系數增廣矩陣變為：其系數增廣矩陣變為：)()()()()()()()()()(3434423234233121241231221141312113babaabaaabaaaaA)0()2(33 a第三輪消元：第三輪消元：u計算計算: m43=a43(2)/a33(2)u用用-m43乘矩陣第三行后加到矩陣第四行乘矩陣第三行后加到矩陣第四行;其系數增廣矩陣變為：其系數增廣矩陣變為：)()()()()()()()()(3443442342343233124124

8、312321221414313212111bxabxaxabxaxaxabxaxaxaxa(1)(1)( )(1)(1)( )(1)(1)(1, )( ,1, )(1, )kikikkkkkkkijijikkjkkkiiikkamiknaaam ai jknbbm bikn）（記ijijaa)0(回代過程回代過程(解上三角方程組解上三角方程組)nnnnnnnnbxabxaxabxaxaxa2222211212111) 1 , 2 , 1()(1niaxabxabxiinijjijiinnnnn回代過程的計算公式：回代過程的計算公式：t總任務量：總任務量：S1=n(n-1)/2+ (n3-n)/

9、3t總任務量：總任務量：S2=n+ n(n-1)/2= n(n+1)/2(1)(1)( )(1)(1)( )(1)(1)(1, )( ,1, )(1, )kikikkkkkkkijijikkjkkkiiikkamiknaaam ai jknbbm bikn1(),(1,2,1)niijjj inninniiba xbxxinaa 假設用克萊姆法那么求解，那么任務量為：假設用克萊姆法那么求解，那么任務量為：“：n+1個個n階行列式的值階行列式的值(n+1)(n-1)n!“：n故總任務量為：故總任務量為： (n+1)(n-1) n!+n如當如當n=6時時, Gauss消元法任務量為消元法任務量為1

10、06 ；而克萊姆法；而克萊姆法那么求解任務量為那么求解任務量為25206。01111 kkkkkaaaaD111( )1,11/(2,3,.,1)kkkkkaDaDDkn 32563034253432321xxx 32303452536432 202230445 . 006432 4200445 . 006432 446245 . 0432321xxxx1= -13, x2 = 8, x3 = 2m21=3/2m31=4/2m32=-3/0.5 1111211nmmL 11111,1,knkkkmmLn將將A 分解為單位下三角矩陣分解為單位下三角矩陣 L 和上三角矩陣和上三角矩陣 U的乘積的算

11、法稱為矩陣的乘積的算法稱為矩陣A的三角分解算法。的三角分解算法。)1(1221 nnnAALLLL)1(1221 nnnbbLLLLLUULLLAn 111211 1112121nnmmmL2131321231111nnnmLmmmmm由由Gauss消元過程可推得消元過程可推得U即為即為Gauss消元后所得的上三角方程的系數矩陣。消元后所得的上三角方程的系數矩陣。111041221Am21=0， m31=2， m32= -1故故112010001200140111A=LU假設曾經有假設曾經有 A = L U 那么那么 AX = b = L U X = b，1求解方程組求解方程組 LY = b

12、得向量得向量 Y 的值；的值； L 是下三角矩陣，用順代算法是下三角矩陣，用順代算法2求解方程組求解方程組 UX = Y 得向量得向量 X 的值。的值。U 是上三角矩陣，用回代算法是上三角矩陣，用回代算法2.3 Gauss列主元素消元法列主元素消元法 000. 3000. 2000. 1643. 5072. 1000. 2623. 4712. 3000. 1000. 3000. 2001. 0321xxx0.0012.0003.0001.0001.0003.7124.6232.0002.0001.0725.6433.000 20036006400101002300520040000. 1000

13、. 3000. 2001. 0A/b= 000. 2000. 5001002300520040000. 1000. 3000. 2001. 0 000. 1000. 3000. 2001. 0000. 2623. 4712. 3000. 1000. 3643. 5072. 1000. 2 002. 1003. 3001. 205000. 0801. 1176. 30000. 3643. 5072. 1000. 2A/b= 6870. 0868. 1005000. 0801. 1176. 30000. 3643. 5072. 1000. 2n解為解為x1=-0.4900，x2=-0.05113，

14、x3=0.3678nnnnnnnbaaabaaabaaa21222221111211在在A的第一列中選絕對值最大的元素作主元，設該元素的第一列中選絕對值最大的元素作主元，設該元素所在行為第所在行為第i1行，交換第一行與第行，交換第一行與第i1行，進展第一次消行，進展第一次消元；再在第元；再在第2n行的第二列中選絕對值最大的元素作行的第二列中選絕對值最大的元素作主元，設該元素所在行為第主元，設該元素所在行為第i2行，交換第二行與第行，交換第二行與第i2行，行，進展第二次消元，進展第二次消元，直到消元過程完成為止。直到消元過程完成為止。 6745150710623321xxx 6515707104

15、623 6515462370710 5 . 255 . 21 . 661 . 070710 2 . 62 . 65 . 255 . 270710列主元矩陣的三角分解：列主元矩陣的三角分解：2) 11032(EEE 3) 11053(EEE 21EE 5150710623A 51562307101APA 55 . 261 . 00710111APFAP 1000010101P 100100131211mmF 0101000012P 10010001322mF3)25 . 21 . 03(EEE 61 . 055 . 2071011211APFPAPF32EE 2 . 655 . 20710112

16、2APFPF2121PFPF 1211 FFL若若記記APFFU 12LUPA 可得可得12PPPPn)(1122111nnFPFPFPPL)(nAU LUPA )1()1(max kijnjknikkjiaakk)1( kjikka)1( kjikka)1( kkka2.4 矩陣三角分解法矩陣三角分解法一一 Doolittle分解法分解法二二 Crout分解法分解法三三 對稱正定矩陣的對稱正定矩陣的Cholesky分解法分解法四四 三對角方程組的數值解法三對角方程組的數值解法 nnnnnnaaaaaaaaa212222111211 1112121nnmmm nnnnuuuuuu2221121

17、1設設A非奇特，且非奇特，且ALU，L為單位下三角陣，為單位下三角陣，U為為上三角陣，即上三角陣，即一、一、Doolittle分解法：分解法：u11=a11, , u1n=a1nm21u11=a21, , mn1u11=an1對對A的元素的元素aij ，當，當 jk 和和 ik時時 kjkrrjkrkjuuma 11kkikkrrkirikumuma 11m21u12+ u22=a22, , m21u1n+ u2n=a2nm21=a21/u11, , mn1=an1/u11 u22=a22 - m21u12, , u2n=a2n- m21u1n m31u12+ m32u22=a32, , mn

18、1u12+ mn2u22=an2 m32=(a32- m31u12)/u22, , mn2=(an2- mn1u12)/u22 11krrjkrkjkjumaukkkrrkirikikuumam/ )(11 矩陣矩陣L和矩陣和矩陣U的元素計算公式：的元素計算公式：計算次序如下圖：計算次序如下圖：) 1 , 2, 1(/ )(/)5(), 3 , 2()4(), 1,;, 3 , 2(/ )()3(), 1,;, 3 , 2()2(), 3 , 2(/), 2 , 1() 1 (111111111111111nniuxuyxuyxniymbybynkkinkuumamnkkjnkumauniua

19、mnjauiinirririinnnnirririikkkrrkirikikkrrjkrkjkjiijjju1ja1 1111uamii 725101391444321131243301024321xxxx14131211uuuuTmmm413121130102T25 . 05 . 112423220uuu5 . 812110 11krrjkrkjkjumau343300uu11211300Tm43100T910044000u4000Tyyyy4321T1611/172010Txxxx4321nnnnuyx iinirririiuxuyx1T4321Tmm423210T11611310kkkr

20、rkirikikuumam1111irririiymby11by 432144434241343332312423222114131211bbbbaaaaaaaaaaaaaaaa 4535251544434241343332312423222114131211aaaaaaaaaaaaaaaaaaaa 4321444342413433323124232221141312111bbbyaaamaaamaaamuuuuk 4321444342413433323124232221141312112bbyyaammaammuuumuuuuk 4321444342413433323124232221141

21、312113byyyammmuummuuumuuuuk 4321444342413433323124232221141312114yyyyummmuummuuumuuuukULy 4535251544434241343332312423222114131211aaaaaaaaaaaaaaaaaaaaA 72510139144432113124330102 725101391444321131243301024321xxxx22100310317111220221334221162913711k 164911621117112113113212021712112310301024kULy 725

22、10139144432113124330102 71391422432215131242310301021k 713911621117112113113212021712112310301023k 44911623112113113212217121123130102yUx解解x 43214321xxxxx 164911621117112113113212021712112310301024kULy 111211221222111nnnnnnuuullllllULAnrriulalrkkrikirir, 1,11 nrrjlulaurkrrkjrkrjrj, 2, 1,/ )(11 例：例：

23、求矩陣求矩陣A的的Crout分解：分解： 303021112A)21() 1(222l)21(3032l)23/()21() 1(023u)31()23()21(3333l 1003/1102/12/1152/3302/31002A所以所以 52/333/12/312/12/12 nnnnnnnnuuuuuummmA, 122112111,121111故故A可進展可進展LU分解：分解： nnnnuuuuuuU222112111/1/11, 1, 111111122211nnnnnnnuuuuuuuuu1DUTTTLLLDLDLDLDUDLDLDULUA)(21212121121211 ),(2

24、211nnuuudiagD 22211),(nnuuudiag 那么那么 A=LDU1121211UDDDUU = AT = U1TDLT = U1=LT對對A的第的第i行元素行元素, 有有 ( j = 1,2,i ) ijjkjkikall 1jjjkjkikijijlllal/ )(11 21112)( ikikiiiilal( j =1,2,i-1 ) nnnnnnnnnnnnnnaaaaaaaaallllllllllll21222211121122212111212221111111al 對于對于 i = 2,3, n 1111112121212121nnnnnllldddlllA n

25、nnnnnnnnnnaaaaaaaaalllddldlddld2122221112112121221121211111jjkjkkikijijdldlal/ )(11 11ikikkikiiildladd1 = a11, l21= a21/d1, , ln1 = an1/d1( i = 2,3, n； j = 1,2,i-1 )yxLbLyTyDxLbLyT1 151531140150523121A 81621bd1 = a11, l21= a21/d1, , ln1 = an1/d1jjkjkkikijijdldlal/ )(11 11ikikkikiiildlad 132130121001

26、20001L 1000090000100001D 1150134324214144232131331212211ylylylbyylylbyylbyby 1/1/1/1/441331221111442332222443333444xlxlxldyxxlxldyxxldyxdyx四、四、 三對角方程組的數值解法三對角方程組的數值解法 nnnnnnnnnffffxxxxbacbacbacb121121111222110011|nniiiabcabcb 且 ai ci 0， (i = 2，3，n - 1) 其中其中三對角方程組方式如下：三對角方程組方式如下： 111112212211nnnnnbab

27、acb )1, 3 , 2(/), 3, 2(,/,111111nkcnkabacbkkkkkkkkk ), 2(,/ )(,/1111niyafyfyiiiii nnnyyyxxx212111111 )1 , 1(,1nixyxyxiiiinn nnnnfffyyy2121221 例：例： 求解三對角方程組求解三對角方程組 101653242231124321xxxx 5324223112A 2/536/55/1225/42/512/12 2/535/1222/512L 16/515/412/11U 1016Ly y=3, 8/5, 4/3, 2T 23/45/83Uxx = 5, 4, 3, 2Tnnnyyyy1111112211nnnnnnnfbafcbafcbafcb11112222111 nnnnnnnfbafcbafcbay11112222111 1 = c1/ b1 y1 = f1 / b1 nnnyyyy1111112211 例：例： 用追逐法求上例三對角方程組用追逐法求上例三對角方程組 101653242231124321xxxx 15302421231612其系數增廣矩陣為其系數增廣矩陣為 15302421231612 153024