现代数值分析课件第2章线性方程组的直接法

1、第二章第二章 解线性方程组的直接方法解线性方程组的直接方法 高斯高斯(Gauss)(Gauss)消元法消元法 矩阵的三角分解法矩阵的三角分解法 矩阵的条件数与方程组的性态矩阵的条件数与方程组的性态11112211211222221122 nnnnnnnnnnna xa xa xba xa xa xba xa xa xb阶线性方程组： , X ,11Tn nijAXbTA nn , , ) , , )(x(bxbab矩阵表示记为这里解线性方程组的两类方法解线性方程组的两类方法: :直接法直接法: : 经过有限次运算后可求得方程组精确解的方经过有限次运算后可求得方程组精确解的方法法( (不计舍入误

2、差不计舍入误差!)!)迭代法：从解的某个近似值出发，通过构造一个无穷迭代法：从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法序列去逼近精确解的方法( (一般有限步内得不到精确解一般有限步内得不到精确解) )一、一、高斯消去法高斯消去法思思路路首先将方程组首先将方程组Ax=b 化为上三角方程组化为上三角方程组，此过程称为此过程称为消去过程消去过程，再求解上三角方程，再求解上三角方程组，此过程称为组，此过程称为回代过程回代过程. .2.1 高斯消去法和选主元高斯消去法高斯消去法和选主元高斯消去法将增广矩阵将增广矩阵的的第第 i 行行 + li1 第第1 1行行，得到：，得到：,)(记)

3、1()1(nnijaAA(1 )(1 )(1 )()1Tbbbbn消去过程：消去过程：第一步第一步：设设 ，计算因子，计算因子0)1(11a)1(11)1(11aalii)1(1)1(1)1(12)1(11.baaan) 2(A) 2 (b其中其中) 1(11) 1()2() 1(11) 1()2(, 3 , 2, ,blbbnjialaaiiijiijij0)( kkka第第k步：步：设设 ，计算因子，计算因子且计算且计算( )( )/(1,., )kkikikkklaaikn (1)()()(1)()()( ,1, .,)kkkijijikkjkkkiiikkaal abbl bi jkn

4、共进行共进行 n 1步，得到步，得到 )()2(2)1(121)()2(2)2(22)1(1)1(12)1(11.nnnnnnnnbbbxxxaaaaaa)()(/nnnnnnabx ) 1., 1()(1)()( niaxabxiiinijjiijiii定理定理2.1：若若A的所有的所有顺序主子式顺序主子式 均不为均不为0，则高斯，则高斯消去法能顺序进行消元，得到唯一解。消去法能顺序进行消元，得到唯一解。iiiiiaaaaA.)det(1111 回代过程：回代过程：例例2.1.1 高斯消元法高斯消元法 32303452536432321321321xxxxxxxxx 323034525364

5、32A 202230445 . 006432A 4200445 . 006432 42445 . 06432332321xxxxxxx1 = -13x2 = 8x3 = 2二、二、选主元消去法选主元消去法为避免这种情况的发生为避免这种情况的发生, 可通过交换方程的次序，选可通过交换方程的次序，选取绝对值大的元素作主元取绝对值大的元素作主元. 基于这种思想导出了主元基于这种思想导出了主元素法素法在高斯消去法消去过程中可能出现在高斯消去法消去过程中可能出现 的情况，的情况，这时高斯消去法将无法进行；即使主元素这时高斯消去法将无法进行；即使主元素 但很小，其作除数但很小，其作除数 ，也会导致其它元素

6、数量级的严，也会导致其它元素数量级的严重增长和舍入误差的扩散重增长和舍入误差的扩散( )0kkka( )0kkkav 列主元消去法列主元消去法在第在第k步消元前，在系数矩阵第步消元前，在系数矩阵第k列的对角线以下列的对角线以下的元素中找出绝对值最大的元。的元素中找出绝对值最大的元。|m ax |0pkkkikkinaa 若若pk, 交换第交换第k个与第个与第p个方程后个方程后, 再继续消去计算再继续消去计算. . 这种方法称为这种方法称为列主元列主元GaussGauss消去法。消去法。 列主元列主元GaussGauss消去法保证了消去法保证了lik1 1 ( ( i = k +1, k +2,

7、，n ).).例例2.1.2 列主元法列主元法 6745150710623321xxx第一列中绝对值最大是第一列中绝对值最大是10，取，取10为主元为主元 6515707104623 6515462370710n阶方程组第阶方程组第k 轮消元时，选第轮消元时，选第k 列的后列的后( n k + 1 )个个元素中绝对值最大作主元。元素中绝对值最大作主元。 5 . 255 . 21 . 661 . 070710 2 . 62 . 65 . 255 . 270710 x3=6.2/6.2=1x2=(2.5-5x3)/2.5= -1x1=(7+7x2-0 x3)/10=0 x1=0 x2= -1x3=

8、1第二列的后两个数中选出主元第二列的后两个数中选出主元 2.5v 全主元消去法全主元消去法在第在第k 步消去前，步消去前， 在系数矩阵右下角的在系数矩阵右下角的n- k+1阶阶主子阵中，主子阵中，选绝对值最大的元素作为主元素选绝对值最大的元素作为主元素。(1) If p k then 交换第交换第 k 行与第行与第p 行行; If q k then 交换第交换第 k 列与第列与第 q 列列;(2) 消元消元列交换改变了列交换改变了xi 的顺序，须记录的顺序，须记录交换次序交换次序，解完后再换回来。解完后再换回来。,|m ax |0pqkkijkijnaa12122131233219/401/1

9、015.25/4910192340201420401419102314554155402014034.94.601/21.5254.940120404.934.601.5250.54.94012rrccllcclA 解：1320404.934.6001.433666.331614.41640,1.76514,2.35233xxx 回代：12312312310192320404455xxxxxxxxx例例2.1.3 全主元解方程组全主元解方程组: 运算量运算量 （Amount of Computation）（1 1）用克莱姆用克莱姆（CramerCramer）法则求解法则求解n n阶线性方程组阶

10、线性方程组每个行列式由每个行列式由n!项相加项相加, 而每项包含了而每项包含了n个因子相乘个因子相乘,乘法运算次数为乘法运算次数为(n-1) n!次次.,1, 2,.,iiDxinD仅考虑乘仅考虑乘(除除)法运算法运算, 计算解向量包括计算计算解向量包括计算n+1个行个行列式和列式和n次除法运算次除法运算, 乘乘(除除)法运算次数法运算次数N=(n+1)(n-1)n!+n.（2） 高斯消去法高斯消去法:在第在第1个消去步个消去步, 计算计算 li1(i=2,3,n), 有有n-1次除法运算次除法运算. 使使aij(1)变为变为 aij(2) 以及使以及使bi(1)变为变为bi(2)有有n(n-

11、1)次乘法运算次乘法运算和和 n(n-1)次加次加(减减)法运算法运算.在第在第k个消去步个消去步, 有有n-k次除法运算次除法运算, (n-k+1)(n-k)次乘次乘法运算和相同的加法运算和相同的加(减减)法运算法运算.首先统计乘法运算总次数首先统计乘法运算总次数. . 将每个消去步的乘法运算将每个消去步的乘法运算次数相加次数相加, ,有有 n(n-1)+(n-1)(n-2)+32+21=n(n-1)(n+1)/3加加(减减)法运算次数总计也为法运算次数总计也为n(n-1)(n+1)/3.除法运算总次数为除法运算总次数为n+(n-1)+1= n(n-1)/2回代过程的计算回代过程的计算除法运

12、算次数为除法运算次数为n次次. 乘法运算和加法运算的总次数乘法运算和加法运算的总次数都为都为n+(n-1)+1= n(n-1)/2次次 Gauss Gauss消去法消去法( (顺序消去法顺序消去法) ) 除法运算次数为除法运算次数为: n(n-1)/2+n= n(n+1)/2， 乘法运算次数为乘法运算次数为: n(n-1)(n+1)/3+n(n-1)/2= n(n-1)(2n+5)/6， 加加(减减)法运算次数为法运算次数为: n(n-1)(2n+5)/6通常也说通常也说Gauss消去法的运算次数与消去法的运算次数与n3同阶同阶,记为记为O(n3)u全主远消去法全主远消去法:比比 高斯消去法多

13、出高斯消去法多出 , 保证稳定保证稳定, 但费时但费时. 33nOu 列主元消去法列主元消去法: 比比 高斯消去法只多出高斯消去法只多出 的的 , 略省时略省时. 32nO2.22.2 三角分解法三角分解法 高斯消元法的矩阵形式：高斯消元法的矩阵形式：L U 分解分解 每一步消去过程相当于左乘初等变换矩阵每一步消去过程相当于左乘初等变换矩阵L1(2)(1)(2)(1)11112111113111 , 1101 2 3001L()()iin i, ,nbbAL ALllaall记：其中(3)(2)(3)(2)222222223222 , 10101 3 4001L()()iini, ,nbbAL

14、 ALlaall记：(3)(1)(3)(1)2121 , bbAL L AL L11,1,11010111= 01010101 iiiiiininiLiiLllll列列( )(1)121( )(1)121nnnnnn AL LL AbbL LL(1)111( )( )121nnnLLUAL LL AAA 的的 LU 分解分解( LU factorization )定理定理2.2.1： 若若A的所有顺序主子式的所有顺序主子式 均不为均不为0，则，则 A 的的 LU 分解唯一（其中分解唯一（其中 L 为为单位单位下三角阵）。下三角阵）。证明：证明：由由1中定理可知，中定理可知，LU 分解存在。下面

15、证明唯一性。分解存在。下面证明唯一性。若不唯一，则可设若不唯一，则可设 A = L1U1 = L2U2 ，推出推出 121UU211122211LLUULL 上三角矩阵上三角矩阵对角线上为对角线上为1 1的的下三角矩阵下三角矩阵I L 为单位下三角阵而为单位下三角阵而 U 为为一般一般上三角阵的分上三角阵的分解称为解称为Doolittle 分解分解 （2）L 为一般下三角阵而为一般下三角阵而 U 为为单位单位上三角阵的分解称为上三角阵的分解称为Crout 分解分解。 Doolittle分解法分解法 ： 利用矩阵乘法，通过比较直接导出利用矩阵乘法，通过比较直接导出L 和和 U 的计算公式。的计算

16、公式。思思路路 nnnnnnnnuuullaaaa.1.11.1111211111 ),min(1jikjkkiul jia一般计算公式11i1111 , 1, , , 2, , jji jn inuala uL第一行乘第一行乘U 每一列：每一列： a11= u11, , a1n= u1nL每一行乘每一行乘U 第一列：第一列： a21 = l 21u11, , an1 = ln1u11l 21u12+ u22=a22, , l21u1n+ u2n=a2n l31u12+ l32u22=a32, , ln1u12+ ln2u22=an2L第二行乘第二行乘U 每一列：每一列：L每一行乘每一行乘U

17、第二列：第二列：注：注：该公式的特点：该公式的特点：U 的元素按行求的元素按行求, L 的元素按列求的元素按列求; 先求先求 U 的第的第k 行行, 再求再求 L 的第的第k 列列, U 和和 L 一行一列交叉计算一行一列交叉计算. 计算计算量与量与 Gauss 消去法同消去法同.1kjkr11ikirkk1 3, , j k , , n ( )/ ik 1 , , nkkjrjrkikrkrknual ulal u u 对计算 u22=a22 - l21u12, , u2n=a2n- l21u1n l32=(a32- l31u12)/u22, , ln2=(an2- ln1u12)/u22一

18、般地：一般地：LU 分解求解线性方程组 LY b , UXY AXb112221313212(1)1111121n22222nnn111 1 UXnnnnn nnnybybLYbybxyxyYxylllllluuuuuu 11-11 , 2, 3, , (1) Y iiijijjinybyybl解：ijii1( , in1, , 1 2) x nnnnnijij iyxuyxu xu 解：例例2.2.1 求矩阵的求矩阵的Doolittle分解分解 11242142612332442A 1122214161232/32442 112221013632/32442 95/192205013632/

19、32442 9000050036302442U 15/192201010012/30001L 旧元素减去左边行与顶上列向量的点积旧元素减去左边行与顶上列向量的点积 计算行不用除法计算行不用除法 计算列要除主对角元计算列要除主对角元注：注： 求解正定方程组的求解正定方程组的Cholesky方法方法(平方根法平方根法)回顾：回顾：对称正定阵对称正定阵A的几个重要性质的几个重要性质（1）A 1 亦对称正定，且亦对称正定，且 aii 0（2）A 的顺序主子阵的顺序主子阵 Ak 亦对称正定亦对称正定（3）A 的特征值的特征值 i 0 （4）A 的全部顺序主子式的全部顺序主子式 det ( Ak ) 0定

20、理定理2.2.2: 设矩阵设矩阵A对称正定，则存在唯一的对角对称正定，则存在唯一的对角元全为正的下三角阵元全为正的下三角阵G 使得使得 AGGT计算格式为计算格式为:12111,1, 2,1,2,.,kkkkkkjjkikijkjjikkkgagknagggikkng平方根法的优点：平方根法的优点：(1) 乘除法运算量比一般乘除法运算量比一般 LU分解要小得多；分解要小得多；不选主元的平方根法是数值稳定的。不选主元的平方根法是数值稳定的。(2)缺点：缺点：有有n 个开方运算。个开方运算。 为避免过多的开方运算，在更多情况下是将A分解为：A=LDLT，其中L为下三角，D为对角阵。事实上，由前面的

21、讨论可知：A=LU，其中L 为单位下三角阵，U为上三角阵。记U的元素为uij，D=diag(uii)，由于A对称正定，必有uii0(i=1,2,n)，所以U=DU*，其中*1*TTUD UA = AL=UA= LDL为单位上三角阵,由,知,因此。注：注：将矩阵将矩阵 A 作作Doolittle分解或分解或Crout分解分解, 由矩阵乘由矩阵乘法可得法可得 A 的的 LDLT 分解。分解。 解三对角方程组的解三对角方程组的追赶法追赶法 nnnnnnnfffxxxbacbacbacb212111122211定理定理2.2.3：若若 A 为为对角占优对角占优 的三对角阵，且满足的三对角阵，且满足 则

22、方程组有唯一的则方程组有唯一的LU分解。分解。0,0, 0|, 0|11 iinncaabcb第一步第一步 : 对对 A 作作Doolittle 分解分解追赶法公式的推导：追赶法公式的推导：（以四阶为例）以四阶为例）11112222223333334444112121223232334343411222133321111(1,2,3);iibcudlabcudALUlabcudlabuudl ul dudl ul dudl ul duububl ddciubl d比较对应位置的元素，得递推公式：第二步：2213324434443/;/laulaulauubl d1111212222233332

23、334444344112221333244431111yfyfl yyflyfLyflyfl yyfyfll yyfyfyfl yyfl yLyyfl yyf第三步：由：先求该过程称为该过程称为“追追”的过程。的过程。1111222233334444443334333432223211121/()/()/()/()/udxyudxyUxyudxyuxyxyuxyd xuyc xuxyc xuxyc xUxyxu第四步：再由回代解得解得：该过程称为该过程称为“赶赶”的过程。的过程。一般情形的三对角方程组计算公式：一般情形的三对角方程组计算公式：1111,1,2,1/,2,3,1,2,3,iiii

24、iiiiidcinublauinubl cin计算次序为：计算次序为：112233(1)(2)(3)iinndcublululu最最好好牢牢记记1212323434112223334412312233442244844154420121144114414411;2,4,2,2,2,2,22 2 2xxxxxxxxxxudludALUludludddulululu追赶例 . . ：法解方程解：(2,8,15,20)(2,4,7,6)(0.5,1,2,3)TTTLyyUxyx由，解得：再由，解得：直接比较等式两边的元素，可得到计算公式：直接比较等式两边的元素，可得到计算公式：第一步第一步: 对对

25、A 作作Crout 分解：分解：1122111nnnA 注注: 也可通过对也可通过对 A 作作Crout 分解进行求解分解进行求解111,(2,3, )/,(1,2,1)iiiiiiiiibaincbain第二步第二步: 追追即解：即解：fyL 111,fy1()(2,.,)iiiiifyyin第三步第三步: 赶赶即解：即解：yxU 1,(1,.,1)nniiiixyxyxin 2.3 矩阵的条件数与方程组的性态矩阵的条件数与方程组的性态预备知识预备知识范数范数 向量范数向量范数 ( vector norms )nRyx,对任意对任意定义定义1：Rn空间的空间的向量范数向量范数 | | ,对任

26、意对任意 满足下列条件满足下列条件00|;0|) 1 (xxx|) 2 (xxR|) 3 (yxyx常用向量范数：常用向量范数： 向量向量 的的Lp范数定义为范数定义为: :12( ,)TnXx xx1/1|(| ),1npppiiXxp 主要主要性质性质性质性质1:1:-x=x性质性质2 2: :x-yx-y性质性质3 3: : 向量范数向量范数x是是Rn上向量上向量x的连续函数的连续函数. .范数等价范数等价: :设设A 和和B是是R上任意两种范数，若存在上任意两种范数，若存在 常数常数 C1、C2 0 0 使得使得 , ,则称则称 A 和和B 等价等价。定理定理1 Rn 上上一切范数一切

27、范数都等价都等价。12212111|max |nniiiiniixxxxxx 向量向量 x 的的1范数范数 向量向量 x 的的2范数范数 向量向量 x 的的范数范数 矩阵范数矩阵范数 ( matrix norms )nmRBA ,定义定义3:3:对任意对任意 , ,称称| | 为为Rm n空间的空间的矩阵矩阵范数范数, 指指| |满足满足(1)-(3)(1)-(3)：00|;0|) 1 ( AAA(2) | | | |AAR对任意对任意|) 3(BABA (4) | AB | | A | | B |若还满足若还满足(4),(4),称为相容的矩阵范数称为相容的矩阵范数相容性相容性（1 1）矩阵范

28、数与）矩阵范数与矩阵矩阵范数的相容性范数的相容性: : ABAB ABAB（2 2）矩阵范数与）矩阵范数与向量向量范数相容性范数相容性设设AM，AM 是矩阵范数是矩阵范数, xR n，xv 是向量范是向量范数数. .如果满足不等式如果满足不等式: :Axv AM xv则称矩阵范数则称矩阵范数AM 与向量范数与向量范数xv 相容相容. .常用的算子范数常用的算子范数： 由向量范数由向量范数 | |p 导出关于矩阵导出关于矩阵 A Rn n 的的p范数范数: :pxpppxAxxAApx|max|max|10| 则则ppppppxAxABAAB| njijaAni1|max|1（行和范数行和范数）

29、 niijaAnj11|max|1（列和范数列和范数）)(|max2AAAT （谱范数谱范数 ( spectral norm ) ）( operator norm ),又称为从属的矩阵范数又称为从属的矩阵范数:算子算子范数范数注：注： 称为称为A的谱半径。的谱半径。1( )max|ii nA 定理定理2对任意算子范数对任意算子范数 | | | | 有有: :|)(AA 证明：证明： 由算子范数的相容性，得到由算子范数的相容性，得到|xAxA 将任意一个特征根将任意一个特征根 所对应的特征向量所对应的特征向量 代入代入u|uAuA |uu 命题命题若若A A对称，则有对称，则有: :)(|2AA

30、 证明：证明：)()(|2maxmax2AAAAT 若若 是是 A 的一个特征根，则的一个特征根，则 2 必是必是 A2 的特征根。的特征根。又：对称矩阵的特征根为实数，即又：对称矩阵的特征根为实数，即 2(A) 为非负实数，为非负实数，故得证。故得证。)()(22maxAA 对某个对某个 A 的特征根的特征根 成立成立定理定理3若矩阵若矩阵 A 对某个算子范数满足对某个算子范数满足 |A| 1，则必有则必有.IA可逆；可逆；.111 |IAA证明证明: 若不然，则若不然，则 有非零解，即存在非零向有非零解，即存在非零向量量 使得使得 ()0IA x0 x00Axx00|1|Axx|1A 1(

31、)()IA IAI11()()IAA IA11()()IAIA IA11|()|1 | |()|IAAIA 病态方程组病态方程组线性方程组的性态（误差分析）线性方程组的性态（误差分析） ( Error Analysis for Linear system of Equations )例：方程组例：方程组1226826.000018.00001xx有精确解有精确解(1,1),(1,1),对系数矩阵和右端常熟项作微小变化对系数矩阵和右端常熟项作微小变化, ,则则1226825.999998.00002xx的解为的解为(10, -2). 扰动后方程组的解面目全非。扰动后方程组的解面目全非。Def :

32、如果方程组如果方程组 Ax = b 中中, 矩阵矩阵A和右端项和右端项 b 的变化的变化|A|和和|b |微小微小, 引起解向量引起解向量 x 的变化的变化|x |很大很大, 则称则称A为关于解为关于解方程组和矩阵求逆的病态矩阵方程组和矩阵求逆的病态矩阵, 称相应的方程组为病态方程组称相应的方程组为病态方程组.反之反之, 如果如果|A|和和|b |微小微小, |x |也微小也微小, 则称则称 A为良态矩为良态矩阵阵, Ax=b 为良态方程组为良态方程组. 思考思考:求解求解 时时, A 和和 的误差对解的误差对解 有何影响有何影响？Axbbx 设设 A 精确，精确， 有误差有误差 ，得到的解为

33、，得到的解为 ，即，即bbxx1xAb1| | |xAb绝对误差放大因子绝对误差放大因子| |bAxAx又又1|Axb1| |xbAAxb相对误差放大因子相对误差放大因子()A xxbb矩阵和方程组病态程度的刻划矩阵和方程组病态程度的刻划 设设 精确，精确，A有误差有误差 ，得到的解为，得到的解为 ，即，即bA xx()()AAxxb()()A xxA xxb1()xAA xx 11| | |xAAxxAAAA()()AA xAAxb()AAxAx 1()A IAAxAx 111()xIAAAAx (只要只要 A充分小，使得充分小，使得11| | | 1)AAAA1111| | |1 | |1 | |AAAxAAAAxAAAAA 是关键是关键的误差放大因子，称为的误差放大因子，称为A的状态数的状态数(条件数条件数)，记为记为cond (A) ,|1 AA注注: : cond(cond(A) ) 与与 所取的范数有关所取的范数有关常用条件数有：常用条件数有：cond 2(A)