數(shù)值分析課件解線性方程組的直接法

1、15第五章線性代數(shù)方程組的數(shù)值解法線性方程求解問題是科學(xué)研究和工程計算中 最常見的問題。如電學(xué)中的網(wǎng)絡(luò)問題、工程力 學(xué)中求解連續(xù)力學(xué)體（微分方程）問題的差分 方法、有限元法、邊界元法及函數(shù)的樣條插值、 最小二乘擬合等，都包含了解線性方程組問題。 因此，線性方程組的解法在數(shù)值計算中占有極 其重要的地位。Cramer）法則，）其中A為用向量對于n階線性方程組 Ax=b，若det（A）* 0,貝0 方程組有惟一解。由克萊姆（丿、）x de址）（i= 1,2| ndet（A）每個n階行列式nb代替A中第i列向量所得矩陣。共有n!項，每項都有n個因子，所以計算一個 階行列式需做（n-1）x n!次乘法，

2、我們共需要計 算n + 1個行列式，要計算出Xi，還要做n次除法， 因此用Cramer法則求解要做（n2 - 1廣n!+ n次乘 除法（不計加減法），計算量十分驚人。如n = 30 時，就需作約2.38X 1035次乘法。可見Cramer法 則在理論上是絕對正確的，但當(dāng)n較大時，在 實際計算中卻是不可行的。因此尋求有效的數(shù) 值計算方法就成為非常必要的課題。線性方程組的類型很多，若按其系數(shù)矩陣階 數(shù)的高低和含零元素多少，大致可分為兩類： 一類是低階稠密線性方程組，即系數(shù)矩陣階數(shù) 不高，含零元素很少。另一類是高階稀疏線性 方程組，即系數(shù)矩陣階數(shù)高，零元素占絕對優(yōu) 勢(比如占70%以上)。線性方程組

3、的數(shù)值解法也可分為兩大類：直 接法和迭代法。直接法是在沒有舍入誤差的情況下，通過有 限步運算可以得到方程組精確解的方法。但是， 在實際計算時，由于初始數(shù)據(jù)變?yōu)闄C器數(shù)而產(chǎn) 生的誤差以及計算過程中所產(chǎn)生的舍入誤差等 都要對解的精確度產(chǎn)生影響，因此直接法實際 上也只能算出方程真解的近似值。常用的有效 算法是Gauss消去法和矩陣的三角分解法。迭代法是用某種極限過程去逼近準(zhǔn)確解的方 法。如對任意給定的初始近似解向量 x(0)，按照 某種方法逐步生成近似解序列 x(0), x，川,x(k) ,111,使極限lim x=x*為方程組的解。因為在實際計 kT處算時，只能做到有限步，所以得到的也是近似 解。常

4、用的迭代法主要有Seidel迭代法、逐次超松弛法及共軛斜量法等。 直接法的優(yōu)點是運算次數(shù)固定，并且可以事 先估計Jacobi 迭代法、Gauss-。它的缺點是通常需要存儲系數(shù)矩陣和 常數(shù)項的所有元素，因而所需存貯單元較多， 編寫程序較復(fù)雜，一般適用于求解低階線性方 程組；迭代法的優(yōu)點是原始系數(shù)矩陣始終不變, 且只需存儲原方程組系數(shù)矩陣中的非零元素， 因而算法簡單，編寫程序較方便，且所需存貯 單元也較少，所以被廣泛用于求解高階稀疏線 性方程組。缺點是存在收斂性和收斂速度問題, 因而只能對具有某些性質(zhì)的方程組才適用。1消元法1.1三角方程組的解法 形如UiiXU12X2 + III+ UinXn

