矩阵A紧凑格式的Doolittle分解为课件.ppt

1、第三章第三章 线性方程组的解线性方程组的解n 概述概述n 消元法消元法n 三角分解法三角分解法n 平方根法平方根法n 向量与范数向量与范数n 方程组的性态与误差分析方程组的性态与误差分析n 迭代法迭代法3.1 概述概述n 大量实际计算问题大量实际计算问题一组线性方程组一组线性方程组如何求解？如何求解？n 1. 术语术语v 非齐次线性方程组：非齐次线性方程组：v 若记：若记： 则：则： (1.1)式可表示为式可表示为 Amnx=b线性方程组的矩阵表示线性方程组的矩阵表示) 1 . 1 ( 22112222212111212111 = =+ + + += =+ + + += =+ + + +mnm

2、nmmnnnnbxaxaxabxaxaxabxaxaxaLLLLLLLLLLLLLLLLLLLLLLLLLLLL = =aaaaaaaaa m n m mnnLLM MM MM MLLLL 2 1 2 22 21 1 12 11A = =xxxnM M21x = =bbbmM M21bn 分类：据方程的个数与未知数的个数间关系，分为：分类：据方程的个数与未知数的个数间关系，分为：v mn，即方程的个数大于未知数的个数，称为，即方程的个数大于未知数的个数，称为超定方程组超定方程组，或或矛盾方程组矛盾方程组。它没有一般意义下的解，但可以求其广义。它没有一般意义下的解，但可以求其广义解。解。v m=

3、n，一般意义下的方程组，本章主要讨论的重点，一般意义下的方程组，本章主要讨论的重点n 面临的问题：面临的问题：v 方程组方程组Ax=b有没有解？有没有解？v 有多少解？有多少解？v 如何求解？如何求解？n 2. Crammer法则求解法则求解Annx=bv 方程组方程组Ax=b有唯一解的有唯一解的充要条件充要条件是：是： |A|0 A可逆可逆 A是非奇异矩阵是非奇异矩阵v Ax=b在在A可逆时，存在唯一解可逆时，存在唯一解v 结论：结论：Crammer法则计算量非常大，需算法则计算量非常大，需算n+1个个n阶行列阶行列式。如：式。如：n100，1000次次/秒的计算机要算秒的计算机要算1012

4、0年年nnninninniiiiiaabaaaabaaDADDDxLLLLM MM MM MM MM MLLLL11111111111|,|,+ +- -+ +- -= = = =n 3 线性方程组的数值解法：线性方程组的数值解法：v 直接法：直接法：1) 思想：经有限步运算思想：经有限步运算求得精确解求得精确解 的方法：舍入误差，的方法：舍入误差，因此也是近似解因此也是近似解 高斯消元法及其变形高斯消元法及其变形2) 特点：特点：它可靠且效率高，但它只适用于中小型方程组。它可靠且效率高，但它只适用于中小型方程组。v 迭代法：迭代法：1) 思想：构造适当的初始近似解向量思想：构造适当的初始近似

