线性方程组的四种数值解法

1、线性方程组的四种数值解法 (电子科技大学物理电子学院，四川 成都 610054)摘要：本文介绍了四种求解线性方程组的数值解法: 雅克比迭代法、高斯赛德尔迭代法、高斯消去法和改进的平方根法的基本原理和算法流程，通过求解具体方程，对四种求解方法进行了对比。对于雅克比迭代法和高斯赛德尔迭代法，研究了两种算法对求解同一方程组的迭代效率差异，结果表明高斯赛德尔迭代法达到同样精度所需迭代次数较少。对于高斯消去法，通过选择列主元的方法提高算法的准确度，计算结果表明高斯消去法计算精确，且运算复杂度也不是很高。对于改进的平方根法，其运算复杂度低，但对于给定的方程组有着严苛的要求。关键词：雅克比迭代法；高斯赛德尔

2、迭代法；高斯消去法；改进的平方根法；线性方程组引言线性方程组的求解在日常生活和科研中有着极其重要的应用，但在实际运算中，当矩阵的维数较高时，用初等方法求解的计算复杂度随维数的增长非常快，因此，用数值方法求解线性方程组的重要性便显现出来。经典的求解线性方程组的方法一般分为两类：直接法和迭代法。前者例如高斯消去法,改进的平方根法等，后者的例子包括雅克比迭代法，高斯赛德尔迭代法等。这些方法的计算复杂度在可以接受的范围内，因此被广泛采用。一般来说，直接法对于阶数比较低的方程组比较有效；而后者对于比较大的方程组更有效。在实际计算中，几十万甚至几百万个未知数的方程组并不少见。在这些情况下，迭代法有无可比拟

3、的优势。另外，使用迭代法可以根据不同的精度要求选择终止时间，因此比较灵活。在问题特别大的时候，计算机内存可能无法容纳被操作的矩阵，这给直接法带来很大的挑战。而对于迭代法，则可以将矩阵的某一部分读入内存进行操作，然后再操作另外部分。本文使用上述四种算法求解对应的方程组，验证各种算法的精确度和计算速度。1 算法介绍1.1 雅克比迭代法1.1.1 算法理论设线性方程组 (1) 的系数矩阵A可逆且主对角元素 均不为零,令 并将A分解成 (2)从而(1)可写成 令 其中 . (3)以B1为迭代矩阵的迭代法(公式) (4)称为雅克比(Jacbi)迭代法(公式),用向量的分量来表示,(4)为 (5)其中为初

4、始向量.1.1.2 算法描述1给定迭代初始向量0以及误差要求delta2根据雅克比迭代公式计算出下一组向量3判断是否满足误差要求，即|k+1 k| < delta4若误差满足要求，则停止迭代返回结果；若否，则返回第二步进行下一轮迭代1.2 高斯赛德尔迭代法1.2.1 算法理论由雅克比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i个分量时,已经计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯塞德尔（Gauss-Seidel）迭代法.把矩阵A分解成 (

5、6) 其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 即其中 (7)以为迭代矩阵构成的迭代法(公式) (8)称为高斯塞德尔迭代法(公式),用变量表示的形式为 (9)1.2.2 算法描述1给定迭代初始向量0以及误差要求delta2根据高斯赛德尔迭代公式计算出下一组向量3判断是否满足误差要求，即|k+1 k| < delta4若误差满足要求，则停止迭代返回结果；若否，则返回第二步进行下一轮迭代1.3 高斯消去法1.3.1 算法理论下面三种变换称为初等行变换：1.对调两行；2.以数k0乘某一行中的所有元素；3.把某一行所有元素的k倍加到另一行对应的元素上去。 利用

6、初等变换将增广矩阵矩阵化为阶梯形式，此阶梯矩阵所代表的方程组与原方程等价，然后利用回代法求此阶梯矩阵的解，与原方程解相同。在转化为阶梯矩阵的过程中可以使用选择列主元的方式减小误差。1.3.2 算法描述以4阶为例：第1步消元在增广矩阵（A，b）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A，b）做初等行变换使原方程组转化为如下形式：第2步消元在增广矩阵（A，b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：第3步消元在增广矩阵（A，b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，

