線性方程組的四種數(shù)值解法

1、線性方程組的四種數(shù)值解法 (電子科技大學(xué)物理電子學(xué)院，四川 成都 610054)摘要：本文介紹了四種求解線性方程組的數(shù)值解法: 雅克比迭代法、高斯賽德爾迭代法、高斯消去法和改進(jìn)的平方根法的基本原理和算法流程，通過求解具體方程，對(duì)四種求解方法進(jìn)行了對(duì)比。對(duì)于雅克比迭代法和高斯賽德爾迭代法，研究了兩種算法對(duì)求解同一方程組的迭代效率差異，結(jié)果表明高斯賽德爾迭代法達(dá)到同樣精度所需迭代次數(shù)較少。對(duì)于高斯消去法，通過選擇列主元的方法提高算法的準(zhǔn)確度，計(jì)算結(jié)果表明高斯消去法計(jì)算精確，且運(yùn)算復(fù)雜度也不是很高。對(duì)于改進(jìn)的平方根法，其運(yùn)算復(fù)雜度低，但對(duì)于給定的方程組有著嚴(yán)苛的要求。關(guān)鍵詞：雅克比迭代法；高斯賽德爾

2、迭代法；高斯消去法；改進(jìn)的平方根法；線性方程組引言線性方程組的求解在日常生活和科研中有著極其重要的應(yīng)用，但在實(shí)際運(yùn)算中，當(dāng)矩陣的維數(shù)較高時(shí)，用初等方法求解的計(jì)算復(fù)雜度隨維數(shù)的增長非常快，因此，用數(shù)值方法求解線性方程組的重要性便顯現(xiàn)出來。經(jīng)典的求解線性方程組的方法一般分為兩類：直接法和迭代法。前者例如高斯消去法,改進(jìn)的平方根法等，后者的例子包括雅克比迭代法，高斯賽德爾迭代法等。這些方法的計(jì)算復(fù)雜度在可以接受的范圍內(nèi)，因此被廣泛采用。一般來說，直接法對(duì)于階數(shù)比較低的方程組比較有效；而后者對(duì)于比較大的方程組更有效。在實(shí)際計(jì)算中，幾十萬甚至幾百萬個(gè)未知數(shù)的方程組并不少見。在這些情況下，迭代法有無可比擬

3、的優(yōu)勢(shì)。另外，使用迭代法可以根據(jù)不同的精度要求選擇終止時(shí)間，因此比較靈活。在問題特別大的時(shí)候，計(jì)算機(jī)內(nèi)存可能無法容納被操作的矩陣，這給直接法帶來很大的挑戰(zhàn)。而對(duì)于迭代法，則可以將矩陣的某一部分讀入內(nèi)存進(jìn)行操作，然后再操作另外部分。本文使用上述四種算法求解對(duì)應(yīng)的方程組，驗(yàn)證各種算法的精確度和計(jì)算速度。1 算法介紹1.1 雅克比迭代法1.1.1 算法理論設(shè)線性方程組 (1) 的系數(shù)矩陣A可逆且主對(duì)角元素 均不為零,令 并將A分解成 (2)從而(1)可寫成 令 其中 . (3)以B1為迭代矩陣的迭代法(公式) (4)稱為雅克比(Jacbi)迭代法(公式),用向量的分量來表示,(4)為 (5)其中為初

4、始向量.1.1.2 算法描述1給定迭代初始向量0以及誤差要求delta2根據(jù)雅克比迭代公式計(jì)算出下一組向量3判斷是否滿足誤差要求，即|k+1 k| < delta4若誤差滿足要求，則停止迭代返回結(jié)果；若否，則返回第二步進(jìn)行下一輪迭代1.2 高斯賽德爾迭代法1.2.1 算法理論由雅克比迭代公式可知,在迭代的每一步計(jì)算過程中是用的全部分量來計(jì)算的所有分量,顯然在計(jì)算第i個(gè)分量時(shí),已經(jīng)計(jì)算出的最新分量沒有被利用，從直觀上看,最新計(jì)算出的分量可能比舊的分量要好些.因此，對(duì)這些最新計(jì)算出來的第次近似的分量加以利用,就得到所謂解方程組的高斯塞德爾（Gauss-Seidel）迭代法.把矩陣A分解成 (

5、6) 其中,分別為的主對(duì)角元除外的下三角和上三角部分,于是,方程組(1)便可以寫成 即其中 (7)以為迭代矩陣構(gòu)成的迭代法(公式) (8)稱為高斯塞德爾迭代法(公式),用變量表示的形式為 (9)1.2.2 算法描述1給定迭代初始向量0以及誤差要求delta2根據(jù)高斯賽德爾迭代公式計(jì)算出下一組向量3判斷是否滿足誤差要求，即|k+1 k| < delta4若誤差滿足要求，則停止迭代返回結(jié)果；若否，則返回第二步進(jìn)行下一輪迭代1.3 高斯消去法1.3.1 算法理論下面三種變換稱為初等行變換：1.對(duì)調(diào)兩行；2.以數(shù)k0乘某一行中的所有元素；3.把某一行所有元素的k倍加到另一行對(duì)應(yīng)的元素上去。 利用