5、= yiI IIIIII(1.1)十_IUn_1,n -1Xn 一1Un 一 1,nX - Vn一 1丨Unn Xn 二 Vn的方程組稱為上三角形方程組。寫成矩陣形式U1, n 一1U1n1X11y11Ii11 -11 -i111 *丨=1 1Un-1,n -1Un T,n丨 IXn_11Wn111 1111unn nnL XnJL ynJU11III簡寫為Ux = y。若det(U p 0，則(1.1)有惟一解* = yn/Unnf n、/(k=n Tll,1)iXk = I yk -藝 UkjXj 丨/UkkLVj二k+1丿/我們稱求解上三角方程組(1.1)的過程為回代過程。同時也看到求解

6、這類方程組是件容易的 事，所以對一般形式的方程組，應(yīng)設(shè)法將它化 為(1.1)式的形式，然后再求解。1.2 Gauss消去法考慮方程組ax, + 2X2+ 111+ ainXn = III III 牯門必 + an2XP annXn其矩陣形式為blbn(1.2)Ax = bA =*w =A 1 * 1， x，b =總 III ann jMn丿1En丿IIIaiiXia1n種方法。 Gauss消去o化線性方程組(1.2)為等價的三角形方程組的 方法有多種，由此導(dǎo)出不同的直接方法，其中 Gauss消去法是最基本的 法分消元計算和回代求解兩個過程為后面的符號統(tǒng)一起見，記方程組(1.2)為A b(1)苴

7、中 Ab(1) = bo消元計算第一步：就是要將A的第一列主對角元以下 的元素全約化為0。設(shè)aj * 0，計算mi1 =玄丁/事(i = 2,3,|,n)用-mi乘(1.2)的第1個方程，加到第i個方程上， 完成第一步消元，得(1.2)的同解方程組a；, 11 X2 丨 I bP 丨:丨I ：I :丨axjb町a(chǎn)12)a22)IIIIII(1.3)I 0簡記為an22IIIAX = b其中A,b的元素的計算公式為 卜(2)i,j = 2,|,nibi(2) = bi-口曲，i = 2,111 ,n假設(shè)前k-1步消元完成后，得(1.2)的同解方程 組為G(1)1 a11IllIN事ib111a2

8、?illa2?IIIa2n 丨 1 X2 丨 1；II ： 11lb22) !1 : 1111IIIakn1 Xk;II : I 1!畀!1 : 111Vank)III1 * 1 aiJ八Xn丿(1.4)簡記為 A(k) X 二 b(k)。第k步：就是要將 A(k)的第k列主對角元以下的元素全約化為0。設(shè)akk 0，計算 mk = ai(kk)/akk)(i = k + 1,111, n) 用- mk乘(1.4)的第k個方程加到第i個 III, n)方程上，完成第k步消元。得同解方程組A 5 x=b(其中A(k), b)元素的計算公式為何-y/mkakk：山、b(k+1) b(k)b(k)(i

9、,J = k+1,|,n)b =b -mkd繼續(xù)上述過程，完成n - 1步消元后，(1.2)化成 同解的上三角方程組 闔aj!a22)(k+1)k+1)(i 二 k+ 1,illIIIain)V；-1n II , a22n !x2i:II ： I aiJ ； X nJ號(1.5)簡記為 A(n)x = b(n)。E解因 det(A) * 0，故 ann)|xn = bnn)/ann)rniXk = I bik)-送aklk)Xl IILV !15左乘矩陣A,b，即LiAbA,b一般地，第k次消元等價于用初等矩陣-g,k左乘矩陣A(k), b(k)，即A(k1), b(k1)b L1 A(k),

10、 b(k) 經(jīng)過n-1次消元后得到Ln_1LllL1A b = A(n),b(n)n1円丄/吐小二A(n) 丄1Ln_2ll|L1b(1)= b(n) 將上三角矩陣A(n)記為U，得到A二 L11L21|Ln1_1U = LUL 二 LIIlL爲(wèi)了 1II叫11m32mn2 Hi為單位下三角陣。這說明，消元過程實際上是把系數(shù)矩陣A分解為單位下三角陣與上三角矩陣的乘積的過 程。這種分解稱為杜利特爾(Doolittle )分解， 也稱為LU分解。定理2 (矩陣的LU分解)設(shè)A為n階方陣，如果A的順序主子式Dj* 0 (i = 1,111, n-1),則存 在惟一的單位下三角陣L和上三角陣U，使A

