1、第四章 线性代数方程组的数值解法 4.1 概 述 线性代数方程组(System of Linear Algebra Equations)的求解是数值计算方法中的一个重要课题。现代工程技术或科研过程中所遇到的一些实际问题，常常直接或间接地归结为求解一个线性代数方程组。例如有分支水流的流速分布、建筑结构中的设计计算和应力分析、仪器分析中的质谱分析、光谱分析、色谱分析、电力系统中的电网分析等等；另外有些数值计算方法本身也是以线性代数方程组的数值解、有限元法和边界元等等，其中间过程或者最后都会导致求解线性代数方程组。这些线性代数方程组的解法是十分必要的，有着十分重要的意义。4.1 概 述本章介绍求解n
2、阶线性代数方程组的一般形式是：或简写成矩阵形式：其中A A称为方程组(4.1)式的系数矩阵；B B称为方程组的右端列向量；x x称谓方程组解的列向量。它们分别为 11 11221331121 1222233221 12233nnnnnnnnnnna xa xa xa xba xa xa xa xba xa xa xa xb Axb(4.1)(4.2)4.1 概 述1112131121222322123123jnjniiiijinnnnnjnnaaaaaaaaaaaaaaaaaaaaA12inbbbbb12inxxxxx 求解线性代数方程组的数值计算方法很多，大致可分为两大类：消去（元）法和迭代
3、法。消去法是直接从方程组的系数矩阵入手，经过有限步运算求出方程组的精确解（假如没有舍入误差的话）。迭代法则是将求方程组的问题化为构造一组递推计算结构，从一组近似解出发，用这组递推结构逐步算出精度更高的近似解。上述这两类算法各有其优点和缺点，消去法的计算量小，但程序复杂；迭代法计算量大，精度不高，但程序结构简单。本章将主要介绍求解线性代数方程组的简单高斯消去法和三角分解法，迭代法将在以后相关章节中介绍。4.2 高斯消去法 4.2.1高斯消去法的基本步骤我们以三元线性代数方程组为例，叙述简单高斯消去法（以下简称为消去法）的基本步骤，这各方法是理解其它方法的基础，消去法分为消元和回代两个过程。例4-
4、1：消去过程实际上是对增广矩阵作行初等变换。上例可表示为 1232414312954625xxx24143129546252132rr2152rr2414071/2150147/235322rr2414070.515002.554.2 高斯消去法 这是上三角方程组，它极容易求解：由第三个方程得 ，代入第二个方程得 ，再代入第一方程 。此一过程称为回代过程。从上述消元过程可以看出，三元议程组要经过两次消元（因为不用消）才能把增广矩阵化为拟上三角矩阵。对于一般的n元线性代数方程组，经经过(n-1)次消元才能把相应的增广矩阵变换为拟上三角矩阵。n元线性代数方程组的增广矩阵的一般形式如(4.3)式。3
5、2x 22x 11x 3x4.2 高斯消去法 以下是消元步骤（共分n-1部）：第一步(k=1)：从第一列中消去第2n行中x1的系数，即消去a11下方元素；第二步(k=2)：从第一列中消去第3n行中x2的系数，即消去a22下方元素；(1)(1)(1)(1)(1)(1)1112131111(1)(1)(1)(1)(1)(1)2122232221(1)(1)(1)(1)(1)(1)3132333331(1)(1)(1)(1)(1)(1)1231jnnjnnjnniiiijininaaaaaaaaaaaaaaaaaaaaaaaa(1)(1)(1)(1)(1)(1)1231nnnnjnnnnaaaaaa
6、(4.3)4.2 高斯消去法第j步(k=j)：从第j列中消去第j+1n行中xj的系数；即消去ajj下方元素;第n-1步(k=n-1)：从第n-1列中消去第n行中xn-1的系数；即消去an-1n-1下方元素因此，高斯消元步骤的计算通式可写为：()()(1)()()1,2,3,11,2,1,2,1kikkkkkkkijijkjknikknaRajkknaaRa(4.4)4.2 高斯消去法经过n-1步消元后，增广矩阵(4.3)式变为：(1)(1)(1)(1)(1)(1)1112131111(2)(2)(2)(2)(2)22232221(3)(3)(3)(3)333331()()()1()()1000
7、0000000jnnjnnjnniiiijininnnnnnnaaaaaaaaaaaaaaaaaaaa(4.5)4.2 高斯消去法高斯回代过程的通用递推计算公式可写为：()1()()()11()()1,2,3,2,1nnnnnnnnkkknkjjj kkkkkaxaaaxxaknn(4.6)4.2 高斯消去法4.2.2高斯消去法的计算量分析我们从(4-4)式可以看出，高斯消去法消去过程共分n-1步，第k步变换n-k行：对这些行先乘系数，再从n+1-k元素减去第k行相对应列的元素；因此共需乘、除次数为回代过程求xn需1次除法，求xn-1需1次乘法、1次除法，求x1需n-1次乘法、1次除法，因此共
8、需乘除次数两过程共需要乘除次数为 。当n=20时，N=3060.1111()(11)(1)(25)6nkNnk nkn nn 2(0 1)(1 1)(1 1)12(1)2Nnnnn 3212/3/3NNNnnn4.2 高斯消去法4.2.2 高斯消去法的程序框图与通用程序 STARTINPUT N,A(),BCALL eliminationCALL backOUTPUT X(N)END4.2 高斯消去法10 OPEN(8,&file=EliminationMetod.dat,&status=new)20 DIMENSION a(n,n+1),x(n),b(n)30C 输入n,a(),b40 RE
9、AD(8,(1x,i2,)n50 DO i=1,n60 READ(8,(1x,e12.6)b(I)70 END DO80 DO i=1,n90 DO j=1,n100 READ(8,(1x,e12.6)a(i,j)110 END DO 120 END DO 120C 消元过程130 CALL elimination(a,b,n)140C 回代过程150 CALL back(a,x,n)160C 输出结果x170 DO i=1,n180 WRITE(9,(1x,e12.6)x(i)190 END DO 200 END2102202302402504.2 高斯消去法SUBROUTINE elimi
10、nationK=1,n-1i=k+1,nR=a(i,k)/a(k,k)j=k+1,n+1a(i,j)=a(i,j)-R*a(k,j)end jend jend iend kEND SUBROUTINE elimination4.2 高斯消去法10 SUBROUTINE&elimination(a,b,n)20 DIMENSION a(n,n+1),b(n)30C 生成增广矩阵40 DO i=1,n50 a(i,n+1)=b(i)60 END DO70C 消去过程开始80 DO k=1,n-190 DO i=k+1,n100 R=a(i,k)/a(k,k)110 DO j=k+1,n+1120
11、a(i,j)=a(i,j)-R*a(j,k)130 END DO140 END DO150 END DO160 END SUBROUTINE elimination4.2 高斯消去法SUBROUTINE backX(n)=a(n,n+1)/a(n,n)k=n-1,1,-1s=0j=k+1,ns=s+a(k,j)*x(j)end jX(k)=(a(k,n+1)-s)/a(k,k)end kEND SUBROUTINE back4.2 高斯消去法10 SUBROUTINE back(a,x,n)20 DIMENSION a(n,n+1),x(n)30 x(n)=a(n,n+1)/a(n,n)40
12、DO k=n-1,1,-160 s=0.0 60 DO j=k+1,n70 s=s+a(k,j)*x(j)80 END DO90 X(k)=(a(k,n+1)-s)/a(k,k)100 END DO110 END SUBROUTINE back4.2 高斯消去法4.2.3 选主元法高斯消去法消去过程中，第k步求n-k个系数aik/akk用到的除数akk,称为主元。如果akk为接近于零或零，则消元递推计算通式(4-4)将出现被零除的情况，从而使计算机溢出或无法进行下去。例4-2准确到九位小数的解是x1=0.250001875，x2=0.499998749，若在4位计算机上按高斯消去法求解，则51
13、2121021232xxxx510212325212 10rr 555102104 102 10 4.2 高斯消去法 回代得解 ，显然产生严重失真。造成这种结果的原因，就是小主元10-5的出现。用它作除数会产生大乘数，数乘大乘数变大数，大数吃小数产生舍入误差,如()，从而使 也有误差。因此，为了保证获得正确的结果，在消元过程中可在第k步的第k列的元素akk,ak+1,k,，ank中选主元，即在其中找出绝对值最大的apk，然后交换第k行和第p行，继续进行消去过程。这种消去法为列主元消去法列主元消去法。也可在第k步的第k行的元素akk,ak,k+1,，akn中选主元，即在其中找出绝对值最大的akq
14、，交换第k行和第q行;或在kn行、kn列选绝对值最大的元素 apq，交换k、p行和k、p列，然后继续进行消去过程。这两种消去法为行行主元消去法主元消去法和和全主元消去法全主元消去法。本节讨论应用广泛且易于在计算机上实现的列主元素消去法。列主元素消去法的计算步骤可归纳如下：选主元素；把主元素所在的行与第k行互换；消元；回代 5634 10399997100.4000 10.5x 20 x 2x 4.2 高斯消去法其中、两步与简单高斯消去法中的消元与回代过程完全一样，也就是说，列主元消去法是在简单消去法的基础上引入选列主远过程，这个过程并不影响获得正确结果。列主元高斯消去法的程序框图与通用程序 为
15、了简单和说明问题，将上面的和两步单独用一个子程序。因此，列主元高斯消去法的程序框图只需在简单高斯消去法的流程图中增加一个选列主元子程序gaussp，具体如下：4.2 高斯消去法SUBROUTINE gausspk=1,n-1 a(p,j)=ti=k+1,nc=abs(a(i,k);p=iend it=a(k,j);a(k,j)=a(p,j)end jEND SUBROUTINE gausspc=abs(a(k,k);p=kIf abs(a(i,K)c?YNj=1,n+1 end k4.2 高斯消去法10 SUBROUTINE gaussp(a,b,n)20 DIMENSION a(n,n+1)
16、,b(n)30C 生成增广矩阵40 DO i=1,n50 a(i,n+1)=b(i)60 END DO70C 选主元开始80 DO k=1,n-190 c=abs(a(k,k);p=k100 DO i=k+1,n110 IF(abs(a(i,k)c)THEN120 c=a(i,k);p=i130 END IF140C p,k行交换开始150 DO j=1,n+1160 t=a(k,j);a(k,j)=a(p,j)170 a(p,j)=t180 END DO 190 END DO200 END SUBROUTINE gaussp4.2 高斯消去法4.2.4 高斯-约当法前面所讲的是如何将一个线性
17、方程组的系数矩阵改变为一个上三角矩阵，然后回代求解。本节所讲的是如何将系数矩阵改变为一个对角矩阵，即，除对角线元素外，其他元素均为零(1)(1)1112131111111(2)(2)2122232222221()()1231()()12310000000000000000jnnjnniiiiiijiniijinnnnnnnjnnnnnnnaaaaabaaaaaaabaaaaaaabaaaaaaabaa 4.2 高斯消去法列主元高斯-约旦消元过程的递推计算通式结构为：高斯-约旦消去法消元过程的计算量比高斯消去法有所增加（约为1.3倍），但是总的计算量增加并不大，另一方面，这一方法的程序结构简单，
18、所以在实际中得到了广泛的应用。高斯-约旦消去法的程序框图与通用程序自行设计。()()(1)()()1,2,3,11,2,(),1,1kikkkkkkkijijkjknikknikaRajk knaaRa()()(1)()()1,2,3,11,2,()1,2,1kikkkkkkkijijkjknikknikaRajkknaaRa改进后(4.8)(4.7)4.3 三角分解法 线性方程组的另一直接解法是三角形分解法，即方程组（4-2）的系数矩阵A分解成两个形式简单的三角形矩阵L和U的乘积：A=LU。从而求解Ax=b的问题转化为三角形方程组 Lx=y Uy=b 其中(4.9)(4.10)2112100
19、101nnlllL11121222000nnnnuuuuuuU4.3 三角分解法比较等式左边和右边乘积矩阵LU的第r行主角元右边（含主角元）的对应元素，得再比较等式两边第r列主角元以下（不含主角元）的对应元素，得当r=1时，有1(,;1,2,)rrirkkikal uirn rn(4.11)1(1,;1,2,1)ririkkrkal uirn rn(4.12)11(1,)iiuain 1111(2,)iialinu(4.13)(4.14)4.3 三角分解法假定已求出U的第1至第n-1行和 L的第1至第n-1列，由（4-11）式计算出U的第r行元素和由（4-12）式计算出L的第r列元素（4-13
20、）式至（4-16）式所表示的矩阵分解称Doolittle分解。类似地，若U为单位上三角矩阵，而L为下三角矩阵，则有11(,;2,)rririrkkikual uirn rn(4.15)11(1,;2,1)rriikkrkirrral ulirn rnu(4.16)11(,;1,2,)riririkkrklal uirn rn(4.17)4.3 三角分解法规定 。称上述分解为Crout分解。实现了系数矩阵A的Doolittle或Crout分解后，方程组Ax=b可以分解成为Ly=b和Ux=b这两个三角矩阵来进行求解。设U为上三角矩阵，而L为单位下三角矩阵，并且 ，根据公式（4-4）和（4-6）得与
21、Doolittle分解相对应的方程组解为11(1,;1,2,1)rrirkkikrirral uuirn rnl010ikkrkl u(4.18)0(1,)iiuin1111(2,)rrrriiiybybl yrn(4.19)4.3 三角分解法类似地，与Crout分解相应的求解公式1(1,1)nnnnnrriii rrrryxuyu xxrnu(4.20)111111(2,)rrriiirrrbylbl yyrnl(4.21)4.3 三角分解法无论Doolittle分解还是Crout分解，其运算量均小于高斯消去法。为保证它们能顺利、稳定进行，也可选列主元。用Crout分解法求解方程组Ax=b的
22、步骤如下：用（4-17）式求lir;确定最大列主元；用（4-18）式求uir,进行分解；回代。通用公式为1(1,1)nnnrrriii rxyxyu xrn(4.22)1111,1,1,111,2,1,1,()1,1max(,iijijikkjkijijijkkiknpijiij ni ni nikk nk iinjinjinaaa aiaaa aiinaaaaa apip i Rxy分解,算 行）算 列（选主元）解回代）当时交换行(4.23)4.3 追赶法求解三对角占优方程组若作Crout分解，则有111222233331111nnnnnnnbcfabcfabcfabcfabfA f 4.3
23、 追赶法比较Crout解，易知回代得 。1223311nnnnlalalalalL1122331111111nnnuyuyuyuyyU 11111,/lbyfl11111,(2,)iiiiiiiiiiiiculbaulfa yyinl1,(1,1)nniiiixyxyu xin(4.24)(4.25)4.3 追赶法按照以上公式求解Ax=b的方法就叫做追赶法，其中ui、li 和yi 称追，回代称赶，乘除法运算量仅为5n-4，远比一般方程组的高斯消去法或三角分解法的运算量少。不用选主元，就可保证顺利、稳定进行计算。4.4 舍入误差对解的影响 本节研究线性代数方程组解的误差。衡量数的误差用绝对值。方
24、程组的解是向量，衡量其误差，自然也会想到绝对值的概念。4.4.1 向量与矩阵范数 众所周知，数x的绝对值x是x的函数:x=(x),具有以下三性质：（1）（2）对任意实数（3）推广到向量 ，具有如下类似性质的函数：。该函数成为向量 x的范数或模或长度，它也应该具有以下性质：（1）非负性（2）齐次性：对任意实数（3）三角不等式：00,00 xx时,xx xyxyT12(,)nx xxx12(,)nx xxx0,0 xx时,xxx+yxy4.4 舍入误差对解的影响几何上三角不等式表示：在以向量x、y及其和x-y构成的三角形中，一边之长不超过另两边边长之和，即从而表示以x、y和x-y为边的三角形中，一
25、边之长不小于另两边边长之差。常用向量范数有3种，即：分别称为1范数、2范数和无穷大范数。不难验证它们具有性质（1）和（2）。对于性xxyyxyyxyxy121nxxxx222 1/2122()nxxxx1maxii nx x4.4 舍入误差对解的影响性质（3），这里已2范数为例，给出证明，以供参考。最后不等式称Cauchy不等式。注意对任意实数 2222222222111221222111()2()()()nnniiiiiiiniiinnniiiiiiixyxyx yxyx y xyx+yxyxy22221111()2()0nnnniiiiiiiiiixyxx yy4.4 舍入误差对解的影响由
26、 二次式的判别式 ，便可推出Cauchy不等式，因而2范数的性质（3）得证。对于n阶方阵，具有如下4种性质的函数 称为矩阵范数。（1）（2）对任意实数（3）（4）容易验证矩阵函数具有上述4种性质，因而，它是一种范数，称为矩阵的F范数。常用矩阵范数是如下三种范数：0()AA0,00OAA时,AAA+BA+BABA B21/211()nnijFijaA4.4 舍入误差对解的影响其中 表示矩阵 的谱半径，即 特征值绝对值的最大值。三种范数分别称为1范数或列范数，无穷大范数或行范数，2范数或谱范数。可以证明，矩阵的这三种范数还具有以下3种性质：（满足性质（1）的矩阵范数成为相容范数）（1）（2）单位矩
27、阵I的范数（3）时 可逆且 111max()nijj nia A11max()niji nja A2()TAA ATA ATA ATA AAxA x1I1BI-B1()1/(1)I-BB(4.26)4.4 舍入误差对解的影响例4-3 设则1201210111max(1 10),(22 1),(0 1 1)5 Amax(120),(12 1),(0 1 1)4 A110120209221121091011011112TA A4.4 舍入误差对解的影响它的三个根为4.4.1 误差分析 设线性方程组Ax=bAx=b中，A A为非奇异矩阵，x x为方程组的精确解。假定b b有误差b b,则解为x+x+
28、x x,即32209det()0911121338250TIA A2max9.1423.0237iA1239.142,2.921,0.96312222222 1/2(12(1)2(1)11)133.6056F A()A x+xb+b4.4 舍入误差对解的影响 即两边取范数又因由（4-30）和（4-31）式，得另外，若矩阵A A有误差A A，则解为x+x+x x，即AxbAxb1xA b11xA bAbA xAxb(4.27)(4.28)(4.29)(4.30)(4.31)1xbA Axb(4.32)()()AA x+xb(4.33)4.4 舍入误差对解的影响即有虽然A A为非奇异，但不能保证A
29、+A+A A也为非奇异，为保证A+A+A A为非奇异，需对A A作些限制。因为由（4-26）式，若 ，则 存在，因此（4.35)式变为由（4-34）和（4-35）式，得两边取范数()AA xAx(4.34)1()AAA IA A(4.35)11A A11()IA A1111()()AAIA AA(4.36)111()xIA AA Ax(4.37)1111111()11xIA AAA xAA xAA xA AAA(4.38)4.4 舍入误差对解的影响设 ，则得从（4-32）和（4-39）式可以看出，如果线性方程组Ax=bAx=b中，A或的微小变化，都将引起方程解x为 倍，若这个倍数巨大，则方程组
30、定义为“病态方程组”，矩阵A A为“病态矩阵”，否则定义为“病态方程组”或矩阵A A为“良态矩阵”。此时，我们称 为条件数，即常用的条件数有11AA111AA AxAAxA AA(4.39)1A A1A A1cond()AA A(4.40)1cond()AAA1max222mincond()TTA AAAAA A(4.41)4.4 舍入误差对解的影响例4-4 试计算 ，并取6位有效数字。解（1）直接计算若使用Gauss列主元消去法，经消元则可得77121 1010112xxcond()A17771 101 1011011111-1A77177(101)(101)cond()10101AAA12
31、0,1xx4.4 舍入误差对解的影响（2）间接计算，在上线性方程组两边同乘以矩阵则有故用Gauss主元消去法得710001D7121101211xx 171771110111 101 1011A172 2cond()41 10AAA121,1xx习 题4-1 用Gauss消去法求解下列方程组Ax=b4-2 求cond(A)2和cond(A).1111(1)122,02111 Ab26(1)26.00001A上机实习题4-3 用选列主元的高斯消去法求解下列方程组：123412341234123441973052068401366.871.347.45211.588.676.410.880225.71.455.96.1336.56.6xxxxxxxxxxxxxxxx
