Ch 解线性方程组的直接方法

1、ch5 解线性方程组的直接方法解线性方程组的直接方法求解求解bxa 高斯消元法：高斯消元法：思思路路首先将首先将a化为上三角阵化为上三角阵 /* upper-triangular matrix */，再回代求解再回代求解 /* backward substitution */。=消元消元记记,)()1()1(nnijaaa )1()1(1)1(.nbbbbstep 1：设设 ，计算因子，计算因子0)1(11 a).,2(/)1(11)1(11niaamii 将增广矩阵将增广矩阵/* augmented matrix */ 第第 i 行行 mi1 第第1 1行行，得到，得到)1(1)1(1)1(

2、12)1(11.baaan)2(a) 2(b).,2,()1(11)1()2()1(11)1()2(njibmbbamaaiiijiijij 其中其中step k：设设 ，计算因子，计算因子且计算且计算0)( kkka)., 1(/)()(nkiaamkkkkikik )., 1,()()()1()()()1(nkjibmbbamaakkikkikikkjikkijkij 共进行共进行 ？ 步步n 1 )()2(2)1(121)()2(2)2(22)1(1)1(12)1(11.nnnnnnnnbbbxxxaaaaaa回代回代)()(/nnnnnnabx 0)( nnnawhat if ?no

3、unique solution exists.) 1., 1()(1)()( niaxabxiiinijjiijiii0)( iiiawhat if ?then we must find the smallest integer k i with , and interchange the k-th row with the i-th row. 0)( ikiawhat if we cant find such k ?no unique solution exists.定理定理 若若a的所有的所有顺序主子式顺序主子式 /* determinant of leading principal su

4、bmatrices */ 均不为均不为0，则高斯消元无需换行即可，则高斯消元无需换行即可进行到底，得到唯一解。进行到底，得到唯一解。iiiiiaaaaa.)det(1111 事实上，只要事实上，只要 a 非奇异，即非奇异，即 a 1 存在，则可通过逐存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出次消元及行交换，将方程组化为三角形方程组，求出唯一解。唯一解。1 gaussian elimination the method0(1,2,1)定理（矩阵的 分解）定理（矩阵的 分解）设 为 阶矩阵，如果 的顺序主子式，设 为 阶矩阵，如果 的顺序主子式，则 可分解为一个下三角矩阵 和一

5、个上三角矩阵 的乘积，则 可分解为一个下三角矩阵 和一个上三角矩阵 的乘积，且这种分解是唯一的.且这种分解是唯一的.iluanadinalu 选主元消去法选主元消去法例：例：单精度解方程组单精度解方程组 211021219xxxx/* 精确解为精确解为 和和 */.1000.00. 1101191 x8个个.8999.99. 0212 xx8个个用用gaussian elimination计算：计算：911212110/ aam999212210101010.0 . 011 ma8个个92121012 mb 9991010011100, 112 xx小主元小主元 /* small pivot

6、element */ 可能导致计可能导致计算失败。算失败。 全主元消去法全主元消去法 /* complete pivoting */每一步选绝对值最大的元素为主元素，保证每一步选绝对值最大的元素为主元素，保证 。1| ikmstep k: 选取选取;0|max|, ijnjikjiaakk if ik k then 交换第交换第 k 行与第行与第 ik 行行; if jk k then 交换第交换第 k 列与第列与第 jk 列列; 消元消元列交换改变了列交换改变了 xi 的顺序，须记录的顺序，须记录交换次序交换次序，解完后，解完后再换回来。再换回来。 列主元消去法列主元消去法 /* parti

7、al pivoting, or maximal column pivoting */省去换列的步骤，每次仅选一列中最大的元。省去换列的步骤，每次仅选一列中最大的元。0|max|, iknikkiaak例：例： 211111091,112 xx 110211 11102119 列主元法没有全主元法稳定。列主元法没有全主元法稳定。0,112 xx例：例： 2111010199 99991010010101注意：这两个方程组注意：这两个方程组在数学上在数学上严格等价严格等价。 标度化列主元消去法标度化列主元消去法 /* scaled partial pivoting */对每一行计算对每一行计算 。

8、为省时间，。为省时间，si 只在初始时计只在初始时计算一次。以后每一步考虑子列算一次。以后每一步考虑子列 中中 最大的最大的 aik 为主元。为主元。|max1ijnjias nkkkaa.iiksa稳定性介于列主元法和全主元法之间。稳定性介于列主元法和全主元法之间。1 gaussian elimination pivoting strategies 高斯高斯- -若当若当消去法消去法 /* gauss-jordan method */与与 gaussian elimination 的主要区别：的主要区别： 每步不计算每步不计算 mik ，而是先将当前主元，而是先将当前主元 akk(k) 变为

9、变为 1; 把把 akk(k) 所在列的上、下元素全消为所在列的上、下元素全消为0;baxibxa1 hey! isnt it better than gaussian elimination?what makes you say so?obviously we no longer need the backward substitution!youd better wait till we go through the next section to draw your conclusion1 gaussian elimination gauss-jordan method 运算量运算量 /

10、* amount of computation */1 gaussian elimination amount of computation由于计算机中乘除由于计算机中乘除 /* multiplications / divisions */ 运算的时运算的时间远远超过加减间远远超过加减 /* additions / subtractions */ 运算的时间，故运算的时间，故估计某种算法的运算量时，往往只估计估计某种算法的运算量时，往往只估计乘除的次数乘除的次数，而且通，而且通常以乘除次数常以乘除次数的的最高次幂最高次幂为运算量的为运算量的数量级数量级。 gaussian eliminatio

11、n:step k：设设 ，计算因子，计算因子且计算且计算0)( kkka)., 1(/)()(nkiaamkkkkikik )., 1,()()()1()()()1(nkjibmbbamaakkikkikikkjikkijkij 共进行共进行n 1步步 )()2(2)1(121)()2(2)2(22)1(1)1(12)1(11.nnnnnnnnbbbxxxaaaaaa)()(/nnnnnnabx ) 1., 1()(1)()( niaxabxiiinijjiijiii(n k) 次次(n k)2 次次(n k) 次次(n k) (n k + 2) 次次nnnknknnk6523)2)(2311

12、 消元乘除次数：消元乘除次数：1 次次(n i +1) 次次22)1(1211nninni 回代乘除次数：回代乘除次数：gaussian elimination 的总乘除次数为的总乘除次数为 ，运算量为，运算量为 级。级。nnn31323 33n1 gaussian elimination amount of computation complete pivoting:比比 gaussian elimination多出多出 比较，保证稳定，但费时。比较，保证稳定，但费时。 33no partial pivoting:比比 gaussian elimination只多出只多出 比较，略省时，但不

13、保比较，略省时，但不保证稳定。证稳定。 32no scaled partial pivoting:比比 gaussian elimination多出多出 除法和除法和 比较，比列主比较，比列主元法稳定。但若逐次计算元法稳定。但若逐次计算 si(k)，则比全主元法还慢。，则比全主元法还慢。)(2no 32no gauss-jordan method:运算量约为运算量约为 。故通常只用于求逆矩阵，而不用于解方。故通常只用于求逆矩阵，而不用于解方程组。求逆矩阵即程组。求逆矩阵即 。 23no 1| aiiahw: p.42 #1 p.43 #62 三角分解法三角分解法 /* matrix facto

14、rization */ 高斯消元法的矩阵形式高斯消元法的矩阵形式 /* matrix form of g.e. */：step 1:)0(/111111 aaamii记记 l1 =1.11121nmm ，则，则 )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 =1.11, 1knkkmm 2 matrix factorization matrix form of g.e. 1kl1.11, 1knkk