6、初等變換將增廣矩陣矩陣化為階梯形式，此階梯矩陣所代表的方程組與原方程等價(jià)，然后利用回代法求此階梯矩陣的解，與原方程解相同。在轉(zhuǎn)化為階梯矩陣的過程中可以使用選擇列主元的方式減小誤差。1.3.2 算法描述以4階為例：第1步消元在增廣矩陣（A，b）第一列中找到絕對(duì)值最大的元素，將其所在行與第一行交換，再對(duì)（A，b）做初等行變換使原方程組轉(zhuǎn)化為如下形式：第2步消元在增廣矩陣（A，b）中的第二列中（從第二行開始）找到絕對(duì)值最大的元素，將其所在行與第二行交換，再對(duì)（A，b）做初等行變換使原方程組轉(zhuǎn)化為：第3步消元在增廣矩陣（A，b）中的第三列中（從第三行開始）找到絕對(duì)值最大的元素，將其所在行與第二行交換，

7、再對(duì)（A，b）做初等行變換使原方程組轉(zhuǎn)化為：按4 à 3à 2à 1 的順序回代求解出方程組的解1.4 改進(jìn)的平方根法1.4.1 算法理論當(dāng)方程的系數(shù)矩陣為對(duì)稱正定時(shí)，可以直接做高斯消去法。也就是說對(duì)稱陣正定矩陣保證能直接作LU分解。由LU分解公式：uli = ali (i = 1,2,3, ,n)lil = ail / all (i = 1,2,3, ,n)因?yàn)锳對(duì)稱： ail = ali (i = 1,2,3, ,n)所以：lil = ali / ull = uli / ull (i = 1,2,3, ,n)綜上所述的如下結(jié)論：若A為對(duì)稱正定矩陣，則A一定能作

8、LU分解，且：lil = uli / ull (k = 1,2,3, , n-1; i = k+1, k+2, k+3, ,n)亦即若A為對(duì)稱正定矩陣，則其單位下三角矩陣不必按正常LU分解公式求得，而只需將求得的U的第k行元素除以u(píng)kk 即得相應(yīng)L的第k列元素。1.4.2算法描述對(duì)于給定的線性方程組A=b，首先對(duì)A進(jìn)行LU分解。根據(jù)公式：uki=aki-q=1nlkquqj計(jì)算出U矩陣的第一行，然后根據(jù)lil=uliull計(jì)算L矩陣第一列，而后以同樣的方法依次將L和U矩陣所有行和列的非零元算出得：LU=b利用前向替換法對(duì)方程組LY=b求解Y，然后利用回代法對(duì)方程組U=Y求解。2 結(jié)果分析2.1

9、 雅克比迭代法現(xiàn)在用雅克比迭代法計(jì)算方程組A=b，用以驗(yàn)證其收斂速度以及計(jì)算準(zhǔn)確性其中：A= b=給定精度為delta = 110-5，迭代初值為3.0, 3.0, 3.0, 3.0，則得到如下表的雅克比迭代法的計(jì)算過程，其中err為所估計(jì)的誤差。迭代次數(shù)共11次，計(jì)算最終結(jié)果為-1.467403, -2.358649, 0.657596, 2.842408。可以看出，雅克比迭代所計(jì)算出的解與預(yù)測(cè)相符，但為了達(dá)到精度要求所需迭代次數(shù)較多。表1雅克比迭代法計(jì)算結(jié)果n1234err0-1.400000-2.1250000.3000002.0000007.3427251-1.372500-2.087

10、5000.6675002.7522730.8385302-1.442250-2.3236650.6657502.7779540.2475913-1.465516-2.3335140.6560832.8358630.0639154-1.464568-2.3564380.6597522.8355550.0232375-1.467594-2.3558640.6572702.8422270.0077566-1.467040-2.3586760.6579322.8415700.0030147-1.467454-2.3583470.6575402.8424470.0010968-1.467342-2.35

11、87250.6576562.8422840.0004419-1.467403-2.3586490.6575962.8424080.00016810-1.467384-2.3587030.6576162.8423760.00006811-1.467394-2.3586900.6576062.8423950.00002712-1.467390-2.3586970.6576092.8423900.00001113-1.467392-2.3586950.6576082.8423920.0000042.2 高斯賽德爾迭代法用高斯賽德爾迭代法進(jìn)行計(jì)算，使用與雅克比迭代法相同的數(shù)據(jù)（方程組，初始值和精度要求

