第三章线性代数方程组的数值解

1、一般地说，差分格式被写成如下的形式：一般地说，差分格式被写成如下的形式： nnnnnnnnnnbTaTaTabTaTaTabTaTaTa22112222212111212111（3.1） 其中，其中，n是节点数，也即方程个数，每个方是节点数，也即方程个数，每个方程对应一个节点。程对应一个节点。 和和bi（i=1，2，n；j=1，2，n）都是常数）都是常数 矩阵，矩阵，ija 第三章：线性代数方程组的数值解法第三章：线性代数方程组的数值解法1 引言引言方程组（方程组（3.1）可被进一步写成矩阵的形式）可被进一步写成矩阵的形式 : , I JIATB（3.2） 111211212222,12nnI

2、 JnnnnnaaaTaaaTATTaaa nIbbbB21其中其中n解线性代数方程组的解线性代数方程组的直接法直接法。n大家知道，线性方程组大家知道，线性方程组 BTA（3.3） 中只要矩阵的行列式中只要矩阵的行列式 ，方程组（，方程组（3.3）就有唯一解，其表达式为就有唯一解，其表达式为:0detAnjAATjj, 2 , 1det/det（3.4） 其中其中 为用右端向量为用右端向量 替换行列式替换行列式 的第的第 j 列而得的，这一公式为著名的列而得的，这一公式为著名的Cramer法则。法则。 jAdet BAdetn显然，按显然，按Cramer法则求解方程组（法则求解方程组（3.3）

3、需要计算）需要计算n+1个个n阶行列式。每个阶行列式。每个n阶行列式按直接展开办法阶行列式按直接展开办法来算，需作（来算，需作（n-1）n！次乘法和！次乘法和n!次加法运算。次加法运算。当当n=30时，共约需完成时，共约需完成 次乘法和加法运算，这次乘法和加法运算，这是一个十分惊人的数字，即使在一台每秒作一亿次是一个十分惊人的数字，即使在一台每秒作一亿次运算的计算机上完成这一计算也是不可能的。所以，运算的计算机上完成这一计算也是不可能的。所以，尽管这种办法也是一种直接法，并且理论上可行，尽管这种办法也是一种直接法，并且理论上可行，但实际上是无法进行求解的。即使采用其它办法来但实际上是无法进行求

4、解的。即使采用其它办法来计算行列式，按计算行列式，按Cramer法则求解的工作量也比通法则求解的工作量也比通常的直接法大得多。因而，常的直接法大得多。因而，Cramer法则对于数值法则对于数值计算来说是没有什么用处的，仅在一些特殊场合才计算来说是没有什么用处的，仅在一些特殊场合才有用。有用。3310n另外，矩阵求逆也是人们熟悉的求解线性方程另外，矩阵求逆也是人们熟悉的求解线性方程的一种方法。的一种方法。1BAT(3.5) 但求逆矩阵但求逆矩阵 时，要计算时，要计算 个个( n -1)阶行列式，阶行列式，和一个和一个n 阶行列式，它的计算工作量也是相当可阶行列式，它的计算工作量也是相当可观的，对

5、于观的，对于n 较大的情况也没有什么现实意义。较大的情况也没有什么现实意义。 2n1An现在计算机上常用的现在计算机上常用的直接解法直接解法大多数是以大多数是以系数矩阵的系数矩阵的三角形化三角形化为基础的。就是说，为基础的。就是说，先对方程组进行变换，使其化为等价的先对方程组进行变换，使其化为等价的（即具有相同解的）三角形方程组。由于（即具有相同解的）三角形方程组。由于三角形方程组三角形方程组的求解十分容易，原方程的的求解十分容易，原方程的求解问题即告解决。求解问题即告解决。为讨论方便起见，我们首先叙述为讨论方便起见，我们首先叙述三角形方程的解法三角形方程的解法，然后讨论将原方程化为等价三角形