5、解向量x，按一定的法则，使，按一定的法则，使之逐步逼近精确解，得到一个满足精度要求的近似解。之逐步逼近精确解，得到一个满足精度要求的近似解。即是用即是用某种极限过程某种极限过程逐步逼近逐步逼近准确解的方法。准确解的方法。2) 主要方法：主要方法：Jaccobi迭代法，迭代法，Gauss-Sidel迭代法等迭代法等3.1.1 GAUSS 消元法消元法n 一一. 几种可以直接求解的线性方程组几种可以直接求解的线性方程组 对角矩阵：对角矩阵： n次除法次除法)0a(abx:iiiiii2211 = = = = =的解为的解为bAxaaaAnn 下三角矩阵：下三角矩阵： 即：即： l11x1=b1 l

6、21x1+l22x2=b2 l31x1+l32x2+l33x3=b3 li1x1+li2x2+liixi=bi ln1x1+ln2x2+lnnxn=bnbAxllllllAnnnn= = = =即：即：LMMML212221110000代代入入iiijjjiiilxlbxlxlbxlbx,11,22121221111)()( - -= =- -= =- -= = =LL 上三角矩阵：上三角矩阵：即：即：u11x1+u12x2+u1nxn=b1 u22x2+u2nxn=b2 uiixi+ui,i+1xi+1+ui,nxn=bi un-1,n-1xn-1+un-1,nxn=bn-1 un,nxn=

7、bniijjinijiinnnnnnnnnnnuxubxuxubxubx,11, 1, 111/ )(/ )(/ + += =- - - - - - -= =- -= = =LL回回代代 = =nnnnuuuuuuALMMMLL00022211211二、同解变换二、同解变换1. 初等方阵：初等方阵： = =1000001)(LMMLLMMLkkiE = =100010010001),(LLLjiE: ri rj: ri *kn 注注v 对矩阵对矩阵A实施初等实施初等行变换行变换 对对A左乘左乘初等方阵初等方阵v 初等方阵是初等方阵是可逆可逆的，多个初等方阵的的，多个初等方阵的积积仍然仍然可逆可

8、逆v 可逆阵可逆阵A经有限次初等经有限次初等行变换行变换单位矩阵单位矩阵 = =1000110001),(MLMMkjkiE: rj +k*ri - -11AEEA)AE),(列变换列变换行变换行变换，（EAn2. 对方程组对方程组Ax=b作如下的变换，解不变作如下的变换，解不变:交换两个方程的次序交换两个方程的次序一个方程的两边同时乘以一个非一个方程的两边同时乘以一个非0数数一个方程的两边一个方程的两边，再将之加到另一，再将之加到另一个方程个方程将方程组将方程组Ax=b对应对应 的增广矩阵（的增广矩阵（A, b)，作如下变换，解不变，作如下变换，解不变交换矩阵的两行交换矩阵的两行某行乘以一个

9、非某行乘以一个非0数数某行乘以一个非某行乘以一个非0数，加到另一行数，加到另一行增广矩阵第增广矩阵第i行行第第i个方程个方程对方程组对方程组Ax=bAx=b的增广矩阵作同解变换的增广矩阵作同解变换 对其增广矩阵对其增广矩阵左乘一系列的初等矩阵左乘一系列的初等矩阵三、三、Gauss消元法消元法n 消元法基本思想：消元法基本思想：对对Ax=b的增广矩阵（的增广矩阵（A|b)进行初等变换，变进行初等变换，变成可直接求解的三种形式之一，再求解成可直接求解的三种形式之一，再求解v Gauss消元法：消元法：将将A化为上三角阵化为上三角阵 回代求解回代求解= nnnnnnnbaaabaaabaaaLMMM

10、MLL21222221111211 - - - - - - -nnnnnnnnnnnnbabaabaaabaaaa000000111122122211111211LLMMMMLLn 用高斯消去法解下列线性方程组用高斯消去法解下列线性方程组n 原方程组对应的增广矩阵为：原方程组对应的增广矩阵为： = =+ + += =+ + += =+ + +52262342321321321xxxxxxxxx 6 . 0446 . 0005 . 15 . 201125642212311121. 消元步骤消元步骤n方程组方程组AX=b的矩阵表示为：的矩阵表示为：A(1) x=b(1)n初始状态：初始状态： nn

11、nnnnaaaaaaaaaLMMMLL212222111211 = = nnbbbxxxMM2121AaaaaaaaaaAnnnnnn= = )1()1(2)1(1)1(2)1(22)1(21)1(1)1(12)1(11)1(LMMMLL = =nbbbbM21)1( 第一次消元：第一次消元：消去对角线下第消去对角线下第1列元素（列元素（a110) 方法：方法：a) 消去对角线下第消去对角线下第1列元素为列元素为0 ，即，即ri-r1ai1/a11： （i=2, 3, AbaaabaaabaaaAnnnnnnn= = )1()1()1(2)1(1)1(2)1(2)1(22)1(21)1(1)1