11、= LU。o性。如果A非奇異，設(shè)A有證明 以上的分析已證明了 A可作LU分解 下面證明分解的惟 兩個分解式A= LU =屮！L,Li都是單位下三角陣，U ,Ui都是上三角 陣。因A非奇異，所以L,Li,U,Ui都可逆。于是 1；工=UQT 左邊為單位下三角陣，右邊為上三角陣，所以L2 uum I01u u2il,U =丨I1L 0UnnL1為n T階單位下陽L1I,L= IannL L31階非奇異陣,LU，得即有L, = L,Ur U，惟一性得證。 若A為奇異陣，因其n-1階順序主子式不等 于零，故可記A AlA= ILA3苴中A為n- 角陣，由AL3U 1，ann L3U2 + unnA =

12、 LUi,Ar 屮2,Ar由已證明的結(jié)果知 A二LU1的分解是惟一。所 以L；1A2,L3= A3U；1亦是惟一確定的。進(jìn)而，unn也是惟一確定的定理2的逆命題是：定理3設(shè)n階方陣A非奇異，若A有惟一的分 解A=LU（ L為單位下三角陣，U為上三角陣）， 試證A的順序主子式Di主0,i = 1,2,111,n-1。（證明留給讀者）2主元素法前面已指出Gauss消去法必須在各個約化主 元素ak：） * 0 （k =n-1）下才能進(jìn)行?，F(xiàn)在還需指出的是：即使akk 0，但當(dāng)|akk很小時，也 不能做主元素的，因為作第k次消元時，需將第 k 個方程乘以-miik = - ai（kk）/ a：）（i=

13、k+1）ll,n）,因 此當(dāng)akk）很小時，乘數(shù)-gk很大，用-mk去乘第 k行的元素akj （j = K +1,111,n），將導(dǎo)致Fkakj的數(shù) 量級迅速增長，這樣在消元計算aij - miik akj(i = k+1,|, n)時，會出現(xiàn)大數(shù)吃掉小數(shù)的現(xiàn)象， 因而導(dǎo)致最后計算結(jié)果的精度很差甚至失真。-2例1解方程組+ X2 = 1L X1 + X2 = 2(1)求精確解；(2)在 P法求解。解(1),1為=1 +99容易求得方程組的精心1,。確解為=10, n = 2的浮點機上用 Gauss消去(2)在所給浮點機上原方程組為I0.10X 10Tx1 + 0.10X 101x 0.10X

14、101i0.10x 101X 0.10X 10,2 = 0.20X 101由于m21 = 0.10 咒 101/(0.10 咒 10一1) = 0.10 103 消去第二式的X,，得(0.10咒 101 - 0.10咒 103)x2 二 0.20咒 101 - 0.10咒 103 對階，出現(xiàn)大數(shù)“吃掉”小數(shù)，結(jié)果有0.10咒 103X廠 0.10汽 103代第一式得X廠0。解得X廠0.10T01，m2i很大。改變 二兩式對換，與精確解相比較，已無精確度可言。產(chǎn)生不 準(zhǔn)確的原因是主元素太小，致使 上述狀況的辦法是將方程組的一、 得101x0.10 咒 101x210-1x1+0.10 滅 101

15、x2Gauss消去法，此時0R 101二 0.10咒 101I 0.1010.10然后再使用m21 = 0.10咒 17(0.10 10)= 0.1 10 (0.10X 10- 0.10X 10 1)x 0.10X 10- 0.20 10 1 于是得近似解 為=0.10X 101,x 0.10X 101。此例表明，高斯消去法解方程組時，小主元 可能帶來嚴(yán)重的后果，因此應(yīng)盡量避免小主元 的出現(xiàn)；另一方面，通過對換方程組中方程的 次序或改動變量次序，選擇絕對值大的元素做 主兀，可減少舍入誤差，提咼計算精度。2.1列主元與全主元消去法k),然后再進(jìn)1.列主元素消去法。在第k步消元時，在A 的第k列元

16、素ai(kk) (i 2 k)中選取絕對值最大者作 為主元，并將其對換到(k,k)位置行消元計算。用方程組(1.2)的增廣矩陣a111a12IIIa1nbl、11_ _ - _ 1 a21a22IIIa2nb2 1A,b =!：*I 4 4*:14-1lan1an2IIIann1 bn丿(2.1)表示它，并直接在增廣矩陣上進(jìn)行計算。具體 步驟如下：第一步：首先在上述矩陣的第一列中選取絕 對值最大的，比如，則 aiiimax ai11n。將(2.1)中第一行與第h行互換。為方便起見，記行互換 后的增廣矩陣為A,b，然后進(jìn)行第一次消 元，得矩陣ajIIIaTba2l)IIIa22n)b22)!HI

