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

1、線性方程組的四種數(shù)值解法摘要：本文介紹了四種求解線性方程組的數(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)鍵詞：雅克比迭代法；高斯賽德爾迭代法；高斯消去法；改進(jìn)的平方根法；線性方程組引言算法介

2、紹雅克比迭代法算法理論設(shè)線性方程組Ax b (1)A 11,a22,.,annD g 11,a22,.,ann 并將A 分解成A A D D (2)從而(1)可寫成令Dx D Ax bx 1x f1其中11B I D1,f11D.(3)以 B1 為迭代矩陣的迭代法(公式)1xk1 Bxk 11 (4)稱為雅克比(Jacobi)迭代法(公式),用向量的分量來表示,(4)為ix ( k ) i1a njjiai ( k )xjxik(5)12其中x0 x0,x0.x012其中為初始向量.算法描述X0delta根據(jù)雅克比迭代公式計(jì)算出下一組向量判斷X是否滿足誤差要求，即|Xk+1 Xk|delta若

3、誤差滿足要求，則停止迭代返回結(jié)果；若否，則返回第二步進(jìn)行下一輪迭代高斯賽德爾迭代法算法理論xx由雅克比迭代公式可知在迭代的每一步計(jì)算過程中是用的全部分量來計(jì)算x1x1k 1的所有分,顯然在計(jì)算第i個(gè)分量 i時(shí),已經(jīng)計(jì)算出的最新分量1i1沒有被利用，從直觀上看,最新計(jì)算出的分量可能比舊的分量要好些.因此，對(duì)這些最新計(jì)算出來k 1的第k次近似xkxj的分量加以利用,就得到所謂解方程組的高斯塞德爾（Gauss-Seidel）迭代法.A 分解成A D L U (6)其中D g 11,a22,.,ann, LU A(1)便可以寫成D Lx Ux b即x B2x f2其中2B DL1U,22f DL12以

4、B22xk1 B xk 22 (8)x稱為高斯塞德爾迭代法(公式),用變量表示的形式為xx( k ) iii1x(k) jjn ( k ) j11j1ji1算法描述ik (9)X0delta根據(jù)高斯賽德爾迭代公式計(jì)算出下一組向量判斷X是否滿足誤差要求，即|Xk+1Xk|delta若誤差滿足要求，則停止迭代返回結(jié)果；若否，則返回第二步進(jìn)行下一輪迭代高斯消去法算法理論下面三種變換稱為初等行變換： 1.對(duì)調(diào)兩行；k0乘某一行中的所有元素；k倍加到另一行對(duì)應(yīng)的元素上去。算法描述4 階為例：a11 x1aaa x b 212n 2 2 aaa xb n2nn n n第 1 步消元在增廣矩陣（A，b）第一

5、列中找到絕對(duì)值最大的元素，將其所在行與第一行交換，再對(duì)（A，b）做初等行變換使原方程組轉(zhuǎn)化為如下形式：* * x10* * x 2 0* * x30* * 4 第 2 步消元在增廣矩陣（A，b）中的第二列中（從第二行開始）找到絕對(duì)值最大的元素，將其所在行與第二行交換，再對(duì)（A，b）做初等行變換使原方程組轉(zhuǎn)化為：* x10* x 2 00* x300* x 4 第 3 步消元在增廣矩陣（A，b）中的第三列中（從第三行開始）找到絕對(duì)值最大的元素，將其所在行與第二行交換，再對(duì)（A，b）做初等行變換使原方程組轉(zhuǎn)化為：* x10*x 2 00* x3000 x 4 按 x4 x 3x 2x 1 的順序回

6、代求解出方程組的解改進(jìn)的平方根法算法理論分解。由分解公式： =(i = 1,2,3, ,n) = /(i = 1,2,3, ,n)因?yàn)?A 對(duì)稱：=(i = 1,2,3, ,n)所以：=/=/ (i=1,2,3,n)綜上所述的如下結(jié)論：若為對(duì)稱正定矩陣，則一定能作 LU 分解，且：=/(k = 1,2,3, n-1; i = k+1, k+2, k+3,n)kk Lk算法描述對(duì)于給定的線性方程組 = b，首先對(duì)進(jìn)行分解。根據(jù)公式： = =1計(jì)算出矩陣的第一行，然后根據(jù)= 計(jì)算矩陣第一列，而后以同樣的方法依次將 L 和 U 矩陣所有行和列的非零元算出得： = b利用前向替換法對(duì)方程組 = 求解，