12、(1)1(12)1(11)1(LMMMLL )2()2()2(2)2(2)2(2)2(22)1(1)1(1)1(12)1(11)2(00nnnnnnbaabaabaaaALMMMLL), 3 , 2(), 3 , 2,(), 3 , 2()1(11)1()2()1(11)1()2()1(11)1(11niblbbnjialaaniaaliiijiijijiiLLL= =- -= = =- -= = = = 第二次：第二次：若若a220，则消去对角线下第二列为，则消去对角线下第二列为0，即：即：ri-r2ai2/a22 (i=3, 4, .) )2()2()2(2)2(2)2(2)2(22)1(

13、1)1(1)1(12)1(11)2(00nnnnnnbaabaabaaaALMMMLL )3()3()3(3)3(3)3(3)3(33)2(2)2(2)2(23)2(22)1(1)1(1)1(13)1(12)1(11)3(00000nnnnnnnbaabaabaaabaaaaALMMMMLLL), 4 , 3(), 4 , 3,(), 4 , 3()2(22)2()3()2(22)2()3()2(22)2(22niblbbnjialaaniaaliiijiijijiiLLL= =- -= = =- -= = = =n 经过经过k-1次消元后，增广矩阵变为：次消元后，增广矩阵变为： + + +

14、+ + + + +- - - -+ +- - - - - - -+ +- -)()()(1,)()(1)(, 1)(1, 1)(, 1)()()(1,)()2(2)1(, 1)1(1, 1)1(, 1)1(1, 1)1(1)1(1)1(1, 1)1(1)1(1, 1)1(11)(000000knknnkknknkkkknkkkkkkkkkkknkkkkkkknkkkkkkkkkknkkkkbaaabaaabaaabaaaabaaaaaALLMMMMMLLLLLMMMMLLn 第第k次消元：次消元： 若若akk 0，则消去，则消去第第 k 列元素列元素n 即：即：ri-rkaik/ akk (

15、i=k+1, k+2, ) + + + + + + + +- - - -+ +- - - - - - -+ +- -)()()(1,)()(1)(, 1)(1, 1)(, 1)()()(1,)()2(2)1(, 1)1(1, 1)1(, 1)1(1, 1)1(1)1(1)1(1, 1)1(1)1(1, 1)1(11)(000000knknnkknknkkkknkkkkkkkkkkknkkkkkkknkkkkkkkkkknkkkkbaaabaaabaaabaaaabaaaaaALLMMMMMLLLLLMMMMLL + + + + + + + + + + + + + + +- - - - -+ +

16、- - - - - - -+ +- -+ +)1()1()1(1,)1(1)1(, 1)1(1, 1)1()()(1,)()1(1)1(, 1)1(1, 1)1(, 1)1(1, 1)1(1)1(1)1(1, 1)1(1)1(1, 1)1(11)1(00000000knknnkknkkknkkkkkkkknkkkkkkkkknkkkkkkkkkknkkkkbaabaabaaabaaaabaaaaaALLMMMMMLLLLLMMMMLL), 1(), 1,(), 1()()1()1()()1()1()()(nkiblbbnkjialaankiaalkkikkikikkjikkijkijkkkki

17、kikLLL+ += =- -= =+ += =- -= =+ += = =+ + + + +2.经经n次消元后得到同解方程组次消元后得到同解方程组A(n)x=b(n) ，且且A(n)为上三角矩阵，可逐步为上三角矩阵，可逐步回代求解回代求解 )()3(3)3(33)2(2)2(23)2(22)1(1)1(13)1(12)1(11nnnnnnaaaaaaaaaaMLLL = = )()3(3)2(2)1(1321nnnbbbbxxxxMM1 , 1,)()(,1)(,)()()(L- -= =- -= = = + += =nniaxabxabxiiinijjijiiiinnnnnn即：UX=YU

18、Y3. 算法步骤算法步骤 将方程组用增广矩阵（将方程组用增广矩阵（A|b)（aij)n(n+1)表示，注意表示，注意c语言中数组下标从语言中数组下标从0开始，故开始，故i=0,1,n-1，j=0,1,n 消元：消元： 对对k0，1，n-2（消去第（消去第k列对角线下的元素）列对角线下的元素）a) 计算计算 l i,k= a i, k /a k,kb) 第第i行第行第k行行 l i, k （i=k+1, k+2, n-1) 回代求解：回代求解： xi=(ai,n-xjai,j )/ai, i . i=n-1,.,0, j=i+1,i+2,.n-1n,.,1,j1-n.,2, 1iiijij+=+

19、=-=kkkkkkj，其中：其中：laaa4.时间复杂度时间复杂度n 消元过程中：消元过程中：v 第第k次消元次消元1) 每行均需计算第每行均需计算第k行要乘以的常数：行要乘以的常数： lik=aik/akk ，2) 消去第消去第 i 行（第行（第k 列）时，列）时，i 行各列元素加上行各列元素加上 k 行各列元素乘以行各列元素乘以常数常数 l i k ：a i j=a i j -a k j* l i k (j=k, k+1, , n)，故需做，故需做n-k+1次次乘法、加减法乘法、加减法3) 一共要做一共要做 n-k 行行4) 综上：综上：a) 第第 k 次消元，除法：次消元，除法：n-k次