12、），得到計(jì)算結(jié)果為-1.467391, -2.358696, 0.657609, 2.842391，與雅克比迭代法所計(jì)算出的結(jié)果相同。但是，高斯賽德爾迭代法所用的迭代次數(shù)明顯比較少，可見高斯賽德爾法是一種更加高效的計(jì)算方法。如下表為高斯賽德爾迭代法的計(jì)算過程：表2高斯賽德爾迭代法計(jì)算結(jié)果n1234err0-1.4-2.1250.6675002.7856817.1492731-1.446000-2.3361930.6555802.8380140.2227092-1.464735-2.3573080.6572162.8422190.0285873-1.467174-2.3586800.657566

13、2.8424030.0028264-1.467381-2.3587050.6576062.8423950.0002135-1.467392-2.3586970.6576092.8423920.0000136-1.467391-2.3586960.6576092.8423910.0000012.3 列主元高斯消去法用列主元高斯消去法計(jì)算線性方程組A=b,消去時(shí)先選擇列主元以減小誤差，其中：A= b=計(jì)算結(jié)果為-3.413793,5.206897,1.448275，與預(yù)測(cè)結(jié)果相符。2.4 改進(jìn)平方根法用改進(jìn)平方根發(fā)計(jì)算線性方程組A=b，當(dāng)然A為對(duì)稱正定矩陣：A= b=計(jì)算結(jié)果為2.0,1.0,-1

14、.0，計(jì)算結(jié)果與預(yù)測(cè)值相同。3 結(jié)論雅克比迭代法和高斯賽德爾迭代法，作為常用的求解線性方程組的迭代算法，都具有思路簡單，容易實(shí)現(xiàn)的特點(diǎn)，且都有很好的收斂特性。但相比之下，高斯賽德爾迭代法作為雅克比迭代法的改進(jìn)，逼近速度更快，需要更少的迭代次數(shù)，所以應(yīng)用廣泛。高斯消去法和改進(jìn)的平方根法則具有各自的特點(diǎn)。高斯消去法運(yùn)算復(fù)雜程度較高，但作為通用解法有著不可忽視的重要性；相反改進(jìn)的平方根法雖然所需計(jì)算量小，但只能計(jì)算其系數(shù)矩陣為對(duì)稱正定的情況，其使用范圍受到了嚴(yán)格的限制。參考文獻(xiàn)： 1Jhn H.Mathews, Kurtis D.Fink. 數(shù)值方法M. 北京：電子工業(yè)出版社，2010.2孫志忠,

15、吳宏偉, 袁慰平, 聞?wù)鸪? 計(jì)算方法與實(shí)習(xí)M. 第5版. 南京:東南大學(xué)出版社, 2011.附錄（pythn語言代碼）雅克比和高斯賽德爾迭代法程序：frm math imprt frm cpy imprt A = 10.0, -1.0, 2.0 , 0.0 , 0.0 , 8.0, -1.0 , 3.0 , 2.0 , -1.0, 10.0, 0.0 , -1.0, 3.0, -1.0 , 11.0B = -11.0, -11.0, 6.0, 25.0P = 3.0, 3.0, 3.0, 3.0 the initial valuesma1 = 1000 ma iteratin timedel

16、ta = 1e-15 precisin demanddef distance(A, B): dstance = 0 fr i in range(len(A): dstance = dstance + (Ai - Bi) 2 dstance = sqrt(dstance) return dstancedef jacbi(A, B, P, delta, ma1): = cpy(P) N = len(B) fr k in range(ma1): fr i in range(N): temp = 0 fr j in range(N): temp = temp + Aij Pj temp = temp

17、- Aii Pi i = (Bi - temp) / Aii err = distance(, P) P = cpy() print(k, , err) if(err < delta): return return Nnedef gseid(A, B, P, delta, ma1): = cpy(P) N = len(B) fr k in range(ma1): fr i in range(N): temp = 0 fr j in range(N): temp = temp + Aij j temp = temp - Aii i i = (Bi - temp) / Aii err = d

18、istance(, P) P = cpy() print(k, , err) if(err < delta): return return Nneif name = 'main': print(jacbi(A, B, P, delta, ma1) print(gseid(A, B, P, delta, ma1)列主元高斯消去法程序：frm cpy imprt A = 2.0, -1.0, 9.0, 4.0, 2.0, 5.0, 1.0, 2.0, 0.0b = 1.0, 4.0, 7.0def eliminatin(A, B, inde): scalar = Binde

19、/ Ainde fr i in range(len(A): Bi = Bi - Ai scalardef gauss(A, b): M = deepcpy(A) length = len(A) = 0.0 length generate the augmented matra fr i in range(length): Mi.append(bi) transfrm t an upper trangular matri fr i in range(length-1): fr j in range(i, length): if Mji > Mii: Mj, Mi = Mi, Mj fr j in range(i+1, length): eliminatin(Mi, Mj, i) rlling back t the answers fr i in range(length-1, -1, -1): i = Milengt