7、然后利用回代法對(duì)方程組 = 求解。結(jié)果分析雅克比迭代法現(xiàn)在用雅克比迭代法計(jì)算方程組 = b，用以驗(yàn)證其收斂速度以及計(jì)算準(zhǔn)確性其中：1012011 081311 = b = 21 1006131 11 -5，迭代初值為3.03.03.03.0，則得到如下表的雅克比迭代法的err 11-1.467403, -2.358649, 0.657596, 2.842408?？梢钥闯觯趴吮鹊?jì)算出的解與預(yù)測(cè)相符，但為了達(dá)到精度要求所需迭代次數(shù)較多。表 1 雅克比迭代法計(jì)算結(jié)果nx1x2x3x4err0-1.400000-2.1250000.3000002.0000007.3427251-1.37250

8、0-2.0875000.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.4673

9、42-2.3587250.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.000004高斯賽德爾迭代法, , , 是一種更加高效的計(jì)算方法。如下表為高斯賽德爾迭代法的計(jì)算過程：表

10、 2 高斯賽德爾迭代法計(jì)算結(jié)果nx1x2x3x4err0-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.6575662.8424030.0028264-1.467381-2.3587050.6576062.8423950.0002135-1.467392-2.3586970.6576092.8423920.0000136-1.467391-2.3586960.65

11、76092.8423910.000001列主元高斯消去法用列主元高斯消去法計(jì)算線性方程組 = ,消去時(shí)先選擇列主元以減小誤差，其中：219 =42 = 4 1207計(jì)算結(jié)果為-3.413793,5.206897,1.448275，與預(yù)測(cè)結(jié)果相符。改進(jìn)平方根法用改進(jìn)平方根發(fā)計(jì)算線性方程組 = b，當(dāng)然為對(duì)稱正定矩陣： 424 =21710 = 3 41097計(jì)算結(jié)果為2.0,1.0,-1.0，計(jì)算結(jié)果與預(yù)測(cè)值相同。結(jié)論的情況，其使用范圍受到了嚴(yán)格的限制。參考文獻(xiàn)：JohnH.Mathews,Kurtis D.Fink. 數(shù)值方法M. 北京：電子工業(yè)出版社，2010.孫志忠, 吳宏偉, 袁慰平,

12、聞?wù)鸪? 計(jì)算方法與實(shí)習(xí)M. 5 版. 南京:東南大學(xué)出版社2011.附錄（python 語言代碼）雅克比和高斯賽德爾迭代法程序：from math import * from copy import *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 valuesmax1 = 1000# max iteration timedelta = 1e-15# pr

13、ecisiondemanddef distance(A, B): dstance = 0for i in range(len(A):dstance = dstance + (Ai - Bi) * 2dstance = sqrt(dstance) return dstancedef jacobi(A, B, P, delta, max1): X = copy(P)N = len(B)for k in range(max1): for i in range(N):temp = 0for j in range(N):temp=temp+Aij*Pj temp=temp-Aii*PiXi = (Bi

14、- temp) / Aii err = distance(X, P)P = copy(X)# print(k, X, err) if(err delta):return X return Nonedefgseid(A,B,P,delta,max1): X =copy(P)N = len(B)for k in range(max1): for i in range(N):temp = 0for j in range(N):temp=temp+Aij*Xj temp=temp-Aii*XiXi = (Bi - temp) / Aii err = distance(X, P)P = copy(X)#

15、 print(k, X, err) if(err Mii:Mj, Mi = Mi, Mjfor j in range(i+1, length): elimination(Mi, Mj, i)# rolling back to the answers for i in range(length-1, -1, -1):xi = Milengthfor j in range(i+1, length): xi = xi - Mij * xi = xi / Mii return xif name = main : print(gauss(A, b)改進(jìn)平方根法程序from copy import *from copy import *A = 4.0, -2.0 , -4.0 ,-2.0, 17.0, 10.0,-4.0, 10.0, 9.0b = 10.0, 3.0, -7.0def imprv_sqrt(A, b): length = len(A) M =deepcopy(A)# the decomposition for k in range(length):for j in range(k, length): for q in range(k):Mkj-=Mkq*Mqjif j != k:Mjk = Mkj / Mkk # roll