20、，乘法次，乘法：（：（n-k)*(n-k+1)次，加减次，加减法：法： （n-k)*(n-k+1)次次v 整个消元过程：整个消元过程：1) 除法：除法：2) 乘法、加减法：乘法、加减法：2/ )1(11- -= =- - - -= =nnknnk3/ )1()1)(211-=+-=nnknknnkn 综上：算法时间度量为综上：算法时间度量为 O(n3)5. Gauss顺序消元法的适用条件顺序消元法的适用条件n GAUSS直接消元法的适用条件：直接消元法的适用条件：akk(k) 0v A的各阶顺序主子式均不为的各阶顺序主子式均不为 0时，有时，有 akk(k) 0 （证明略）（证明略）v A是严

21、格对角占优矩阵，则：是严格对角占优矩阵，则：akk(k) 0 （证明略）（证明略） 举例说明举例说明 = = niijijiiaa1|A的每行主对角元的绝对值的每行主对角元的绝对值同行其余元素的绝对值之和同行其余元素的绝对值之和n 例例2：用：用Gauss消元法解方程组（消元法解方程组（4位十进制机）位十进制机） 0.0001x1+x2=1 x1+x2=2 因要对阶，损失有效数字，求得因要对阶，损失有效数字，求得x1=0,x2=1 显然显然x1,x2不是方程的解，怎么办？不是方程的解，怎么办？ x1+x2=2 0.0001x1+x2=1 解之，得到解之，得到x1=1,x2=1，接近与精确解，接

22、近与精确解 3.1.2 列主元列主元Gauss消元法消元法n 当采用当采用Gauss消元法求解时，若消元法求解时，若|akk|aik|,会导致误会导致误差扩散差扩散v 第第k步消元时，步消元时，ri-rk*aik/akk=aij-akj*aik/akkv 因为因为|akk|aik|，若，若rk行的误差为行的误差为e=(ek+1,.,en+1)，则，则e将将被扩大（被扩大（aik/akk)倍倍 Gauss消元法在特殊情况下是不稳定的消元法在特殊情况下是不稳定的n 解决办法：解决办法：v 选择列主元：即选出第选择列主元：即选出第k列对角线下绝对值最大者列对角线下绝对值最大者1) 设设|a t k|

23、=max |a i k| k i n2) rk rt3) 以新的以新的 akk 作为主对角元进行消元作为主对角元进行消元v 故列主元消元法只是在每步消元之前选出主对角元故列主元消元法只是在每步消元之前选出主对角元v 举实例说明列主元举实例说明列主元Gauss消元法的计算过程消元法的计算过程n 列主元算法，只需在消元的第列主元算法，只需在消元的第1步前增加以下几行代码步前增加以下几行代码 i=k; /选择第选择第k列的主元列的主元for(j=k+1;jfabs(aik) i=j;if(i!=k)for(j=k;j0n 对称正定阵对称正定阵A A的几个重要性质的几个重要性质vA A- -1 1 亦

24、对称正定，且亦对称正定，且 a aiiii 0 0vA A 的顺序主子阵的顺序主子阵 A Ak k 亦对称正定亦对称正定vA A 的特征值的特征值 i i 0 0 vA A 的全部顺序主子式的全部顺序主子式 detdet ( ( A Ak k ) ) 0 0n 对称正定矩阵的乔累斯基分解对称正定矩阵的乔累斯基分解v 对称正定矩阵对称正定矩阵A=Rnn，存在，存在的单位下三角矩阵的单位下三角矩阵L和对角矩阵和对角矩阵D，使，使A=LDLT1) 证明：前面已经证明了对称正定矩阵证明：前面已经证明了对称正定矩阵A做做Doolittle分解的唯一性，分解的唯一性，设设A=LU（其中（其中L为单位下三角

25、矩阵，为单位下三角矩阵，U为上三角矩阵）为上三角矩阵）2) 令令，其中，其中di为为U阵对角线元素，故阵对角线元素，故di03) 则：则：A=LD ( D-1U )，AT=(D-1U) TDTLT= (D-1U ) T DLT，又，又A=AT，故：故： L D ( D-1U )=(D-1U ) D LT4) 由上式知：由上式知：L= ( D-1U ）T，LT= D-1U 5) 所以所以A=LDLT=nddD1v对称正定矩阵对称正定矩阵A=Rnn，存在，存在唯一唯一的主对角元素均为正数的下三角矩阵的主对角元素均为正数的下三角矩阵L ，使，使 A=LLT，此即：对称正定矩阵的，此即：对称正定矩阵的

