數(shù)值計(jì)算方法線代方程組直接法

1、1 n階線性代數(shù)方程組的一般形式稱為右端常數(shù)向量。稱為解向量，稱為系數(shù)矩陣，其中bxbbbbxxxxaaaaaaaaaAbxaxaxabxaxaxabxaxaxannnnnnnnnnnnnnnnnnA.212121222211121122112222212111212111解線性代數(shù)方程組的數(shù)值解法有：解線性代數(shù)方程組的數(shù)值解法有：直接直接法，法，迭代法。迭代法。在沒(méi)有在沒(méi)有舍入誤差舍入誤差的情況下，經(jīng)過(guò)有限次的情況下，經(jīng)過(guò)有限次運(yùn)算可以得到方程組的運(yùn)算可以得到方程組的精確解精確解的方法。的方法。線性代數(shù)方程組的直接解法線性代數(shù)方程組的直接解法/*Direct Method for Solvi

2、ng Linear Algebraic Systems*/求解求解,n nAxb AR 0det( )A Cramer法則法則: :1 2, ,iiDxinD所需乘除法的運(yùn)算量大約為所需乘除法的運(yùn)算量大約為( (n+1)!+)!+nn=20時(shí)，每秒時(shí)，每秒1億億次運(yùn)算速度的計(jì)算機(jī)要算次運(yùn)算速度的計(jì)算機(jī)要算30多萬(wàn)年！多萬(wàn)年！直接法直接法 Gauss消去與矩陣LU分解一 Gauss消去1 直接法的關(guān)鍵思想如果方程組是“上三角方程組”或“下三角方程組”就可以很容易地求出方程組的解。因此，直接法的關(guān)鍵思想就是如何把一般方程組約化為上/下三角方程。屬于解方程的直接法屬于解方程的直接法2 上三角方程組與

3、回代過(guò)程nijiijijiinnnnnnnnnnnnnnniaxabxabxbxabxaxabxaxaxaRAbAx12222211212111) 1 , 2,.,1(/ )(/組的解代過(guò)程”可以求的方程則通過(guò)自下而上的“回的上三角方程組非奇異）被約化為如下（如果方程組3 下三角方程組與前推過(guò)程1 - i1i1111221122221121111).,32(/ )(/jiijijinnnnnnnnniaxabxabxbxaxaxabxaxabxaRAbAx，組的解推過(guò)程”可以求的方程則通過(guò)自上而下的“前的下三角方程組非奇異）被約化為如下（如果方程組4 Gauss消去過(guò)程)1(1nn,)1(nn

4、)1(n2)1(n1)1(1n2,)1(2n)1(22)1(21)1(1n1,)1(1n)1(12)1(11)1()1(1,1,21,121aaaaaaaaaaaab|A表示成增廣矩陣形式將記為為了表達(dá)方便，將bAxaaabbbbbnnnnn則Gauss消去過(guò)程如下：) 1,.,1;,.,1(),.,1(/) 1,.,3 , 2;,.,3 , 2(), 3 , 2(/1)()()1()()()1(11)1()2()1(11)1(11nnkjnkialaankiaalknnjnialaaniaalkkjikkijkijkkkkikikjiijijii步：第步：第該該Gauss消去法為順序高斯消去

5、法消去法為順序高斯消去法ijijikkjaaa a ;ikikkkaaa 1 21, ,knfor11,jknfor12,ikknforGauss法的法的消元消元過(guò)過(guò)程程回代回代過(guò)程過(guò)程11,;nnnnnnaaa 1111,()ni ni jj nj ii ni iaa aaa 11,in 65546237710Gauss32132121xxxxxxxx消去法求解方程組例：用順序高斯消去法必須保證第順序高斯消去法必須保證第k 步的步的akk0，因?yàn)樗谙ミ^(guò)程和回代過(guò)程中起關(guān)鍵作因?yàn)樗谙ミ^(guò)程和回代過(guò)程中起關(guān)鍵作用，所以稱它為主元素。用，所以稱它為主元素。例：例：用用基本基本Gauss消元法

