线性代数方程组的直接解法赖志柱

1、第二章线性代数方程组的直接解法教学目标：1. 了解线性代数方程组的结构、基本理论以及相关解法的发展历程；2. 掌握高斯消去法的原理和计算步骤，理解顺序消去法能够实现的条件，并 在此基础上理解矩阵的三角分解（即 LU分解），能应用高斯消去法熟练计算简 单的线性代数方程组；3. 在理解高斯消去法的缺点的基础上，掌握有换行步骤的高斯消去法，从而 理解和掌握选主元素的高斯消去法，尤其是列主元素消去法的理论和计算步骤， 并能灵活的应用于实际中。教学重点：1. 高斯消去法的原理和计算步骤；2顺序消去法能够实现的条件；3. 矩阵的三角分解（即LU分解）；4. 列主元素消去法的理论和计算步骤。教学难点：1.

2、高斯消去法的原理和计算步骤；2. 矩阵的三角分解（即LU分解）；3. 列主元素消去法的理论和计算步骤。教学方法：教具：在自然科学和工程技术中，许多问题的解决常常归结为线性方程组的求解， 有的问题的数学模型中虽不直接表现为线性方程组，但它的数值解法中将问题 “离散化”或 “线性化”为线性方程组。例如，电学中的网络问题、船体数学 放样中建立三次样条函数问题、最小二乘法用于求解实验数据的曲线拟合问题、 求解非线性方程组问题、用差分法或有限元法求解常微分方程边值问题及偏微分 方程的定解问题，都要导致求解一个或若干个线性方程组的问题。目前，计算机上解线性方程组的数值方法尽管很多， 但归纳起来，大致可以

3、分为两大类：一类是直接法（也称精确解法）；另一类是迭代法。例如线性代数 中的Cramer法则就是一种直接法，但其对高阶方程组计算量太大，不是一种实 用的算法。实用的直接法中具有代表性的算法是高斯（ Gauss）消元法，其它算 法都是它的变形和应用。在数值计算历史上，直接法和迭代法交替生辉。一种解法的兴旺与计算机的 硬件环境和问题规模是密切相关的。 一般说来，对同等规模的线性方程组，直接 法对计算机的要求高于迭代法。对于中、低阶（ n : 200 ）以及高阶带形的线性 方程组，由于直接法的准确性和可靠性高， 一般都用直接法求解。对于一般高阶 方程组，特别是系数矩阵为大型稀疏矩阵的线性方程组用迭代

4、法有效。§ 2.1基本定理和问题设具有n个未知数的n个方程的线性方程组为行11为 +ai2X2 +川 +ainXn =bi(2.1)a2lXi +a22X2 +HI +a2nXn =b2min、an1 Xi * an2X2 刊 11 ann X bn或nZ ajXj =b( i =12|il,n)(2.1)j 4其矩阵形式可以表示为Ax = b其中a12IIIa1nA =a21a22IIIa2nIIIHIHIIH&n1an21*1ann丿X =(Xi，X2，|,Xn)T(2.2)b=(b1,b2,|,0)T其增广矩阵为a11a12IHa1na1,n 卅a21a22IHa2na

5、2,n +IHIHIII川IH(an1an2IHannan,n1 )A=A|b-或(ai,n * = b， i = 1,2, I ： I , n )(2.3)则方程组(2.1 )、(2.1)'(2.2 )和(2.3 )是一一对应的。若用r(A)表示矩阵A的秩，则有关于线性方程组(2.2 )的解的存在性的基 本定理(在高等代数或线性代数中有证明)：定理2.1( 1)方程组(2.2 )有解的充分必要条件是r(A)二r(A)；(2) 若r(A)二r(A)二k 5，则方程组(2.2 )具有一族解，其解可表示为n -k(i)X 二 X _CiXi £其中x”为Ax =b的任一特解，X(