26、乔累斯基分解乔累斯基分解v证明：证明：1)因因A可唯一分解为可唯一分解为A=LDLT，且，且D是是di0的对角矩阵的对角矩阵2)故可将故可将D进行拆分进行拆分3)令：令：11T22() ()ALDLD=TLLALDL21= = =，则则n 比较的等式两边，可得比较的等式两边，可得n 第第 j 列元素：列元素：-=-=11211111111jkkjjjjjiilallalal = =nnnjnni ij iiij jjjllllllllllllllLLMLMLLLMLMLLML112121222111 nnnnllllllLLLLL22211211* = =nnnnnnaaaaaaaaaALMM

27、MLL212222111211njjilllaljjjkkjkijiji,.,2, 1)(11+=-=-=n 计算顺序为：计算顺序为：l11 li 1 li 2 li nn |A|=|L|LT|=l112lnn2n 示例：对矩阵做乔累斯基分解，并求示例：对矩阵做乔累斯基分解，并求|A|n 乔累斯基分解的缺点：开方乔累斯基分解的缺点：开方3353595917A= )17()9()5()9()5()3()5()3()3(335322223改进的平方根法：改进的平方根法：n 利用已经证明的结论：对称正定矩阵利用已经证明的结论：对称正定矩阵A，可唯一分，可唯一分解为：解为：A=LDLT = =1112

28、121LLLLnnlll 111*2112LLLLLnnlll = =nnnnnnaaaaaaaaaALMMMLL212222111211 nddd21*n 故：故：d1=u11，d1l12=u12， di li k=ui kn 则实对称正定矩阵则实对称正定矩阵A=LDLT，可视为，可视为A=LU分解，由于分解，由于A的对称的对称性，致使性，致使L、U之间具有以下关系：之间具有以下关系： = =nnnnndldldlddlddlllMLLL2211232212112112*111kkikkikkkiuudldl= = = = =1112121LLLLnnlll nnnnuuuuuuLLLLL2

29、2211211*n 则写成紧凑格式为：则写成紧凑格式为：n 此即为改进平方根法：在此即为改进平方根法：在LU分解过程中，利用系数矩阵的对称分解过程中，利用系数矩阵的对称性，作性，作Doolittle分解，只需计算分解，只需计算U矩阵，矩阵，L矩阵第矩阵第 k 列的值可由列的值可由U矩阵矩阵k行的值除以行的值除以ukk算得，使得算得，使得LU分解的计算量减少了一半分解的计算量减少了一半n 其他可完全按其他可完全按Doolittle分解求解线性方程组的解的方法进行分解求解线性方程组的解的方法进行 = =- - - -+ +- - - - - - - -MMMM1, 1, 12213111, 1,1

30、, 1, 1, 1, 11, 1222311132232211121131211kknknkkknkkkkkknkkkkknnuuuuuuluuuuuuuuuuuuuuuuuuuun 示例示例3-17：用改进平方根法解线性方程组：用改进平方根法解线性方程组n 解：因系数矩阵是对称正定矩阵（需判断），用改进平方根解：因系数矩阵是对称正定矩阵（需判断），用改进平方根法对方程组的增广矩阵做法对方程组的增广矩阵做Doolittle分解分解 - -= = 4201795953533321xxx - -= =)4()17()9()3()2()9()5()3()0()5()3()3(A335 015/3224-22/30yn 则则A=LU，Ax=bLUx=b，令，令y=Ux，则，则Ly=b，A矩阵改进平矩阵改进平方根法方根法LU分解的最后一列即为分解的最后一列即为y向量向量n 所以解所以解Ux=yn 得：得：x3=0，x2=-1，x1=1ULA = =3242533*1235111 - -= = = =0203242533321xxxUx