6、消元法求解下列方程組求解下列方程組123123123223347712457xxxxxxxxx 解：解：增廣矩陣增廣矩陣223347712457()A b 223303150684 223347712457()A b 223303150684 223347712457121121112( )( )ala 131131111( )( )ala 232232222( )( )ala 223323151684 223323151266 321122,xxx 基本基本Gauss消元法的工作量消元法的工作量消元過(guò)程：消元過(guò)程：11111()()()nnkknknk nk 回代過(guò)程：回代過(guò)程：12()n

7、n 325326nnn332333nnnn加減法的次數(shù)加減法的次數(shù)32353263nnnn乘除法的次數(shù)乘除法的次數(shù)3()O n基本基本Gauss消元法的實(shí)現(xiàn)條件消元法的實(shí)現(xiàn)條件全不為全不為零零的充要條件是的充要條件是1 2( )(, , ()iiiaik kn的的順序主子式順序主子式都不等于都不等于零零，即，即A11121212221201 2, , ()iiiiiiiaaaaaaiknaaa 小主元小主元 可能可能導(dǎo)致計(jì)算失導(dǎo)致計(jì)算失敗敗例：例：在在8位制計(jì)算機(jī)上解方程組位制計(jì)算機(jī)上解方程組912121012xxxx 要求用要求用Gauss消去法消去法計(jì)算。計(jì)算。921211110/laa

8、9992221110 001 101010. .al 8個(gè)個(gè)92212110bl 999101101010 2110,xx 解：解：121xx主元素要作除數(shù)，其絕對(duì)值相對(duì)于被除數(shù)主元素要作除數(shù)，其絕對(duì)值相對(duì)于被除數(shù)過(guò)小會(huì)影響到計(jì)算的精度，所以，我們可過(guò)小會(huì)影響到計(jì)算的精度，所以，我們可以采用以采用5、列主元、列主元Gauss消去法解方程組消去法解方程組除了除了“列主元消去法列主元消去法”，還有，還有“全主元消全主元消去法去法”，但后者計(jì)算量過(guò)大，一般不用。，但后者計(jì)算量過(guò)大，一般不用。思思想想每次每次消元之前消元之前，在剩余元素中選擇，在剩余元素中選擇絕對(duì)值最大絕對(duì)值最大的的非零元素作為主元，

9、非零元素作為主元， 然后經(jīng)過(guò)然后經(jīng)過(guò)換行換行換到換到主元位置主元位置 列主元列主元消去法消去法/* Column Pivoting Strategies */Step k：第第k步首先選擇步首先選擇主元主元1 21, ,kn尋求尋求 滿足滿足ki1( )( ),max,kkkiki kk i naaik kn 然后交換矩陣然后交換矩陣 的第的第 行和行和 行，再進(jìn)行消元過(guò)程行，再進(jìn)行消元過(guò)程 ( )kAkki1 21( )( )( )( ),;, ,kkkkkkkjkjijijtaaaat jn 算法算法: : Gauss列主元列主元消去算法消去算法求方程組求方程組Ax=b 的解的解. .輸入

10、輸入: :增廣矩陣增廣矩陣An(n+1)=（A|b）.輸出輸出: 近似解近似解 xk=ak,n+1(k=1,2,n) 或失敗信息或失敗信息.消元消元過(guò)程過(guò)程 for k = 1,2,n-1 do Step 1 - Step 4 Step 1 尋找行號(hào)尋找行號(hào) ik , 使得使得Step 2 如果如果 ，則交換第，則交換第k行和行和ik行；行； 否則轉(zhuǎn)否則轉(zhuǎn)Step 7 1,max,kiki kk i naaik kn 0,kika 算法算法: : Gauss列主元列主元消去算法（續(xù)）消去算法（續(xù)）Step 3 for i=k+1,n 計(jì)算計(jì)算 Step 4 for j=k+1,n+1 計(jì)算計(jì)算

11、回代回代過(guò)程過(guò)程Step 5Step 6 for i=n-1,1 計(jì)算計(jì)算 Step 7 Output (系數(shù)矩陣奇異系數(shù)矩陣奇異); /*不成功不成功 */ STOP.ikikkkala ijijikkjaal a 1,n nnn naxa 11,()/nii ni jji ij ixaa xa 20111 . 031045321321xxxGauss消去法解方程組例：用列主元例：例：用用Gauss列主元消列主元消去法求解下列方程組去法求解下列方程組123123123223347712457xxxxxxxxx 解：解：首先寫(xiě)出首先寫(xiě)出增廣矩陣增廣矩陣223347712457()A b 477

12、1223324 57 12i Step 112rr47711233214 572 21211112ala 31311112ala 消消 元元47713511222215171312222 23i Step 223rr47711517131222235112222 32312215ala 15 消消 元元47711517131222266112555 47715171222611255 12 2 全主元全主元消去法消去法/* Complete Pivoting Strategies */Step k：第第k步選擇步選擇主元主元1 21, ,kn( )( ),maxkkkkiji jk i j n

13、aa 尋求尋求 和和 滿足滿足kikj然后交換矩陣然后交換矩陣 的第的第 行和行和 行，第行，第 列和列和 列列( )kAkkikkj記錄記錄交換信息交換信息二 LU分解1 矩陣的三角分解對(duì)方程組Ax=b的順序Gauss消去過(guò)程的結(jié)果，就是把矩陣A分解成兩個(gè)三角距陣L和U的乘積：A=LU。利用這個(gè)特點(diǎn)可以進(jìn)行線性方程組的直接三角分解法。解方程組的直接三角分解法有3種形式：（1）A=LU（2）A=LU（3）A=LDU單位下三角陣單位下三角陣上三角陣上三角陣下三角陣下三角陣單位上三角陣單位上三角陣單位下三角陣單位下三角陣對(duì)角陣對(duì)角陣單位上三角陣單位上三角陣A的的Doolittle分解分解A的的Cr

14、out分解分解A的的LDU分解分解2 解方程組的直接三角分解法Ax=bLUx=byUxbLyDoolittle分解法分解法本本質(zhì)質(zhì)它是基本它是基本Gauss消元法的一種等價(jià)變形消元法的一種等價(jià)變形 Gauss消元法的矩陣形式消元法的矩陣形式 /* Matrix Form */11111221nnLLL L AU 為為上三角上三角矩陣矩陣U其中其中1,kkl 111, nklkL 1,kkl 111,n kl 1kL 1221nnAL LLLU LU 為單位為單位下三角下三角LL 32l1112nl121l31l1nl141l42l43l3nl1nnl nkl矩陣分解為矩陣分解為單位下三角單位下

15、三角和和上三角上三角矩陣的乘積矩陣的乘積ALU 11121222.nnnnuuuuuuU 710413221232321xxx解方程組例：用直接三角分解法1 21, ,knfor12,jkknfor12,ikknfor111()kikikkrrjrkkaaa aa 1111,kkjkjkrrjraaaa 計(jì)算計(jì)算 的的k列列L計(jì)算計(jì)算 的的k+1行行U矩陣矩陣 分解的實(shí)際計(jì)算公式分解的實(shí)際計(jì)算公式：LUfunction L,U,flag=LU_decom(A)n,m=size(A);if n=m error(The rows and columns of matrix A must be eq

16、ual!); return;endL=eye(n);U=zeros(n);flag=OK;for k=1:n for j=k:n z=0; for q=1:k-1 z=z+L(k,q)*U(q,j); end U(k,j)=A(k,j)-z; end if abs(U(k,k)eps flag=failure; return; end for i=k+1:n z=0; for q=1:k-1 z=z+L(i,q)*U(q,k); end L(i,k)=(A(i,k)-z)/U(k,k); endend課堂練習(xí)134627112413421321321321xxxxxxxxxLUGauss分解法

17、求解線性方程組消去法和、用答案312361225510134236122551013421688255101342131462711241342332321xxxxxxx再由回代過(guò)程得到組于是，得到上三角方程示為解：消去過(guò)程用矩陣表 31236251125103421543121111413221232332322131211323121xbUxybLyAuuuuuulllLUA解得由解得由可得，即令30161017959535332,12341321xxxCholeskyGaussLUA。分解法解下面的方程組消去法和、試用分解。作矩陣的、設(shè)有矩陣習(xí)題： Cholesky分解唯一的。元素為正時(shí)

18、，此分解是的對(duì)角且當(dāng)限定分解）（該分解為使得，實(shí)的下三角陣對(duì)稱正定，則存在一個(gè)定理：設(shè)LCholeskyLLALRATnn,該方法也稱為平方根法。該方法也稱為平方根法。思思想想Axb TLL xbTLybL xy Cholesky分解的計(jì)算公式分解的計(jì)算公式()ijn nAa 設(shè)設(shè)由由 對(duì)應(yīng)元素相等得對(duì)應(yīng)元素相等得TALL L 32l2nl21l31l1nlnkl11l22l33lnnl1jijirjrrijal l 11jijijjjirjrral ll l Cholesky分解公式分解公式11()jijirjrrijjjal lll 1 21, ,jnijn 112211 2() , ,j

19、jjjjjrrlaljn 因因?qū)ΨQ性對(duì)稱性無(wú)需存儲(chǔ)無(wú)需存儲(chǔ)11l21l31l1nlStep122l32l2nlStep233l3nlStep312a13a1na23a2na3nannlStepn的計(jì)算過(guò)程的計(jì)算過(guò)程: :L逐逐 列列 計(jì)計(jì) 算算1 2(;, , )iklik kn元素元素A仍然仍然存放在矩陣存放在矩陣 的相應(yīng)位置上的相應(yīng)位置上function L,flag=chol_factor(A)n=length(A);L=zeros(n);flag=OK;for k=1:n delta=A(k,k); for j=1:k-1 delta=delta-L(k,j)2; end if del

20、ta0，其中|x|=0當(dāng)且僅當(dāng)x=0,稱為正定性； |kx|=|k|x|,kR,稱為齊次性； |x+y|x|+|y|,且x,yRn,稱為三角不等式。), 1 ,)|(|2|maxxxpxxxxxxpnipipniiininii范數(shù)：一般地，定義范數(shù)：范數(shù)：范數(shù)：種向量范數(shù)、常用的例：設(shè)x=(2,-4,3)T,求|x|1,|x|2和|x|的一個(gè)矩陣范數(shù)。上為則稱，）（；，）（；，其中）（滿足：并記為）對(duì)應(yīng)的一個(gè)實(shí)值函數(shù)（若與設(shè)矩陣、定義二、矩陣范數(shù)AR| |A|,|,|)4(;RBA| |B| |A| |BA|3R| |A| |A|20A0| |A|0| |A|1| |

21、A|A,RA1nnnnnnnnRBABAAB)max()(|2|max|max|32211111AAAAAaAaATTniijnjnjijni范數(shù)：列范數(shù)：行范數(shù)：種矩陣范數(shù)、常用的稱為矩陣稱為矩陣ATA的譜半徑的譜半徑例：例：給定矩陣給定矩陣210111012A 求矩陣求矩陣 的的1、2、 范數(shù)。范數(shù)。A 13A 3A 210111012TAA23A 若若 是是實(shí)對(duì)稱實(shí)對(duì)稱矩陣，則矩陣，則A2( )AA 矩陣矩陣 的特征值為的特征值為A0 2 3, , 46. 522115|2 A0430201414102AAIT。，求例：設(shè)21| ,| ,|4321AAAA22115解得7|4|3| |,2| 1max|A解：2014141043