6、i)是齐次方程组Ax = 0的解，且x,x，|",*7线性无关，c为任意常数。（3）若r（A）=r（A）二n，则方程组（2.2 ）有唯一解。定义2.1如果两个方程组的解相同，则称这两个方程组为 同解（等价）方 程组。不难证明：将方程组中任意两个方程交换次序， 所得方程组和原方程组为同 解方程组；将方程组中任一方程的一个倍数加到另一个方程上，所得方程组和原 方程组为同解方程组。用（Ei）表示（2.1 ）的第i个方程或（2.3 ）的第i行，记（2.3 ）的第i行与 第j行互换为（E ）（Ej ），而（2.3 ）的第j行乘以a HO加到第i行记为 （Ej E）（EJ。这是矩阵的初等变换，相

7、当于A|b左乘一个初等矩阵。同样 的运算符号，我们也理解为方程组（2.1 ）作相应的变换。经过这些变换后得到 的方程组与方程组（2.1 ）同解（或等价）。对于线性方程组（2.2 ）的求解，在理论上并不存在困难。若r（A）二n，即A 为非奇异（可逆）矩阵，它的行列式 D=detA = O，贝U应用Cramer法则可求得DXi = D （ i =1,2,川，n ）其中Di是用b代替A中第i列而得到的相应的行列式。然而在实际中，当未知数 的个数n比较大时，按Cramer法则进行计算，其工作量就会大得惊人，因而该 方法在实际操作中并不可行。n阶行列式共有n!项，每项都有n个因子，所以计 算一个n阶行列

8、式需要做（n-1） n!次乘法，我们共需要计算n,个行列式，要计算出X,还要做n次除法，因此用Cramer法则求解线性方程组（2.2 ）就要做2N = （n 1） （n -1） n! n 二（n -1） n! n次乘除法（不计加减法）。如nh0时，N =359251210 ；当n =20时，N : 9.7073 1020，可见，在实际计算中Cramer法则几乎没有什么用处。本章的 主要目的就是研究求解线性方程组（2.2 ）的有效算法。某一算法的效率可以用下列两个主要的准则来判断：（1）该算法的计算速度 如何？即计算中要设计多少次运算？（2）计算所得到解的精度如何？这两个准则是针对在计算机上求解

9、高阶方程组而提出的。由于线性方程组阶数很高时求解所需要的计算量极其巨大，因而很自然地提出准则（1）。准则（2） 的提出，是由于在实际问题当中舍入误差的影响，可能使计算解产生偏离真实解的不可忽视的误差。特别地，在解高阶方程组时涉及大量的运算，舍入误差潜在地积累有可能造成计算解对真实解的严重偏离。后续章节还要详细地研究误差的影响。在研究（2.2 ）的数值方法之前，先考察一下求解中会遇到的一些困难，这有助于理解后面将要提出的一些数值计算方法§ 2.2 Gauss消去法从这一节开始，我们来讨论线性方程组（2.1 ）或（2.2）的直接解法。所谓 直接法，就是它只包含有限次的四则运算，若假定每一

10、步运算过程中都不产生舍 入误差，计算的结果就是方程组的精确解。这种方法中最基本和最简单的就是 Gauss消去法及其变形。2.2.1 Gauss消去法的计算过程设方程组（2.1 ）或（2.2 ）的系数矩阵A非奇异，并记a(1) =aj，i =1,2,1山 n， j =1,2川1, n,n 1A|b=A|b这样方程组（2.2 ）又记为Ax =b。a（1）要完成Gauss消去法的第1步须先假定 事=0，令1和二丢（i =2,Hl,n ），冃1则运算（-li1E1 Ei） > （Ei）（ i=2川l,n ）将A|b变换为A | b(2)二ajIIIaa(1)a22)IIIa22)f2)a2,n

11、卅IIIHIIIIIIIan2)IIIan?6,n4fa：；(2.4 )其中a（1）=2川 |,n 1a（2）=a（1）Ti1a1（1），h =哈，i = 2,|H,n，an对应（2.4 ）的方程组是Ax=b，它与（2.2 ）等价，而其第2至第n个方程 中的X1项已经消去。般地，设消去法已进行了 k-1步，得到方程组A（k）x二b（k）。此时对应的增广矩阵为(1)a)1(1) a12 a2?IllIIIm m(1) a1n d2) a2n口 a爲A(k) |b(k)=III(k)hiIII(k)IIIk)akk)HIakn)ak,nH1IIIIIIHIIIIank)1*1a(k)6na(k)