17、III11IIIan?bn2)丿V 0a(2) ,b (2)+ 0假設(shè)已完成k-1步的主元素消去法，約化為a(1)a1kfe1IIIa(1)a1nIfebTi.1akkk)曲III*a(k) dkn+f.*ank)III*a(k) dnn事III如a；)，使aikk二 maxkEi n第k步：在矩陣A(k) ,b(k)的第k列中選主元， aikk)|。將A(k) ,b(k)的第 k行 與第ik行互換，進(jìn)行第k次消元。經(jīng)過n-1步， 增廣矩陣(1.7)被化成上三角形bIb22)iaiT 事 III ai(n)! a)? Ill a22n)ann)bnn)j算法1 (Gauss列主元素消去法)說明

18、：消元結(jié)果沖掉 A，行乘數(shù)mk沖掉aik， det存放行列式值。輸入 A = (aij),b= (bi,b2,lll,bn)T,n，置 kT，det = 11.消元過程對 k = 1,|, n-1(a)選主元：按列選主元Pk，即確定arkr，使(b)(c)akjbk -Pk三若Pk = 0,停機； 若r * k (進(jìn)行行交換)-arj (p k,k + 1,|,n)br二 max aki k)仍是一個指標(biāo)為k的單位 下三角陣。若令P = pn_j,n-1 pnd,n-2 III pi,1則p是排列陣。于是歸iClIlLi PA= U注意 陣，若令PA=(Lyiii(LL2)j(LnJu到(Li

19、rill Ln2)LiJ仍為單位下則有PA= LU上式表明，帶行交換的 Gauss消去過程所產(chǎn)生 的矩陣分解，相當(dāng)于對系數(shù)矩陣先施行每步消 元時所做行交換后，再將所得矩陣PA進(jìn)行LU分解。以上討論可敘述為定理4 (列主元三角分解定理)如果A為非奇異矩陣，則存在排列矩陣 P，使PA= LU35L為單位下三角矩陣，U為上三角矩陣3矩陣三角分解法3.1直接三角分解法1.不選主元的直接三角分解法 設(shè)A = LU為如下形式戸耳2rI a21 a22為 I = I b 1ainU12 小 U1n JU22 J (3.1)V u1l| 1IIH1AIII Oin 丿 Um In2 川由矩陣的乘法規(guī)則，得a1

20、j = U1j (j = 1,2,111, n)Van1an2n1 n2Unn丿mi n(i,j)aij =無likUkjk=1r ijZ hkUkj + uijk：1(i=2,|,n)I.I 藝 1 ikUkj +IijUjj j V iL k=1由此可得計算lij和Uij的公式I Uj = aljIi-1! Uij = aij -送 likukjI心I1j-1 I廠佝-送 likUkj )iujj2具體步驟如下：1)計算u的第1行, gj = a1j ( j = 1,2, |,n);2 )計算U的第r行,(j = 2,111, n)(i= 2,111,n, j = UII,n) (3.2)

21、(j = 1,111,n,i = j + 1川 n)lir 二丄(airUrr例2求矩陣r -1Urj = arj -送 IrkUkjk=1r -1-藝 IikUkr ) ( i 二k=1L的第1列l(wèi)i廠 aii/uii (i 二 2,III, n)L的第 r列(r = 2,lll,n)(j = r,r + 1,lll, n)r +1,111, n,r H n)2i-2317;5j的三角分解。解按式(3.2)Uiil21U22U23l32U33所以印1 二 2,U12 二 312 二 2,U13 二 313 二 3a21/u11 = 4/2 = 2,131 = a31/u11 = -2/2=

