数值计算：线性方程组求解(2)

1、4 三角分解法三角分解法 /* Matrix Factorization */ 高斯消元法的矩阵形式高斯消元法的矩阵形式 /* Matrix Form of G.E. */：Step 1:)0(/111111 aaalii1.11121nll 记记 L1 =，则，则 )1()1(1bAL)1(1)1(1)1(11.baan)2(A) 2 (bStep n 1: )()2(2) 1(1)()2(2)2(22) 1(1) 1(12) 1(11121.nnnnnnnnnbbbaaaaaabALLL其中其中 Lk =111,1knkkll 4 Matrix Factorization Matrix F

2、orm of G.E. 1kL1.11,1knkkll 111211.nLLL111jil,记为记为L单位下三角阵单位下三角阵/* unitary lower-triangular matrix */记记 U =)()2(2)2(22)1(1)1(12)1(11.nnnnnaaaaaaLUA A 的的 LU 分解分解/* LU factorization */Hey hasnt GE given me enough headache? Why do I have to know its matrix form?!When you have to solve the system for dif

3、ferent with a fixed A.bCould you be more specific, please?Factorize A first, then for every you only have to solve two simple triangular systems and.bbyL yxU 4 Matrix Factorization Matrix Form of G.E.证明：证明：由由2中定理可知，中定理可知，LU 分解存在。下面证明唯一性。分解存在。下面证明唯一性。若不唯一，则可设若不唯一，则可设 A = L1U1 = L2U2 ，推出，推出 121UU21112

4、2211LLUULL Upper-triangularLower-triangularWith diagonal entries 1I L 为一般下三角阵而为一般下三角阵而 U 为为单位单位上三角阵的分解称为上三角阵的分解称为Crout 分解分解。 实际上只要考虑实际上只要考虑 AT 的的 LU 分解，即分解，即 ，则，则 即是即是 A 的的 Crout 分解。分解。ULAT TTLUA 定理定理 设设A为为n阶矩阵，如果阶矩阵，如果A的顺序主子式的顺序主子式 /* determinant of leading principal submatrices */ ，则，则 A 的的 LU 分解唯

5、一（其中分解唯一（其中 L 为为单位单位下三角阵）。下三角阵）。)1, 2 , 1(0 niDi考虑考虑A非奇异，奇异时见教材非奇异，奇异时见教材p170的证明。的证明。4 Matrix Factorization Doolittle 道立特分解法道立特分解法 /* Doolittle Factorization */： LU 分解的紧凑格式分解的紧凑格式 /* compact form */反复计算反复计算,很浪费哦很浪费哦 通过比较法直接导出通过比较法直接导出L 和和 U 的计算公式。的计算公式。思思路路 nnnnnnnnuuullaaaa.1.11.1111211111 ),min(1j

6、ikjkkiul jia4 Matrix Factorization Doolittle ),min(1jikjkkijiula固定固定 i : :对对 j = i, i+1, , n 有有ijkjikikijuula 11lii = 1kjikikijijulau 11a将将 i ，j 对换，对对换，对 j = i, i+1, , n 有有iijikiikjkjiulula 11iiikkijkjijiuulal/ )(11 b Algorithm: Doolittle FactorizationStep 1: u1j = a1j; lj1 = aj1 / u11; ( j = 1, , n

7、 )Step 2: compute and for i = 2, , n 1;Step 3: ab 11nkknnknnnnulau一般采用一般采用列主元列主元法增强稳定性。但注意法增强稳定性。但注意 也必须做相应的也必须做相应的行交换。行交换。b例：例： 131321112A解：解： 131321112A 10542725112157211211再求第一列212 先求第一行213 )25()1(213 27*57211 kjikikijijulau 11iiikkijkjijiuulal/ )(11 4 Matrix Factorization Doolittle例例 用直接三角分解法求解用

8、直接三角分解法求解 201814513252321321xxx解解：用分解公式计算得用分解公式计算得LUA 2400410321153012001求解求解TLY)20,18,14( ，得，得TY)72,10,14( TUX)72,10,14( ，得，得TX)3 , 2 , 1( 4 Matrix Factorization Doolittle4 Matrix Factorization Choleski 平方根法平方根法 /* Choleskis Method */： 对称对称 /* symmetric */ 正定正定 /* positive definite */ 矩阵的分解法矩阵的分解法定

9、义定义一个矩阵一个矩阵 A = ( aij )n n 称为称为对称阵对称阵，如果，如果 aij = aji 。定义定义一个矩阵一个矩阵 A 称为称为正定阵正定阵，如果，如果 对任意非对任意非零向量零向量 都成立都成立。0 xAxTx回顾：回顾：对称正定阵的几个重要性质对称正定阵的几个重要性质 A 1 亦对称正定，且亦对称正定，且 aii 0若不然，则若不然，则0 xA存在非零解，即存在非零解，即0 xAxT存在非零解。存在非零解。1111)()(, AAIAAIAATTT对任意对任意 , 存在存在 , 使得使得 ,即即 。0 x0 yxyA xAy1 011 yAyyAAAyxAxTTT 其中

10、其中0 xAxaTiiTx)0.1.0( 第第 i 位位 A 的顺序主子阵的顺序主子阵 /* leading principal submatrices */ Ak 亦对亦对 称正定称正定对称性显然。对任意对称性显然。对任意 有有 , 其中其中 。kkRx 00 xAxxAxTkkTknkRxx 00 A 的特征值的特征值 /* eigen value */ i 0 设对应特征值设对应特征值 的非零特征向量的非零特征向量为为 ，则，则 。20 xxxxAxTT x A 的全部顺序主子式的全部顺序主子式 det ( Ak ) 0因为因为 niiA1)det( 4 Matrix Factoriza

11、tion Choleski将将对称对称 正定阵正定阵 A 做做 LU 分解分解U =uij=u11uij / uii111u22unn记为记为UD A 对称对称TUL 即即TLDLA 记记 D1/2 =11u22unnuWhy is uii 0?Since det(Ak) 02/1LDL 则则 仍是下三角阵仍是下三角阵TLLA nnRL 定理定理 设矩阵设矩阵A对称正定，则存在非奇异下三角阵对称正定，则存在非奇异下三角阵 使得使得 。若限定。若限定 L 对角元为正，则分解唯一。对角元为正，则分解唯一。TLLA 对于对称正定阵对于对称正定阵 A ，从，从 可知对任意可知对任意k i 有有 。即即

12、 L 的元素不会增大，误差可控，不的元素不会增大，误差可控，不需选主元。需选主元。 ikikiila12iiikal |TTTTDLUULDAALUULD)( 4 Matrix Factorization Choleski Algorithm: Choleskis MethodTo factor the symmetric positive definite n n matrix A into LLT, where L is lower triangular.Input: the dimension n; entries aij for 1 i, j n of A.Output: the en

13、tries lij for 1 j i and 1 i n of L. Step 1 Set ;Step 2 For j = 2, , n, set ;Step 3 For i = 2, , n 1, do steps 4 and 5Step 4 Set ; Step 5 For j = i+1, , n, set ; Step 6 Set ;Step 7 Output ( lij for j = 1, , i and i = 1, , n ); STOP. 1111al 1111/laljj 112ikikiiiilal iiikikjkjijilllal 11 112nknknnnnlal

14、因为因为A 对称，所以只需存半个对称，所以只需存半个 A，即，即其中其中 nnnaaaaannA.,.,2/)1(1222111 jiiAaij 2/)1(运算量为运算量为 O(n3/6)， 比普通比普通LU分解少一半，但有分解少一半，但有 n 次开方。用次开方。用A = LDLT 分解，可省开方时间分解，可省开方时间(p.173)。HW: p.203 #9, #10, #114 Matrix Factorization Tridiagonal System 追赶法解追赶法解三对角三对角方程组方程组 /* Crout Reduction for Tridiagonal Linear Syste

15、m */ nnnnnnnfffxxxbacbacbacb212111122211Step 1: 对对 A 作作Crout 分解分解 111121nnnA 直接比较等式两边直接比较等式两边的元素，可得到计的元素，可得到计算公式算公式。Step 2: 追追即解即解 ：fyL ,111 fy ),.,2()(1niyrfyiiiii Step 3: 赶赶即解即解 ：yxU )1,., 1(,1 nixyxyxiiiinn 与与G.E.类似，一旦类似，一旦 i = 0 则算法中断，故并非任何则算法中断，故并非任何三对角阵都可以用此方法分解。三对角阵都可以用此方法分解。4 Matrix Factoriz

16、ation Tridiagonal System定理定理 若若 A 为为对角占优对角占优 /* diagonally dominant */ 的三对角的三对角阵，且满足阵，且满足 ，则追赶，则追赶法可解以法可解以 A 为系数矩阵的方程组。为系数矩阵的方程组。0,0, 0|, 0|11 iinncaabcbHey, what does diagonally dominant mean? It means that the diagonal entries of the matrix are very LARGE.Well, how large is LARGE? They satisfy the

17、 following inequality: ijijiiaa|F 如果如果 A 是是严格严格对角占优阵，则不要求三对角线上的对角占优阵，则不要求三对角线上的所有元素非零。所有元素非零。根据不等式根据不等式 可知：分解过程中，矩阵元素不会过分增大，算法可知：分解过程中，矩阵元素不会过分增大，算法保证稳定。保证稳定。运算量为运算量为 O(5n)。|,1|1iiiiiiiiabbab 1111,iiiiiiiiiccbaab 5 线性方程组的误差分析线性方程组的误差分析 /* Error Analysis for Linear system of Equations */求解求解 时，时，A 和和

18、 的误差对解的误差对解 有何影响？有何影响？bxA bx 设设 A 精确，精确， 有误差有误差 ，得到的解为，得到的解为 ，即，即bb xx bbxxA )(bAx 1 |1bAx 绝对误差放大因子绝对误差放大因子|xAxAb 又又|1bAx |1bbAAxx 相对误差放大因子相对误差放大因子5 Error Analysis for . bxA 设设 精确，精确，A有误差有误差 ，得到的解为，得到的解为 ，即，即bA xx bxxAA )( bxxAxxA )()( )(1xxAAx |11AAAAAAxxx bxAAxAA )()(xAxAA )(xAxAAIA )(1xAAAAIx 111

19、)( Wait a minute Who said that ( I + A 1 A ) is invertible?|1|1|1111AAAAAAAAAAAAxx 是关键是关键的误差放大因子，称为的误差放大因子，称为A的的条件数条件数，记为，记为cond (A) ,越越 则则 A 越病态，越病态，难得准确解。难得准确解。|1 AA大大(只要只要 充分小，使得充分小，使得)1|11 AAAA | A 5 Error Analysis for . bxA 注注: FF cond (A) 的具体大小与的具体大小与 | | 的取法有关，但相的取法有关，但相对大小一致。对大小一致。FF cond (A

20、) 取决于取决于A，与解题方法无关。，与解题方法无关。 |)(1)(|bbAAAAAcondAcondxx FF 常用条件数有：常用条件数有：cond (A)1cond (A) cond (A)2)(/ )(minmaxAAAATT 特别地，若特别地，若 A 对称，则对称，则|min|max)( 2 Acond条件数的性质：条件数的性质： A可逆，则可逆，则 cond (A)p 1； A可逆，可逆， R 则则 cond ( A) = cond (A) ； A正交，则正交，则 cond (A)2=1； A可逆，可逆，R正交，则正交，则 cond (RA)2 = cond (AR)2 = cond

21、 (A)2 。5 Error Analysis for . bxA 精确解精确解为为.11 x例：例： 97.199.1,98.099.099.01bA计算计算cond (A)2 。 10000990099009800A 1 = 解：解：考察考察 A 的特征根的特征根 0)det(AI 000050504. 0980050504. 121 212)( Acond 39206 1 测试病态程度：测试病态程度：给一个扰动给一个扰动b 3410106.01097.0b ，其相对误差为，其相对误差为%01.010513.0|422 bb 此时此时精确解精确解为为 0203. 13*x 0203.22*xxx 22|x