12、an,n4l /(2.5 )a(k)假设 akk)=O ,令 h =召(i=k 1川,n ),则运算(_hEk EJ > (EJ akk(i=k -1|,n )的结果是方程组A(k1)x=b(k1)，对应增广矩阵为A(5|b(5, 其中的元素为a(k),i=1,il|,k; j=1li, n+1a(k ° Wajk) -hkakjk), i =k 1,n; j k 1, n1(2.6)、0, i=k+1,川，nlik = aik)/ akk), i = k 1J | |, n如果依次有akk"( k=1川I, n-1)，则可进行第门-1步，得到与(2.2 )等 价的方

13、程组A(n)x =b，其中A(n)是一个上三角阵，且'a口ana12IIIHIa1n)d1,n 卅A(n)|b(n) =a2?IIIIIIa22)a2,n 卅IIIIIIIII卅_(n)(n)annan,n十丿这样就完成了消去过程。因为A非奇异，故有an；式0。接下来解A(n)x = b(n)，因为A(n是上三角阵，这只要用逐次向后代入的方法即可，这个过程称为回代过程， 其计算公式为a(n)an,n 1Xn = a(n)(2.7)annna(in+_ 送 a(i)XjX =, i =n-1, n-2,川,1aii以上有消去过程和回代过程合起来求出(2.2 )的解的过程就称为Gauss消

14、 去法，或称为顺序Gauss消去法。从(2.6 )可以看出，消去过程的第k步共含有除法运算n-k次，乘法和减法运算各(n -k)(n / k)次，所以消去过程共含有乘除法次数为5nn丄nj(n-k) '、(n -k)(n 1 -k)二k 4k 4含加减法次数为n 1' (n -k)(n 1 -k)k 4而回代过程含乘除法次数为叮，加减法次数为咛)，所以Gauss消去法3总的乘除法次数为- n23332 厂3nnnn5nn，加减法次数为333263当n =10时，用Cramer法则需要3592512103.6 108乘除法，而用 Gauss消去法仅需430次乘除法运算例2.1用G