6、方程组的方法。然后讨论将原方程化为等价三角形方程组的方法。由于计算机的字长是有限的，每次运算之后还要由于计算机的字长是有限的，每次运算之后还要对结果进行舍入，所以，虽然理论上直接法在有对结果进行舍入，所以，虽然理论上直接法在有限步内可以得到精确解，但计算机上实际得到的只限步内可以得到精确解，但计算机上实际得到的只是是近似解近似解。 2 三角形方程组的解法三角形方程组的解法 所谓三角形方程组是指下面两种形式的方程组所谓三角形方程组是指下面两种形式的方程组：(3.6) nnnnnnbTlTlTlbTlTlTlbTlTlbTl.2211333323213122221211111方程组方程组(3.6)

7、叫作叫作下三角形方程组下三角形方程组或或(3.7) 方程组方程组(3.7)叫作叫作上三角形方程组上三角形方程组。 nnnnnnnnnnnnnnndTudTuTudTuTuTudTuTuTuTu1, 111, 12232322211313212111.n若用矩阵符号可分别写为：若用矩阵符号可分别写为：,DTUBTL其中其中L为方程组为方程组(3.6)的系数所构成的下三角的系数所构成的下三角形矩阵，其元素满足关系：形矩阵，其元素满足关系：jilij 0U为方程组为方程组(2.7)的系数所构成的上三角形的系数所构成的上三角形矩阵，其元素满足关系：矩阵，其元素满足关系： jiuij 0n三角形方程组的

8、求解是很简单的。方程组三角形方程组的求解是很简单的。方程组(3.6)的的计算公式可归结为：计算公式可归结为：)(, 3, 2/ ),(/1122111111jinilTlTlTlbTlbTiiiiiiiii(3.8) 这个计算过程通常也只作这个计算过程通常也只作前推前推过程。过程。n对于方程组对于方程组(3.7)，其计算公式可归结为：，其计算公式可归结为：/,(,)/, 1 1, 221 ,2, ,1 ()Td unn nnTd uTuTu Tuiii iii iiin niii nnn j i (3.9) 这个计算过程通常也叫作这个计算过程通常也叫作回代方程回代方程。由以上分析可以看到，只要

9、把方程组化成了等由以上分析可以看到，只要把方程组化成了等价的三角形方程组，求解就容易了。价的三角形方程组，求解就容易了。 结论：结论：3 Gauss消去法消去法 Gauss消去法（简称消去法）的提出已有相当长消去法（简称消去法）的提出已有相当长的时间了，是一种古老的方法。然而，近年来在的时间了，是一种古老的方法。然而，近年来在计算机上求解线性代数方程组的实践表明，它仍计算机上求解线性代数方程组的实践表明，它仍是是直接法直接法中最常用的一种方法，也是最有效的方中最常用的一种方法，也是最有效的方法之一。其基本思想是：法之一。其基本思想是：用逐次消去一个未知数用逐次消去一个未知数的办法把原来的方程组

10、化为等价的（具有相同解）的办法把原来的方程组化为等价的（具有相同解）三角形方程组三角形方程组。这样，求解就很容易了。这样，求解就很容易了。系数矩阵的系数矩阵的三角形化三角形化n假定把要求的假定把要求的n 阶线性的方程组阶线性的方程组(3.1)改写成如改写成如下形式：下形式： ) 1 () 1 (3) 1 (32) 1 (21) 1 (1) 1 (2) 1 (23) 1 (232) 1 (221) 1 (21) 1 (1) 1 (13) 1 (132) 1 (121) 1 (11nnnnnnnnnnnbTaTaTaTabTaTaTaTabTaTaTaTa(3.10) 用矩阵符号记为用矩阵符号记为

11、)1()1(BTAn其中其中 为为 方阵，方阵， 为为 向量，它们分向量，它们分别为：别为：分别从原方程组的第二个方程减去第一个方分别从原方程组的第二个方程减去第一个方程乘以程乘以 ，第三个方程减去第一个方程，第三个方程减去第一个方程乘以乘以 ，如此等等，即可消去后面，如此等等，即可消去后面n -1个方程中的未知量个方程中的未知量 。 ) 1 (11) 1 (21/aa)1(11)1(31/aa1T)1(Ann)1(B1n) 1 () 1 (2) 1 (1) 1 () 1 () 1 (3) 1 (2) 1 (1) 1 (2) 1 (23) 1 (22) 1 (21) 1 (1) 1 (13)