7、再对（A，b）做初等行变换使原方程组转化为：按4 à 3à 2à 1 的顺序回代求解出方程组的解1.4 改进的平方根法1.4.1 算法理论当方程的系数矩阵为对称正定时，可以直接做高斯消去法。也就是说对称阵正定矩阵保证能直接作LU分解。由LU分解公式：uli = ali (i = 1,2,3, ,n)lil = ail / all (i = 1,2,3, ,n)因为A对称： ail = ali (i = 1,2,3, ,n)所以：lil = ali / ull = uli / ull (i = 1,2,3, ,n)综上所述的如下结论：若A为对称正定矩阵，则A一定能作

8、LU分解，且：lil = uli / ull (k = 1,2,3, , n-1; i = k+1, k+2, k+3, ,n)亦即若A为对称正定矩阵，则其单位下三角矩阵不必按正常LU分解公式求得，而只需将求得的U的第k行元素除以ukk 即得相应L的第k列元素。1.4.2算法描述对于给定的线性方程组A=b，首先对A进行LU分解。根据公式：uki=aki-q=1nlkquqj计算出U矩阵的第一行，然后根据lil=uliull计算L矩阵第一列，而后以同样的方法依次将L和U矩阵所有行和列的非零元算出得：LU=b利用前向替换法对方程组LY=b求解Y，然后利用回代法对方程组U=Y求解。2 结果分析2.1

9、 雅克比迭代法现在用雅克比迭代法计算方程组A=b，用以验证其收敛速度以及计算准确性其中：A= b=给定精度为delta = 110-5，迭代初值为3.0, 3.0, 3.0, 3.0，则得到如下表的雅克比迭代法的计算过程，其中err为所估计的误差。迭代次数共11次，计算最终结果为-1.467403, -2.358649, 0.657596, 2.842408。可以看出，雅克比迭代所计算出的解与预测相符，但为了达到精度要求所需迭代次数较多。表1雅克比迭代法计算结果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 高斯赛德尔迭代法用高斯赛德尔迭代法进行计算，使用与雅克比迭代法相同的数据（方程组，初始值和精度要求

12、），得到计算结果为-1.467391, -2.358696, 0.657609, 2.842391，与雅克比迭代法所计算出的结果相同。但是，高斯赛德尔迭代法所用的迭代次数明显比较少，可见高斯赛德尔法是一种更加高效的计算方法。如下表为高斯赛德尔迭代法的计算过程：表2高斯赛德尔迭代法计算结果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 列主元高斯消去法用列主元高斯消去法计算线性方程组A=b,消去时先选择列主元以减小误差，其中：A= b=计算结果为-3.413793,5.206897,1.448275，与预测结果相符。2.4 改进平方根法用改进平方根发计算线性方程组A=b，当然A为对称正定矩阵：A= b=计算结果为2.0,1.0,-1

14、.0，计算结果与预测值相同。3 结论雅克比迭代法和高斯赛德尔迭代法，作为常用的求解线性方程组的迭代算法，都具有思路简单，容易实现的特点，且都有很好的收敛特性。但相比之下，高斯赛德尔迭代法作为雅克比迭代法的改进，逼近速度更快，需要更少的迭代次数，所以应用广泛。高斯消去法和改进的平方根法则具有各自的特点。高斯消去法运算复杂程度较高，但作为通用解法有着不可忽视的重要性；相反改进的平方根法虽然所需计算量小，但只能计算其系数矩阵为对称正定的情况，其使用范围受到了严格的限制。参考文献： 1Jhn H.Mathews, Kurtis D.Fink. 数值方法M. 北京：电子工业出版社，2010.2孙志忠,

15、吴宏伟, 袁慰平, 闻震初. 计算方法与实习M. 第5版. 南京:东南大学出版社, 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