15、auss消去法解方程组:(1 2(1)12-334 'z15 '-183-1_1X2-151 1111X36<31-1104丿<2<2 432 X21人X3丿(2)解：(1)对增广矩阵进行初等变换f12 3 14、<12314 0 12 8(去1 七3)-3=0 1 2 8? 4113 y0-515丿x-i 2x2 3x3 =14得等价的方程组x2 2x3 8，解得x3=3， %=2，= 1I 5x3=15(2)对增广矩阵进行初等变换12-334、12-334)371537150505221572215(E3 托4E41129211、1129200003

16、6113611735791000000X36)33丿(6e2 E3JE32.12-3 3 415 '-18 3-1 -1-15111161-112z12 -33415、c37150 一 -52 2205321944347770 0<444丿3(空曰 E2) jE21(爲E1引乓12(pE1 E4)E4得等价方程组<12x1-3x?32x2113 x33x3 4x4 = 157 a 152x3 5x4229 仆 x4 = 11691 门x = 033回代得X4=0 , X3=3 , x2 = 2 , X1an设矩阵A的顺序主子式为，即aiiIIIai1IIIIIIIIIali

17、IIIaii，i =1,21,n(2.8 )2.2.2消去法的进一步讨论，矩阵的LU分解从上面的消去过程可以看出，Gauss消去步骤能顺序进行的条件是 £,11扁：爲全不为零。则有下面的定理:定理 2.2 ai(i)(i= 1,2川I,k )全不为零的充分必要条件是A的顺序主子式:i -0( i =1,2川l,k)，其中k - n。步A(m)由A证明：设ai(i-0( i=1,2,Hl,k )，则可以进行消去法的k-1步，每 逐次实行(TjEj Ei) > (Ei)的运算得到，这些运算不改变相应顺序主子式之值,所以有(i) a-11(i)ai2a22IIIIIIIII(i)a1

18、ma2mIH(m) amm这样便有A皿芒0（ m=i,2,IH,k），必要性得证。用归纳法证明充分性。k=i时显然成立。设命题对k-i成立。现设0,1 耳=0,厶k = 0。由归纳假设有a； = 0,1 IhakT=0 , Gauss消去法就可 以进行k -1步，A约化为Ai（ik）A；八< oa22）其中Ak）是对角元为aiT,a2?,ll）,ak：仁的上三角阵。因为A（k）是通过消去法由A逐步得到的，A的k阶顺序主子式等于A（k）的k阶顺序主子式，即A（k）*久=2霁制）川ak粘ak：由也k式0可推出ajk鼻0。定理2.3对方程组Ax二b，其中A非奇异，若A的顺序主子式均不为零， 则

19、可以Gauss消去法求出方程组的解。定义2.2设矩阵A =（a）価每一行对角元素的绝对值都大于同行其它元素n绝对值之和，即 丨' jajjl，i =i,2川|,n，则称A为严格（行）对角占优矩阵jm定理2.4设线性方程组Ax =b的系数矩阵为严格对角占优矩阵，则用顺序Gauss消去法求解时aj（ i =1,2,ltl,k ）全不为零。证明：（略）下面我们来讨论Gauss消去过程用矩阵运算表示的形式。第1步令A二A，作一次运算（-liiErEJ > （Ei） , i=2川l,n，这相当于A左乘矩阵Mi = -li| 1, i=2川|,nk4'、b第1步的全过程相当于LJA|

20、b二A|b，其中L|=MnMn|M2121 1<_l n1L2L2 川 L.A二 A(k)设k -1步后系数矩阵化为A(k)，其分块形式写成其中A(1k)为上三角的k -1阶方阵，A22°为n -k 1阶方阵，设其左上角元素a；：=0 ，则下一步的乘数为 lia(kk) /akk) ， i = k T,|l, n。若记ek =( 0IJ011 T 0是第k个分量为1的单位向量，记l(k(0j|,0,lk+,kJH,lnk)T，其前k个分量为零，从而有eTl(k)=0。第k步中系数矩阵的约化可用矩阵运算描述为LkA(k)A(k 1)fA(k+)A11OA(k+)、A12(k 1)

21、A22其中阳1)是上三角的k阶矩阵，a22 d)是n-k阶方阵，而(2.9)1Tk 1,k1-lnk可以验证(i -l(k)eT)(i l(k)eT"(i l(k)eT)(i -l(k)eT)l(k)(&l(k),I即有 L： =(l 严&)=1 - l(k)eT。这样，经过 n1 步得到 LnJLn|L1A(1 A(n),这里的A(n)是上三角阵，记U =A(n)，又记(LnJll|L)J =(l IeT川 1(1 1(7)I21I3132(2.10)n1n2IIIL是一个对角线元素全为1的下三角阵，这种矩阵称为单位下三角阵。L的对角 线以下元素就是各步消去的乘数。

22、最后我们可以得到 A-LU，其中L是一个单 位下三角阵，U是一个上三角阵。定义2.3将矩阵A分解为一个下三角矩阵 L和一个上三角矩阵U的乘积A二LU，称为对矩阵A的LU分解或三角分解。当L是单位下三角矩阵时称为 Doolittle (杜里克尔)分解。当U是单位上三角矩阵时称为 Crout (克洛特) 分解。由上面的分析过程知，Gauss消去法的实质是将系数矩阵分解为一个下三角 矩阵和一个上三角矩阵相乘，即将系数矩阵进行LU分解。在矩阵A的LU分解A = LU 中,我们将U写成U = DU，其中D是对角矩阵，U是单位上三角阵，进一步记L二LD，它是一个下三角阵，这样有A 二 LU 二 LDU =

23、(LD)U其中匚是一个下三角阵，U是单位上三角阵，此即A的Crout分解。在矩阵A的Doolittle 分解A二LU中，我们将上三角阵U写成DU的形式， 这里的D为对角阵，U为单位上三角阵，这样得到 A = LDU，其中L为单位下 三角阵，D为对角阵，U为单位上三角阵，称其为 A的LDU分解。定理2.5设非奇异矩阵A Rnn，若其顺序主子式冷(i=1,IH, n-1)都不等于 零，则存在唯一的单位下三角阵L和上三角阵U，使A二LU。证明：上面的分析过程已经说明了非奇异矩阵A可作LU分解，下面只需证明分解的唯一性。设A有两个分解式A二屮1和A二L2U2，其中JI都是单位下三角阵，UnU2 都是上

24、三角阵，则有 屮1 =L2U2。因为A非奇异，从而J,L2， U1,U2都可逆。故 在屮产L2U2两边同时左乘L,和右乘U2，这样得到SU'rLrL。因为Uf1仍 为上三角阵，故U1U2'也是上三角阵，同理可得 LL2是单位下三角阵，结合U1U2=L/L2 知只可能 UU2,=L：L2=1，即有Li=L2 ,Ui=U2。证毕可以证明，对于奇异矩阵 A Rnn依然满足定理2.5。而且，从上面的推导过程可以看到，对A作LU分解时，其中的U为还可以验证A的顺序主子式为对应的L和U的顺序主子式的乘积，而A的m阶顺u11 u12 IIIu1nanajIlla(1)a1nu22 HIu2n

25、a2?III2)a2n+I14 1-1IIIIIIunna(n) annU序王子式满足-:m = UiiU22 HlUmm， m=12l|l, n。定理2.6非奇异矩阵A Rnn有唯一的LDU分解的充分必要条件是 A的顺序主子式.(i =1,1山n-1)都不等于零。§ 2.3主元素消去法2.3.1有换行步骤的消去法顺序Gauss消去法可以进行的条件是 a” =0(i =1,2,川,n-1)，如果有某个 a =0，消去法就不能继续使用。例如，若an = 0,贝U第1步就不能进行，但我 们可以在A的第1列中找出一个非零元ai1，进行换行(EJ(巳)后再做消去法, 其他各步可以类似处理。有

26、时候虽然 审工0，但|审|很小，顺序Gauss消去法可以顺利进行下去， 但计算过程舍入误差导致误差增长过大，以致结果不可靠，此时称ai(i)为小主元例2.2用三位十进制浮点运算求解1.00X10* +1.00x2 =1.001.00% 1.00x2 二 2.00解：不难发现这个方程组的准确解应该接近(1.00,1.00/。下面我们用顺序Gauss消去法来求解，则有121=亜=1.00 105，a22)二 a22 T21q2 =1.00一1.00 105a11a23 = a23 T21Q3 = 2.00 T.00 105在三位十进制运算的限制下，得到a2?x2鲁=1.00a；代回第一个方程得花=

27、0。这显然不正确。从上例的计算过程不难发现，计算解产生误差的原因是小主元做除数。因此, 为使这种不稳定现象发生的可能性减至最小，在每一步就要寻找一个远不是零的数作除数，既要寻找ik(ik)使|ai(k)Umax|ai(kk)|，并将第ik行与第k行对换。这 k艺勺样akk)的当前值就远不是零，这时的元素akk)称为k列的主元素，ik行叫做主元素 行，这样就引出了列主元消去法。2.3.2列主元消去法与完全主元消去法列主元消去法也称为按列部分选主元的消去法。其具体过程如下:a；?，即| max | ai1(11“ I 1 :i： ：n I i 1|，其中21。因为a非奇异，故有事式0。若h =1，

28、则消去过进行第1步之前，在 A的第1列中选出绝对值最大的元素程与顺序Gauss消去法一样；若ii 1，则先进行换行(E) (EJ，然后类似顺 序Gauss消去法的运算得到A|b。般地，设进行了 k -1步选主元和消去法的步骤，得到A(k)|b(k)。第k步先选主元a(kk，使| ai*戶max |ai(kk) |，即在A(k)的第k列的对角线及对角线以下的元素中选出绝对值最大值a，显然ik 一 k，且由于A(k)非奇异，从而有 嚟=0。若ik =k，贝U进行顺序Gauss消去法的第k步；若ik k，则对A(k)|b先换行 (Ekh (EJ，然后进行类似Gauss消去法的运算。如上进行n -1步

29、选主元、换行与消去法运算，得到A(n)x = b(n)，其中A(n)是一个上三角阵，再回代进行求解。列主元消去算法：1、消去过程：对k =1,2,|山n-1(1)选主元：找 ik k,川,n，使旧丁 尸max |ai(kk) | ；(2)若 aik；)=0，则停止计算(detA=|A|=0)； (3)若i/k，则换行(E”匕);(4) 消元：对i=k+1,川,n , hk"kk)/akk);对 j=k+1川,n+1 , a(* = a(k)爲；)2、回代过程：(1)则停止计算(det A =| A |=0 );(2)Xnn,n 1 .二 ；ann(3)na(i仃送ajXj-2J|I,

30、1，Xi游a()在列主元消去算法的过程中，不是按列来选主元，而是在壮，即城化爱严，)，然后将完全主元消去法：A(k)右下角的n -k 1阶子阵中选主元a( jA(k)|b(k)的第ik行与第k行、第jk列与第k列交换，再进行消去运算完全主元法比列主元法运算量大得多，可以证明列主元消去法的舍入误差 般已比较小，在实际计算中多用列主元法。例2.3用列主元消去法解方程组Ax =b，计算过程取五位数字，其中"-0.002220.4'A|b =10.7812501.381613.9965.562547.4178 解：-0.002A |b=A|b =1J 3.99620.781255.5

31、6252 0.401.3816 a314 7.4178,3.996= 3.996(曰)尸乍3)>5.56250.781251-0.00227.41781.38160.412厂 1/3.996 =0.25025131 二-0.002/3.996 二-0.00050050(E2 42启)_乍2)_(E3 41E1)畑)A(2) |b(2)二3.99605.5625-0.610772.00284-1.00102.00207.4178-0.474710.40371二 2.0028'3.9965.56250 2.00280-0.6107747.41782.0020 0.40371-1.00

32、10 -0.47471.-0.61077I32 ：2.0028=-0.30496回代得:"3.9965.562547.4178 （E3 丄2E2刍A（3） |b02.00282.00200.40371< 00-0.39047-0.35159 x3 = 0.90043, X2 - -0.69850 , x<i =1.9273。其精确解为x = (1.92730,-0.698496,0.900423$而用顺序Gauss消去法，则可解得x = (1.9300, _0.68695,0.88888)丁这个结果误差比较大，这是因为消去法的第1步中，a(1)按绝对值比其他元素小很多所引

33、起的。§ 2.4直接三角分解法Gauss消去法还有许多变形，有些变形是为了利用特殊技巧减少误差，把 Gauss消去法改写为更紧凑的形式，还有一些变形时根据某类矩阵的特性作一些 修正和简化，这些方法可统称为 直接三角分解法。2.4.1矩阵的三角分解设A的顺序主子式=0(k =1,2J|I,n)，则可建立线性方程组Ax二b的Gauss消去法与矩阵分解的关系，即矩阵 A的LU分解。这个问题前面已经讲的比较详 细了，此处不再赘述。2.4.2 Doolittle 分解法首先假设A的顺序主子式 耳都不为零，则 A可作Doolittle 分解，即A-LU，其中L是单位下三角阵，有h =1，“ j时

34、1厂0 ； U是上三角阵，ij时5 =0。仔细写出为a12IIIa1n 'nS11u12IHu1n 'a21a22HIa2nl211u22IHU2n+ 1-11 + FF 4 +ill+ 1-1* » hqI -1 I-lan1an2IIIann丿n1In2 HI bunn丿-LU （ 2.11）A 二在前面逐步推导L和U的元素公式都要借助于有关的a(k)来表示。现在强调 指出，只要从给定的A通过比较(2.11)式的两边就可能逐步地把L和U构造出 来,而不必利用Gauss消去法的中间结果，这种方法称为Gauss消去法的紧凑格 式。根据矩阵的乘法规则，比较（2.11 ）

35、式的两边可得min（i,j）aj - likUkj , i,j =12|,n（2.12）k=1先算出u的第1行，在（2.12）令i =1，由S ", hj =0（1 :： j乞n）可得Uij 二a“，j =1,2,川,n接着在（2.12）中令j =1，由aM -li1 U11，从而算出L的第1列h = ai1 / u11 = a1 / a11， i = 2,111, n若U的第1至k -1行和L的第1至k -1列已经算出，则由（2.12）可计算出 U的第k行元素k 二Ukj - akj - l kr U rj ， j = k,k 1|,n（2.13 ）r 二接着再计算出L的第k列元素k -1aik 1 ir U