22、-1322 l21u12 = 7 - 2咒 2 = 3323 一 I21U13 = 7 - 2 X 3 = 1(a32 l31U12VU2 - 1 2/3 = 2a33 (l31U1 l32U23(尸 3 十 2咒 “ =6012310| |0 3 1 I1H06j(311) U11(a12)U12(a21)l21(a22)U22(a31) l31fe(a32) l32fe-(an1)ln1(a n2)ln2(an3) ln3(a13)U13(a1n) U1n23)U23 (a2 n) U2n33)U33 (a3 n) U3n(a nn) Unn例2中矩陣A的三角分解按緊湊格式計算， 結(jié)果見下

23、表(2) 23I (2) 24,(4) 7=22(-2)手=-12(7) 7 - 2咒 2 = 3(7) 7-2咒 3=1 144 -(_1)c 2 = 2 (5) 5- 2% 1 (1)% 3= 6 3如果線性方程組Ax=b的系數(shù)矩陣已進(jìn)行三 角分解，A= LU，則解方程組Ax= b等價于求 解兩個三角形方程組Ly = b,Ux = y。(1Ly121UniIIIVII1.1丨III1八yn丿Ibn丿第一步：先求解下三角方程組Ly=b#yr得解4kTk= 2,3,|,n丨yk二bk-無和Ij=1第二步：再求解上三角方程組 Ux= yF U11U12IIU22IIIIIIUin 1f yI.

24、I I.I IU2n I X2 I y2 IUnn丿I Xn丿 In丿47得解*|Xkyn/Unn 丄Wk ukknZj=k+1UkjXj), k 二 n- 1,n- 2,|,12 .選主元的直角三角分解法不選主元的直接三角分解過程能進(jìn)行到底的 條件是 奇異，也可能出現(xiàn)某個主 0(r = 1,111,n-1),實際Urr = 0的情況，這時分解 主0但很 小，會使計算過程中的舍入誤差急劇增大，導(dǎo)Urr過程將無法進(jìn)行下去。另外，如果,即使A非Urr致解的精度很差。但如果 A非奇異，我們可通 過交換A的行實現(xiàn)矩陣PA的LU分解，實際上是 采用與列主元消去法等價的選主元的三角分解 法，即只要在直接三

25、角分解法的每一步引進(jìn)選主元的技術(shù)即可。設(shè)r -1步分解已經(jīng)完成，這時有IIIUiiI l21IIIU22l32IllIIIIIIUinAtUr-1,r-1I r ,r -1arrUr1,n Iarn IIII考慮r步分解，對于r =1,2,III, n做到 (1)計算輔助量S , S沖掉airrjll, n)In1ln2III1 n,r -1anrannr-1air 1 S = air -送 likUkr 0k=1(2) 選主元。即確定行號srl = maxim,并記錄主行ir，使口.,號.1 p(r) ir ( Ip (n)為一維數(shù)組)(3) 交換a的r行與ir行元素 a廠 air,j (p

26、 1,111,n)(4) 計算U的r行與L的第r列元素 arr Ur廠 ar廠 Srarjr 1u廠 arj -三 IrkUkj (j 二 r+ 1,111,n)k =1hr 二 s/urr = dir /(H r + 1,|11, n,r = n) 蘭1)air *（這里|lir上述計算過程完成后，就實現(xiàn)了PA的LU分解，U保存在A的上三角部分，L保存在A的嚴(yán) 格下三角部分，排列陣 得。例如，設(shè)n= 4,1 PP = 1 341 231 12P由Ip (n)的最后記錄可=2, I0!0I0IL1PR)100001003,lp(3)= 4,貝 y010!11I0于是a21a22a23a24 1

27、11 a31a32a331a34 1$41a42a43a44 jLai1a12a13耳4PA =方程組Ax= b的求解等價于解 PAx = Pb，即LUx = Pb可通過Ly = Pb及Ux = y完成。例3用選主元的直接三角分解法解方程組其14-1-4-9316, b= (1,4,1)t,x= (X1,X2,X3)T2間結(jié)果沖 掉A相應(yīng)位置元素，數(shù)組I p（3）記錄主行號 第一步分解：因為S3對A作選主元的LU分解，中a31 1S3 = 4，八13A4；21-9-4-1第二步分解:Sb = max32肖今S，故有4IApq-J 丨114!1L2-9-1-4口亠。印廠 S f 1 , a=ma