15、mm 111211.nlll111jim,记为记为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 different with a fixed a.bco

16、uld you be more specific, please?factorize a first, then for every you only have to solve two simple triangular systems and.bbyl yxu 2 matrix factorization matrix form of g.e.定理定理 若若a的所有顺序主子式的所有顺序主子式 /* determinant of leading principal submatrices */ 均不为均不为0，则，则 a 的的 lu 分解唯一（其分解唯一（其中中 l 为为单位单位下三角阵）。

17、下三角阵）。证明：证明：由由1中定理可知，中定理可知，lu 分解存在。下面证明唯一性。分解存在。下面证明唯一性。若不唯一，则可设若不唯一，则可设 a = l1u1 = l2u2 ，推出，推出 121uu211122211lluull upper-triangularlower-triangularwith diagonal entries 1i l 为一般下三角阵而为一般下三角阵而 u 为为单位单位上三角阵的分解称为上三角阵的分解称为crout 分解分解。 实际上只要考虑实际上只要考虑 a* 的的 lu 分解，即分解，即 ，则，则 即是即是 a 的的 crout 分解。分解。ula* *lua

18、 2 matrix factorization doolittle 道立特分解法道立特分解法 /* doolittle factorization */： lu 分解的紧凑格式分解的紧凑格式 /* compact form */反复计算反复计算,很浪费哦很浪费哦 通过比较法直接导出通过比较法直接导出l 和和 u 的计算公式。的计算公式。思思路路 nnnnnnnnuuullaaaa.1.11.1111211111 ),min(1jikjkkiul jia2 matrix factorization choleski 平方根法平方根法 /* choleskis method */： 对称对称 /*

19、 symmetric */ 正定正定 /* positive definite */ 矩阵的分解法矩阵的分解法定义定义一个矩阵一个矩阵 a = ( aij )n n 称为称为对称阵对称阵，如果，如果 aij = aji 。定义定义一个矩阵一个矩阵 a 称为称为正定阵正定阵，如果，如果 对任意非对任意非零向量零向量 都成立。都成立。0 xaxtx回顾：回顾：对称正定阵的几个重要性质对称正定阵的几个重要性质 a 1 亦对称正定，且亦对称正定，且 aii 0若不然，则若不然，则0 xa存在非零解，即存在非零解，即0 xaxt存在非零解。存在非零解。1111)()(, aaiaaiaattt对任意对任

20、意 , 存在存在 , 使得使得 ,即即 。0 x0 yxya xay1 011 yayyaaayxaxttt 其中其中0 xaxatiitx)0.1.0( 第第 i 位位 a 的顺序主子阵的顺序主子阵 /* leading principal submatrices */ ak 亦对亦对 称正定称正定对称性显然。对任意对称性显然。对任意 有有 , 其中其中 。kkrx 00 xaxxaxtkktknkrxx 00 a 的特征值的特征值 /* eigen value */ i 0 设对应特征值设对应特征值 的非零特征向量的非零特征向量为为 ，则，则 。20 xxxxaxtt x a 的全部顺序主

21、子式的全部顺序主子式 det ( ak ) 0因为因为 niia1)det( 2 matrix factorization 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 对角元为正，则分解唯一。对角

22、元为正，则分解唯一。tlla 对于对称正定阵对于对称正定阵 a ，从，从 可知对任意可知对任意k i 有有 。即即 l 的元素不会增大，误差可控，不的元素不会增大，误差可控，不需选主元。需选主元。 ikikiila12iiikal |2 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

23、 aij for 1 i, j n of a.output: the entries 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 112ikikiiiila

24、l iiikikjkjijilllal 11 112nknknnnnlal因为因为a 对称，所以只需存半个对称，所以只需存半个 a，即，即其中其中 nnnaaaaanna.,.,2/)1(1222111 jiiaaij 2/)1(运算量为运算量为 o(n3/6)， 比普通比普通lu分解少一半，但有分解少一半，但有 n 次开方。用次开方。用a = ldlt 分解，可省开方时间分解，可省开方时间(p.50-51)。hw: p.54 #2, #5, #62 matrix factorization tridiagonal system 追赶法解追赶法解三对角三对角方程组方程组 /* crout re

25、duction for tridiagonal linear system */ nnnnnnnfffxxxbacbacbacb212111122211step 1: 对对 a 作作crout 分解分解 111121nnna 直接比较等式两边直接比较等式两边的元素，可得到计的元素，可得到计算公式算公式（p.52）。step 2: 追追即解即解 ：fyl ,111 fy ),.,2()(1niyrfyiiiii step 3: 赶赶即解即解 ：yxu )1,., 1(,1 nixyxyxiiiinn 与与g.e.类似，一旦类似，一旦 i = 0 则算法中断，故并非任何则算法中断，故并非任何三对角

26、阵都可以用此方法分解。三对角阵都可以用此方法分解。2 matrix factorization 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.

27、well, how large is large? they satisfy the following inequality: ijijiiaa|如果如果 a 是是严格严格对角占优阵，则不要求三对角线上的对角占优阵，则不要求三对角线上的所有元素非零。所有元素非零。根据不等式根据不等式 可知：分解过程中，矩阵元素不会过分增大，算法可知：分解过程中，矩阵元素不会过分增大，算法保证稳定。保证稳定。运算量为运算量为 o(6n)。|,1|1iiiiiiiiabbab 矩阵范数矩阵范数( ) |,(1)| 0,| 00;(2)| | |;(3)| |;(4)| | |( ).定义如果矩阵的某个非负的实数

28、函数定义如果矩阵的某个非负的实数函数满足条件：满足条件：且且则称是上的一个矩阵范数则称是上的一个矩阵范数n nn narn aaaaacacaababababn ar 0,| ,| |max| | |定理设若给出一种向量范数定理设若给出一种向量范数定义非负函数定义非负函数则是上矩阵范数，且满足相容性条件则是上矩阵范数，且满足相容性条件nn nvvvxvn nvvvvxrarxaxaxaraxax 111112max1,2,|max|()|max|()|()(2)定理当时定理当时行范数行范数列范数列范数范数范数niji njnijj nitvaaaaaa a 12,34例设计算 的各种范数例设计

29、算 的各种范数aa 1(1,2, ) ( )max|.定义设的特征值为，称定义设的特征值为，称为 的谱半径为 的谱半径n niii narinaa ,( ) |.定理设则定理设则n naraa 2|( ).定理如果为对称矩阵，则定理如果为对称矩阵，则n naraa 1| 1,1|()|1 |定理如果则为非奇异矩阵，且定理如果则为非奇异矩阵，且bibibb 误差分析误差分析矩阵的条件数矩阵的条件数1212,11211.0001211211.00012.0001考察方程组当 或 有微小扰动时对解的影响考察方程组当 或 有微小扰动时对解的影响如方程组如方程组与与axbabxxxx 精确解精确解 (2

30、, 0)精确解精确解 (1, 1)定义如果矩阵 或常数项 的微小变化，引起方程组定义如果矩阵 或常数项 的微小变化，引起方程组解的巨大变化，则称此方程组为病态方程组，矩阵 称为解的巨大变化，则称此方程组为病态方程组，矩阵 称为“病态”矩阵，否则称为“良态”方程组和“良态”矩阵.“病态”矩阵，否则称为“良态”方程组和“良态”矩阵.abaxba 下面两个定理描述了当系矩阵或是常数项下面两个定理描述了当系矩阵或是常数项有微小变化对解的影响有微小变化对解的影响10,(), | |定理设 是非奇异矩阵，且定理设 是非奇异矩阵，且则则aaxba xxbbxbaaxb 1110,()().| | 1,| |

31、1 | |定理设 为非奇异矩阵，且定理设 为非奇异矩阵，且如果则如果则aaxbaa xxbaaaaaxaaxaaa 1 ( )| | (1,2,).定义设 为非奇异矩阵，称数定义设 为非奇异矩阵，称数为矩阵 的条件数为矩阵 的条件数vvvacond aaava 1max222min1212222()(1) ( )| |;()( ),(2)( )1(3)()( )(4)( )1 ()()( )性质性质当 为对称矩阵时，当 为对称矩阵时，其中 ，为绝对值最大和绝对值最小的特征值.其中 ，为绝对值最大和绝对值最小的特征值.；若 为正交矩阵，则 ；如果 为正交阵，则若 为正交矩阵，则 ；如果 为正交阵

32、，则ttnnvvva acond aaaa aacond acond acond cacond aacond arcond racond arcond a 矩阵的正交三角化及应用矩阵的正交三角化及应用2112122122212112,1,( )1212222122 2212(,) .定义设向量且称矩阵 定义设向量且称矩阵 为初等反射阵或豪斯霍尔德（househol der）变换，为初等反射阵或豪斯霍尔德（househol der）变换，其中 其中 nttnnnntnwrw wh wwwww ww ww www ww ww wwww ww -1112,1,(1)(2);(3).定理设有初等反射阵其中则：定理设有初等反射阵其中则：是对称矩