12、1 (12) 1 (11) 1 (,nnnnnnnnbbbBaaaaaaaaaaaaAn这时方程组这时方程组(3.10)就变为如下等价方程组：就变为如下等价方程组： (1)(1)(1)(1)(1)12311121311(2)(2)(2)(2)0_23222322(2)(2)(2)(2)0_233233330_(2)(2)(2)(2)0_2323aTaTaTaTbnnaTaTaTbnnaTaTaTbnnaTaTaTbnn nnnnn表示成矩阵形式为：表示成矩阵形式为：)2()2(BTA其中其中)2()2(3)2(2) 1 (1)2()2()2(3)2(2)2(3)2(33)2(23)2(2)2(

13、23)2(22) 1 (1) 1 (13) 1 (12) 1 (11)2(b,000nnnnnnnnbbbBaaaaaaaaaaaaaAn若令若令 ，则系数，则系数 和和 的计的计 算公式应为：算公式应为：) 1 (11) 1 (11/aamaii)2(ija)2(ibnjibmbbamaaiiijiijij, 3, 2,) 1 (11) 1 ()2() 1 (11) 1 ()2(类似地，分别从上述等价方程组的第三个方程减类似地，分别从上述等价方程组的第三个方程减去第二个方程乘以去第二个方程乘以 ，第四个方程减去第，第四个方程减去第二个方程乘以二个方程乘以 ，如此等等，即可进一步，如此等等，即

14、可进一步消去后面消去后面n -2个方程中的未知量个方程中的未知量T2，而将方程，而将方程(3.10)变为如下等价形式：变为如下等价形式： )2(22)2(32/ aa)2(22)2(42/ aa)3()3(3)3(3)3(44)3(43)3(43)3(33)3(33)3(33)2(2)2(23)2(232)2(22)1(1)1(13)1(132)1(121)1(11nnnnnnnnnnnbTaTabTaTabTaTabTaTaTabTaTaTaTa表示成矩阵形式为：表示成矩阵形式为：)3()3(BTAn其中其中)3()3(4)3(3)2(2) 1 (1)3()3()3(3)3(4)3(43)3