28、x S10 QCi 工 0,i 二2,111, n-1 (3.4)則它可分解為1A = LU = I 21 ru11 III IIn 1L! (3.5)Cn_1 IIUnCT 0 (i = 2,|,n-1)為矩陣(3.3)中所給出，且分解是惟一的證明將式(3.5)右端按乘法法則展開，并與A的元素進(jìn)行比較，得；b U1Vai = hUiT(i = 2,111, nT)g = Ci + Ui如果 Ui H 0 (H 2川，n- (i = 2,111, n)，即5 = b1Hi = ai/Ui_1(i = 2,|, n) (3.6)I Ui = h - ci-1ii從以上的公式和消元過程我們可看出，

29、要使 得分解得以實施，必須滿足Ui 主 0 (i = 2,111, n1)現(xiàn)在我們用歸納法證明：q當(dāng)i = 1時，顯然5 = b, 成立，即 有UiUiUkUk-1Ck-11，那么我們可計算li,Ui 0(i = 1,2,111,n- 1)。cj成立。假定i = k- 1 ，我們將證明H k時也成立bkUk-1uk1akUk-1bk -Ck二 bk - c-ilkbkUk-1 -q -1 akUk-1CkiakUk-1當(dāng)矩陣(33)按式(3.6)計算進(jìn)行分解后，求解方程組可化為求解方程組 Ly = d及Ux = y。解Ly 二 d 得5)53(y廠 di o,m= 1,HI,n，由此可推出 d

30、m 0,mT,|,n。記1 D2 = diag(Vd?,|,7d7)則有11 1 1A L1D2D2lT = (L1D2)(L1D2) llt1其中L= LD：，它為對角元為正的下三角陣，。由分解l1dlT的惟一性，可得分 性。證畢。所以(4.2)成立解(4.2)的惟一，分解式A = L L稱為正定矩陣的喬列斯基(Cholesky)分解。利用Cholesky分解來求系 數(shù)矩陣為對稱正定矩陣的方程組AX = b的方法稱為平方根法。用比較法可以導(dǎo)出L的計算 公式。設(shè)ii LI21kin11 nn丿Un1I1 n1III1 n,n1比較A與LLt的對應(yīng)元素，可得j-i1(ajj Y ijk)2 (j

31、k=11jj-Gj-w likijk) (i = ijjkT1,2,111, n)(4.3)廠1川，n)0這里規(guī)定送-0。計算順序是按列進(jìn)行的k nI11T Ii1(i = 2,3,111,n)T I22T Ii2(i=3,|,n)oT III當(dāng)矩陣A完成Cholesky分解后，求解方程組Ax= b就轉(zhuǎn)化為依次求解方程組Ly = b, L X = y它們的解分別為fi-1yi = (bi-2 likyQ/lii (i = 1,2,III,n)k；/(4.4)! Xi = Wi - 2 ikiXk) /lii (i = n,n - 1ll,1)Ij =i +1Lt的元素也就求出，所 n3/6次乘

32、除法，大約為一般LU 半。同時還節(jié)省存貯，只當(dāng)L的元素求出后， 以平方根法約需 分解法計算量的 存系數(shù)矩陣A的下三角部分和右端項b，中間結(jié) 果L沖掉A，計算解沖掉b。此外，由式(4.3)的 第一式得aj廠送 l2k, i = 1,2,111,n k=1所以21 jk ajj maxx(ajj)上式表明，在矩陣A的喬列斯基分解過程中 ljk|的平方不會超過A的最大對角元，舍入誤差 的放大受到了控制，且對角元素Ijj恒為正數(shù)，于是不選主元素的平方根法是數(shù)值穩(wěn)定的，計 算實踐也表明了不選主元已有足夠的精度，所以平方根法是目前計算機上解決這類問題的最 有效的方法之一。4.2改進(jìn)的平方根法（LDLt法）利用平方根法解對稱正定線性方程組時，計 算矩陣L的元素Ij時需要用到開方運算。另外, 當(dāng)我們解決工程問題時，有時得到的是一個系 數(shù)矩陣為對稱但不一定是正定的線性方程組， 為了避免開方運算和求解對稱（未必正定）方 程組，我們可以引入下面 改進(jìn)平方根法。L為單位下三角陣，D為由定理5我們知道，A為對稱矩陣時，它可 分解成A = LDLt，其 對角陣。記（1l，D = diag(di,d2,lll,dn)I1丿Ijj = hijk = 0( j k)，得1 l