15、(3)3(33)2(2)2(23)2(22) 1 (1) 1 (13) 1 (12) 1 (11)3(,000000nnnnnnnnbbbbbBaaaaaaaaaaaaaA若令若令 ，则系数，则系数 和和 的计算的计算公式应为：公式应为：)2(22)2(22/aamii)3(ija)3(ibnjibmbbamaaiiijiijij, 4, 3,)2(22)2() 3()2(22)2() 3(n上述的消去步骤还可以进行下去。如此继续之，上述的消去步骤还可以进行下去。如此继续之，重复上述步骤重复上述步骤(n -1)次以后，我们即可得到如次以后，我们即可得到如下等价三角形方程组：下等价三角形方程组：

16、)()()3(3)3(33)3(33)2(2)2(23)2(232)2(22) 1 (1) 1 (13) 1 (132) 1 (121) 1 (11nnnnnnnnnnnnbTabTaTabTaTaTabTaTaTaTa(3.11) 表示为矩阵形式：表示为矩阵形式：)()(nnBTA其中其中 为如下上三角形矩阵，为如下上三角形矩阵， 为为 向量；向量；)(nA)(nB1n)() 3(3) 2(2) 1 (1)()() 3(3) 3(33) 2(2) 2(23) 2(22) 1 (1) 1 (13) 1 (13) 1 (11)(,nnnnnnnnnnnbbbbBaaaaaaaaaaA0(3.12

17、) n三角形方程组三角形方程组 很容易用前述的回很容易用前述的回代过程代过程(3.9)求解，这样就完成了消去法求解求解，这样就完成了消去法求解n阶阶线性代数方程组的过程。从原来方程组线性代数方程组的过程。从原来方程组(3.10)得得出等价三角形方程组出等价三角形方程组(3.11)的过程称之为消去过的过程称之为消去过程。采用前面的记号，我们可将消去过程的计算程。采用前面的记号，我们可将消去过程的计算公式归结为对于公式归结为对于 ，递推地计，递推地计算如下各量：算如下各量：)()(nnBTA1n21k,nijkjabaabbnjknikaaaaanjkibbnkaakijkkkkkkikkikik

18、kjkkkkikkijkijkikikijkij1,1011,) 1,.(2 , 1) 1()()()()() 1()()()()() 1()() 1()() 1(3.13) )( kiiaiia)(kiia)( kiia（i）系数矩阵中对角线上的元素）系数矩阵中对角线上的元素 都不应都不应 为零，因为在消元的过程中，不断用为零，因为在消元的过程中，不断用 作为除数，倘若有一个作为除数，倘若有一个 为零，为零， 计算就无法进行下去。计算就无法进行下去。（ii）在系数矩阵）在系数矩阵A每一行的元素中，每一行的元素中， 的绝对值最好比同一行的其他元素都的绝对值最好比同一行的其他元素都 大，这大，这

19、 样在作除法运算时，引起的舍样在作除法运算时，引起的舍 入误差就比较小。入误差就比较小。0Adet用用Gauss消去法解线性代数方程组时，为了能求到消去法解线性代数方程组时，为了能求到最后的解（尽管已经具备了最后的解（尽管已经具备了 的条件），并的条件），并使解尽可能的精确，应使解尽可能的精确，应注意注意如下两点：如下两点：n对照上节中差分格式，不难看到，由有限差分法对照上节中差分格式，不难看到，由有限差分法得到的系数矩阵是能够满足以上两个条件的。因得到的系数矩阵是能够满足以上两个条件的。因为每个节点方程（第为每个节点方程（第i个方程）都代表着一个单个方程）都代表着一个单元（第元（第i个单元）

20、与周围单元或外界环境的热量个单元）与周围单元或外界环境的热量交换关系。在这些热量交换中，无疑都与该单元交换关系。在这些热量交换中，无疑都与该单元的温度的温度(Ti)有关，有关，Ti 的系数的系数aii 当然不能为零。当然不能为零。 n而且在第而且在第i个方程中，个方程中，Ti 的地位比其周围单的地位比其周围单元温度更为突出，表现在系数上，它的绝元温度更为突出，表现在系数上，它的绝对值总是最大的。对于给定温度的节点方对值总是最大的。对于给定温度的节点方程而言，这种性质更明显。因为方程中除程而言，这种性质更明显。因为方程中除了该节点温度以外再也没有别的节点温度，了该节点温度以外再也没有别的节点温度

21、，它的系数当然也就最大了。它的系数当然也就最大了。n综上所述，由于对稳定导热问题用有限综上所述，由于对稳定导热问题用有限差分法得到的代数方程具有上述性质，差分法得到的代数方程具有上述性质，因此在求解方程组时，可大胆放心使用因此在求解方程组时，可大胆放心使用Gauss消去法。（对于更一般的方程组，消去法。（对于更一般的方程组，目前更多采用主元素消去法，而不用目前更多采用主元素消去法，而不用Gauss消去法。）消去法。）4 迭代法迭代法n前面介绍的解线性代数方程组的直接法对于前面介绍的解线性代数方程组的直接法对于阶数不是很高的问题是非常有效的，这种场阶数不是很高的问题是非常有效的，这种场合一般不使

22、用下面介绍的迭代法。合一般不使用下面介绍的迭代法。n 对于高阶稀疏矩阵，尽管提出了很多特殊的直对于高阶稀疏矩阵，尽管提出了很多特殊的直接法来处理它们，在运算量和存储量的节省接法来处理它们，在运算量和存储量的节省 方面也取得了很大的进展，但仍然难于克服存方面也取得了很大的进展，但仍然难于克服存 储需要量大的缺点，特别在不具备大型计算机储需要量大的缺点，特别在不具备大型计算机 的条件下，采用下面介绍的的条件下，采用下面介绍的迭代法迭代法更为合适。更为合适。 n迭代法的优点：迭代法的优点：由于不需要存储系数矩由于不需要存储系数矩阵的零元素，所以占用的存储单元少。阵的零元素，所以占用的存储单元少。同时

23、程序也比较简单，对于稳定导热用同时程序也比较简单，对于稳定导热用有限差分法所得到的方程组，求解收敛有限差分法所得到的方程组，求解收敛较快，因此广泛地被采用。较快，因此广泛地被采用。n 迭代法的缺点是：迭代法的缺点是：它所得到的是一种近似解，它所得到的是一种近似解， 在运算过程中需要进行多次迭代才能达到收敛在运算过程中需要进行多次迭代才能达到收敛 指标的要求，而迭代次数事先是不知道的，指标的要求，而迭代次数事先是不知道的， 这样，往往要耗费较多的时间。因此，这样，往往要耗费较多的时间。因此， 一般地讲，直接法与迭代法各有优缺点。一般地讲，直接法与迭代法各有优缺点。 n迭代法所要讨论的问题，仍是如

24、下线性代数方程组：迭代法所要讨论的问题，仍是如下线性代数方程组：11 112 21121 122 222.1 12 2a Ta Ta Tbn na TaTaTbn na TaTaTbnnnn nn(3.14) 简写成：简写成：nibTaijijnj, 2, 11(3.15) n迭代法的迭代法的基本思想基本思想: 构造一个由构造一个由 组成的向量序列，使其收敛于某个极限向量组成的向量序列，使其收敛于某个极限向量 ，并且，并且 就是方程组就是方程组(3.14)的精确解。的精确解。n21TTT,*2*1,nTTT*2*1,nTTT根据构造向量序列的方法不同，常用的有简单根据构造向量序列的方法不同，常

25、用的有简单迭代法，迭代法，GaussSeidel迭代法迭代法与与超松弛迭代法超松弛迭代法，下面分别予以介绍。下面分别予以介绍。n1简单迭代法简单迭代法 最简单的迭代法称为简单迭代法，也称同步迭代最简单的迭代法称为简单迭代法，也称同步迭代法，法，Jakobi迭代法。迭代的最终目的是求解方程迭代法。迭代的最终目的是求解方程组组(3.14)中的中的 。 n21TTT,niaii, 2, 1, 0前面已经说到，用有限差分法（包括有限元法）前面已经说到，用有限差分法（包括有限元法）得到的代数方程组能保证系数矩阵得到的代数方程组能保证系数矩阵A对角线上对角线上元素不为零，即元素不为零，即 n则可将式则可将

26、式(3.14)改写成：改写成： )()()(11111232312121222121211111nnnnnnnnnnnTaTabaTTaTaTabaTTaTabaT(3.16) 其中任一方程均可写成：其中任一方程均可写成：niTabaTjijnijjiiii, 2, 1)(11(3.17) (1). 任意给定各节点上的温度值任意给定各节点上的温度值 作为解的第零次近似，把它们代入式作为解的第零次近似，把它们代入式(3.17)的的 右端，由此算得右端，由此算得), 2, 1()0(niTiniTabaTjijnijjiiii, 2, 1)()0(11) 1 (计算步骤：计算步骤：(2). 作为解

27、的第一次近似，把第一次近似得到的解作为解的第一次近似，把第一次近似得到的解 再代入式再代入式(3.17)的右端，得到解的第二次近似。的右端，得到解的第二次近似。 一般地讲，在已得到解的第一般地讲，在已得到解的第k次近似次近似 后，代后，代 入式入式(3.17)右端，得右端，得)(kiTniTabaTkjijnijjiiiki, 2, 1)()(11) 1( 为解的第为解的第k+1次的似。这样得到的序列次的似。这样得到的序列 为方程组的近为方程组的近似解。只要方程组似解。只要方程组(2.43)存在唯一解，则不存在唯一解，则不论零次近似如何选取，当论零次近似如何选取，当 时，此序时，此序列列 必然

28、收敛，且收敛于方程组的必然收敛，且收敛于方程组的解解 。实际计算中。实际计算中k不可能取不可能取 ，但可以说，当但可以说，当k充分大时，序列充分大时，序列 已足够精确地接近方程组的。已足够精确地接近方程组的。 ),(,)()()(210kTTTknk2k1kn21TTT,*2*1,nTTTn21TTT,n通常，对充分大的通常，对充分大的k，其相邻两次迭代解，其相邻两次迭代解 之间的偏差小于预先给定的适当小量之间的偏差小于预先给定的适当小量 ，即满足即满足),(,)()(n21kTTk21k1niTTkiki, 2, 1|)() 1(就结束迭代过程，而取就结束迭代过程，而取 作为方程组作为方程组

29、(3.17)的近似解。的近似解。 ), 2, 1()(niTkin2Gauss Seidel 迭代法迭代法 简单迭代法虽然计算程序简单，但它计算每一个简单迭代法虽然计算程序简单，但它计算每一个 都要用到全部都要用到全部 的值，因此在计的值，因此在计算机上，必须有两套工作单元来存放全部节点的算机上，必须有两套工作单元来存放全部节点的旧值和新值。为了节省工作单元，并加快迭代收旧值和新值。为了节省工作单元，并加快迭代收敛速度，对上述计算作如下修改。敛速度，对上述计算作如下修改。), 2, 1() 1(niTki)(kiTn方程组方程组(3.16)中的第一式仍写成中的第一式仍写成)()(1)(2121

30、111) 1(1knnkkTaTabaT或或)()(121111) 1(1kjjnjkTabaT方程组方程组(3.16)中第二式，其中中第二式，其中 换成第一次换成第一次算得算得 ，)(1kT) 1(1kT)()(23) 1(1212122) 1(2kjjnjkkTaTabaT(3.18) 即即 n按此规律，可得到的一般关系式：按此规律，可得到的一般关系式： )()(1) 1(111) 1(kjijnijkjijijiiikiTaTabaT(3.19) 这种计算方法称为这种计算方法称为Gauss-Seidel 迭代法，也迭代法，也称为异步迭代法。用称为异步迭代法。用Gauss-Seidel迭代

31、法时，迭代法时，必须将节点或单元按顺序排列，并按顺序逐个必须将节点或单元按顺序排列，并按顺序逐个进行迭代。进行迭代。n由式由式(3.19)即可看到，用即可看到，用Gauss-Seidel迭代法进行迭代求解时，只需用一套工作迭代法进行迭代求解时，只需用一套工作单元存放近似值单元存放近似值 或或 ，这样节省，这样节省了工作单元。同时在迭代过程中，有一半了工作单元。同时在迭代过程中，有一半用迭代的新值，加速了收敛速度。因此，用迭代的新值，加速了收敛速度。因此，Gauss-Seidel迭代法是一种常用的方法。迭代法是一种常用的方法。)(kiT) 1( kiTn3超松弛迭代法超松弛迭代法 实际计算表明，当节点个数