1、1 高斯消去法高斯消去法 主元素法主元素法 直接三角分解法直接三角分解法 平方根法与改进平方根法平方根法与改进平方根法 误差分析误差分析第二章第二章 线性方程组的直接方法线性方程组的直接方法 2 .,.,.22112222212111212111nnnnnnnnnnbxaxaxabxaxaxabxaxaxa讨论讨论 n元线性方程组元线性方程组的直接解法的直接解法.3其中其中 ,.212222111211 nnnnnnaaaaaaaaaA,21 nxxxx.21 nbbbb 方程组的矩阵形式为方程组的矩阵形式为,bAx 若矩阵若矩阵A非奇异非奇异,即即det(A)0,则方程组有唯一解则方程组有唯
2、一解.4 直接解法是指直接解法是指,若不考虑计算过程中的舍入误若不考虑计算过程中的舍入误差，经过差，经过有限次算术运算就能求出线性方程组的有限次算术运算就能求出线性方程组的精确解精确解的方法的方法.但由于实际计算中舍入误差的存但由于实际计算中舍入误差的存在在,用直接解法一般也只能求出方程组的近似解用直接解法一般也只能求出方程组的近似解.Cramer法则是一种不实用的直接法，下面介法则是一种不实用的直接法，下面介绍几种实用的直接法绍几种实用的直接法.解线性方程组的方法解线性方程组的方法:直接方法和迭代法直接方法和迭代法.迭代法是从解的某个近似值出发迭代法是从解的某个近似值出发,通过构造一通过构造
3、一个无穷序列去逼近精确解的方法个无穷序列去逼近精确解的方法.一般地一般地,有限有限步内得不到精确解步内得不到精确解.5 Gauss消元法是一种规则化的消元法，其基本消元法是一种规则化的消元法，其基本思想是通过逐次消元计算，把一般线性方程组的思想是通过逐次消元计算，把一般线性方程组的求解转化为等价的上三角形方程组的求解。求解转化为等价的上三角形方程组的求解。1 Gauss消去法消去法6 322413312321321321xxxxxxxxx消去后两个方程中的消去后两个方程中的x1得得 166225123232321xxxxxxx再消去最后一个方程的再消去最后一个方程的x2得得 613121321
4、xxx消元结束消元结束.5754222512332321xxxxxx经过回代得解经过回代得解:例例 考虑线性方程组考虑线性方程组 顺序顺序Gauss消去法消去法7 消元过程消元过程:先逐次消去变量先逐次消去变量 x1,x2,将方程组化将方程组化为同解的上三角形方程组为同解的上三角形方程组.上过方法可推广到一般情况上过方法可推广到一般情况 回代过程回代过程:按方程相反的顺序求解上三角形方程组按方程相反的顺序求解上三角形方程组.8 )1()1()1(3)1(2)1(1)1(3)1(3)1(33)1(32)1(31)1(2)1(2)1(23)1(22)1(21)1(1)1(1)1(13)1(12)1
5、(11)1()1(.),(nnnnnnnnnbaaaabaaaabaaaabaaaabA第一步第一步.设设 依次用依次用,0)1(11 a),.,3,2(,)1(11)1(11niaalii 乘矩阵的第乘矩阵的第1行加到第行加到第 i 行上行上,得到矩阵得到矩阵:iiijijbbaabb,AA )1()1(1)(1),则线性方程组的增广矩阵为则线性方程组的增广矩阵为记记9 )2()2()2(3)2(2)2(3)2(3)2(33)2(32)2(2)2(2)2(23)2(22)1(1)1(1)1(13)1(12)1(11)2()2(.0.0.0.),(nnnnnnnnbaaabaaabaaabaa
6、aabA其中其中njialaajiijij,.,3,2,)1(11)1()2(niblbbiii,.,3,2,)1(11)1()2(第二步第二步.设设 依次用依次用0)2(22 a),.,4,3(,)2(22)2(22niaalii 乘矩阵的第乘矩阵的第2行加到第行加到第i 行行,得到矩阵得到矩阵:10 )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.00.00.0.)(nnnnnnn(3)(3)baabaabaaabaaaab,A其中其中njialaajiijij,.,4,3,)2(22)2()3(niblbb
7、iii,.,4,3,)2(22)2()3(11 )()()3(3)3(3)3(33)2(2)2(2)2(23)2(22)1(1)1(1)1(13)1(12)1(11)()(.)(nnnnnnnnnnbabaabaaabaaaab,A这就完成了消元过程这就完成了消元过程.如此继续消元下去如此继续消元下去,第第n 1步结束后得到矩阵步结束后得到矩阵:12 .,.,.)()()2(2)2(22)2(22)1(1)1(12)1(121)1(11nnnnnnnnnnbxabxaxabxaxaxa 对此方程组进行回代，就可求出方程组的解对此方程组进行回代，就可求出方程组的解.对应的方程组变成：对应的方程组
8、变成：.1,2,1,)(,)(1)()()()(nniaxabxabxiiinijjiijiiinnnnnn13 能用顺序能用顺序Gauss消去法求解的条件是在消元过程中消去法求解的条件是在消元过程中得到的主元必须全不为得到的主元必须全不为0，即，即),2,1(,0)(nkakkk 顺序顺序Gauss消去法通常也简称为消去法通常也简称为Gauss消去法消去法.主元素都不为零主元素都不为零矩阵矩阵A的各阶顺序主子式都不为零的各阶顺序主子式都不为零.顺序顺序Gauss消去法中的消去法中的 称为称为主元素主元素.),.,2,1()(nkakkk)()2(22)1(11)det(kkkkaaaA 14
9、 顺序顺序Gauss消去法求解消去法求解n元线性方程组的乘除运算量元线性方程组的乘除运算量 第第1次消元乘除运算量次消元乘除运算量:消元过程乘除运算量消元过程乘除运算量求求 li1:(n1)求求aij(2):(n1)2求求bi(2):(n1)共共 (n21)次次 111)1(1222 nn15 回代过程乘除运算量：回代过程乘除运算量：求求 xn:1求求 xn1:2求求 x1:n.n 2116 nknkkk112)1()3(3123nnn n=30时时,顺序顺序Gauss消去法只需消去法只需9890次乘除法运算次乘除法运算.nnn 21111)1(1222 顺序顺序Gauss消去法求解消去法求解
10、n元线性方程组的乘除运算量元线性方程组的乘除运算量17 高斯消去法优缺点：高斯消去法优缺点：简单易行简单易行 要求主元均不为零要求主元均不为零,因而适用范围小因而适用范围小 数值稳定性差数值稳定性差 18例：例：单精度解方程组单精度解方程组 211021219xxxx 精确解精确解.,1000.00.1101191 x8个个.8999.99.0212 xx8个个 用顺序用顺序Gauss消去法计算：消去法计算：911212110/aal99921)2(2210101010.0.011 la8个个921)2(21012 lb 9991010011100,112 xx小主元小主元 可能导可能导致计算
11、失败致计算失败.2 主元素法主元素法19 若将方程组改写成若将方程组改写成:110221921xxxx用顺序用顺序Gauss消去法消去法,消元得消元得 12221xxx回代得解回代得解:x2=1,x1=1与准确解非常接近与准确解非常接近.可见可见,第一种算法是第一种算法是不稳定不稳定的的,第二种算法是第二种算法是稳定稳定的的.20 此例说明此例说明,在消元过程中在消元过程中,应避免选取绝对值应避免选取绝对值较小的数作主元较小的数作主元.如例中的第二种解法如例中的第二种解法,通过交换通过交换方程次序方程次序,选取选取绝对值较大的元素作为主元绝对值较大的元素作为主元.基基于这种想法导出了主元法于这
12、种想法导出了主元法.为了提高计算的数值稳定性为了提高计算的数值稳定性,在消元过程中采在消元过程中采用选择主元的方法用选择主元的方法.常采用的是常采用的是列主元消去法列主元消去法和和全主元消去法全主元消去法.21 给定线性方程组给定线性方程组 Ax=b,记记 A(1)=A,b(1)=b,列主元列主元Gauss消去法的具体过程如下消去法的具体过程如下:.1max)1(11)1(1行互换行互换行与第行与第为主元素，第为主元素，第 kaainik 首先在增广矩阵首先在增广矩阵 B(1)=(A(1),b(1)的第一列元素中的第一列元素中,取取 然后进行第一步消元得增广矩阵然后进行第一步消元得增广矩阵 B
13、(2)=(A(2),b(2).列主元消去法列主元消去法22.2max)2(22)2(2行互换行互换行与第行与第为主元素，第为主元素，第kaainik 然后进行第二步消元得增广矩阵然后进行第二步消元得增广矩阵 B(3)=(A(3),b(3).按此方法继续进行下去按此方法继续进行下去,经过经过 n 1步选主元和消元步选主元和消元运算运算,得到增广矩阵得到增广矩阵B(n)=(A(n),b(n).则方程组则方程组 A(n)x=b(n)是与原方程组等价的上三角是与原方程组等价的上三角形方程组形方程组,可进行回代求解可进行回代求解.只要只要|A|0,列主元列主元Gauss消去法就可顺利进行消去法就可顺利进
14、行.再在矩阵再在矩阵 B(2)=(A(2),b(2)的第二列元素中的第二列元素中,取取23 全主元素法全主元素法 每一步选绝对值最大的元素为主元素，保证每一步选绝对值最大的元素为主元素，保证 .1|ikl Step k:选取选取;0|max|)(,)(kijnjikkjiaakk If ik k then 交换第交换第 k 行与第行与第 ik 行行;If jk k then 交换第交换第 k 列与第列与第 jk 列列;消元消元列交换改变了列交换改变了 xi 的顺序，须记录的顺序，须记录交换次序交换次序，解完后再换回来。解完后再换回来。24例例 用主元素法求解线性方程组用主元素法求解线性方程组计
15、算过程保留三位小数计算过程保留三位小数,方程的精确解为方程的精确解为 x1*=1,x2*=2,x3*=3.1515613183312111321xxx25 1513181533126111 解解 1.按列主元素法，求解过程如下按列主元素法，求解过程如下 6111153312151318 167.5944.0167.105333.210151318 5333.210167.5944.0167.10151318 428.9142.300167.5944.0167.10151318 .000.1,000.2,001.3123xxx消元消元回代得回代得26 1513181533126111 解解 2.
16、按全主元素法，求解过程如下按全主元素法，求解过程如下 6111153312151318 167.5944.0167.105333.210151318 167.5167.1944.0051333.20153118 144.3572.10051333.20151318 .000.1,000.3,000.2132xxx回代得回代得3x2x27 全主元素法的精度优于列主元素法全主元素法的精度优于列主元素法,这是由于全这是由于全主元素是在全体系数中选主元主元素是在全体系数中选主元,故它对控制舍入误故它对控制舍入误差十分有效差十分有效.但全主元素法在计算过程中但全主元素法在计算过程中,需同时作行与列的需同
17、时作行与列的互换互换,因而因而程序比较复杂程序比较复杂,计算时间较长计算时间较长.列主元素法的精度虽然稍低于全主元素法列主元素法的精度虽然稍低于全主元素法,但其但其计算简单计算简单,工作量大为减少工作量大为减少,且计算经验与理论实且计算经验与理论实践均表明践均表明,它与全主元素法同样具有它与全主元素法同样具有良好的数值稳良好的数值稳定性定性.列主元素法是求解中小型稠密线性方程组的最好列主元素法是求解中小型稠密线性方程组的最好方法之一方法之一.例例3的计算结果表明的计算结果表明28 Gauss消元法的矩阵表示消元法的矩阵表示 3332312322211312111001001aaaaaaaaab
18、a 133312321131132312221121131211baabaabaaaaaaaaaaaaaa 333231232221131211aaaaaaaaa12arr 13brr 133312321131132312221121131211baabaabaaaaaaaaaaaaaa3 直接三角分解法直接三角分解法两者等价两者等价29 333231232221131211aaaaaaaaaA,112121aal 113131aal)2()2(33)2(32)2(23)2(22131211:00Aaaaaaaa jiijijalaa11)2(2,2 ji其中其中)2(1AAL,1001001
19、31211 llL两者等价两者等价 n=3时时Gauss消元法的矩阵表示消元法的矩阵表示30 )2(33)2(32)2(23)2(22131211)2(00aaaaaaaA)2(22)2(3232aal)3()3(33)2(23)2(22131211:000Aaaaaaa )2(2332)2(33)3(33alaa 其中其中)3()2(2AAL,10010001322 lL两者等价两者等价31ALLALA12)2(2)3()3(1211ALLA 1001001100010001100010001100100131213121llll:11 L 1001000110001000110001000
20、1100100013232ll:12 L 1010011001000110010013231213231211211llllllLL32LUaaaaaalllA )3(33)2(23)2(22131211323121000101001条件条件:.0,0)2(2211 aa,112121aal,113131aal.)2(22)2(3232aal:L:U矩阵矩阵A经经Gauss消元法后得到的上三角矩阵消元法后得到的上三角矩阵.33例例 求矩阵求矩阵 542774322A的三角分解的三角分解.解解 542774322A 860130322,22421 l,12231 l 600130322.2363
21、2 l 600130322121012001542774322LUA34 上述上述 n=3的情形可以推广到一般情形的情形可以推广到一般情形,例如例如 n=4 44434241343332312423222114131211aaaaaaaaaaaaaaaaA,112121aal )2(44)2(43)2(42)2(34)2(33)2(32)2(24)2(23)2(2214131211000aaaaaaaaaaaaa,113131aal,114141aal,)2(22)2(3232aal,)2(22)2(4242aal 35 )3(44)3(43)3(34)3(33)2(24)2(23)2(221
22、413121100000aaaaaaaaaaa,)3(33)3(4343aal )4(44)3(34)3(33)2(24)2(23)2(2214131211000000aaaaaaaaaa36LUaaaaaaaaaallllllA )4(44)3(34)3(33)2(24)2(23)2(22141312114342413231210000001010010001条件条件:.0,0,0)3(33)2(2211 aaa:U矩阵矩阵A经经Gauss消元法后得到的上三角矩阵消元法后得到的上三角矩阵.,112121aal,113131aal,)2(22)2(3232aal:L,114141aal,)2(
23、22)2(4242aal.)3(33)3(4343aal 37 20116126384102785124A 580012600030512424821 l14431 l341241 l 580012000030512423632 l03042 l 100012000030512442843 lLUA 1000120000305124140301210012000120116126384102785124例例38 )()3(3)3(33)2(2)2(23)2(2211312113213231211111)(nnnnnnnnnnnijaaaaaaaaaallllllLUaA条件条件:.0,0,0)
24、1(11)2(2211 nnnaaa,2,1111niaalii .2,1,)()(nknikaalkkkkikik :L:U矩阵矩阵A经经Gauss消元法后得到的上三角矩阵消元法后得到的上三角矩阵.Gauss消元法的矩阵表示消元法的矩阵表示Doolittle分解分解39 矩阵的三角分解矩阵的三角分解定理定理 设设 A 为为n 阶方阵阶方阵,若若 A 的前的前(n 1)阶顺阶顺序主子式序主子式 Ai (i=1,2,n 1)均不为均不为0,则矩阵则矩阵A存在唯一的存在唯一的Doolittle分解分解.)()2(2211iiiiaaaA 存在性：存在性：唯一性：唯一性：反证法反证法40 nnnnn
25、nnnuuuuuuuuuullllllA3332232211312113213231211111 Doolittle分解中分解中 LU 元素的求解次序元素的求解次序 Doolittle分解中分解中 LU 的另一求法的另一求法41 nnnnnnnnuuuuuuuuuullllllA3332232211312113213231211111njuajj 1,11U的第一行的第一行 1111ulaii1111ualii L的第一列的第一列 Doolittle分解中分解中 LU 元素的求解元素的求解 右端用矩阵乘法展开右端用矩阵乘法展开,比较两边的第比较两边的第1行和第行和第1列得列得42 假设已求得假
26、设已求得U 的前的前(r 1)行和行和 L的前的前(r 1)列列,r 1.下面求下面求U 的第的第 r 行和行和 L 的第的第 r 列列.右端用矩阵乘法右端用矩阵乘法,比较两边的第比较两边的第 r行的后行的后(nr+1)个元素个元素arj=L的第的第r行行向量与行行向量与U的第的第j列列向量的内积列列向量的内积(j r)0,0,(11jjrjjrjuuuu rjjrrrjrjrrjuululula 112211njrulululaujrrrjrjrrjrj ),(112211U的第的第r行行)0,0,1,(11 rrrll43.)1()1()(112211的内积的内积个分量个分量列的前列的前的
27、第的第个分量与个分量与行的前行的前的第的第 rjUrrLaulululaurjjrrrjrjrrjrjnjr U的第的第r行行44 右端用矩阵乘法右端用矩阵乘法,比较两边的第比较两边的第 r列的后列的后(nr)个元素个元素air=L的第的第i行行向量与行行向量与U的第的第r列列向量的内积列列向量的内积(i r)0,0,1,(111 iiirriillll)0,0,(11rrrrruuu rrirrrriririirulululula 112211 niruulululalrrrrriririirir ,)(112211L的第的第r列列 45rrirrrrrriririirirurrUriLau
28、ulululal 的内积的内积个分量个分量列前列前的第的第个分量与个分量与行前行前的第的第)1()1()(112211nir L的第的第r列列4611u1 irlrru1 jru1 1 rrl1ilirl1rlrjurruju1ru121l12u22u),2,1(nj),3,2(ni),1,(nrrj ),1(nri 对对r=2,3,n,jjau11 1111ualii)(112211jrrrjrjrrjrjulululau rrrrirririiriruulululal)(112211 矩阵三角分解的紧凑格式矩阵三角分解的紧凑格式47)2()2()3()4()7()7()4()2()5(例例
29、 用紧凑格式求矩阵用紧凑格式求矩阵 542774322A的三角分解的三角分解.解解223224 122 3227 1327 232)1(4 6123)1(5 LUA 60013032212101200148bLUxbAx 先求先求Ly=b,得得 y;再求再求 Ux=y,得得x.直接三角分解法直接三角分解法或或Doolittle分解法分解法.yDoolittle分解法分解法49 nnnnbbbyyy212121211.11lll Ly=b nkybybykiikikk,3,2,1111l 乘除运算量为乘除运算量为.2)1(nn50 1,2,1,1nniuxuyxuyxiinijjijiinnnn
30、 nnnnnnyyyxxxuuuuuu212122211211.Ux=y 乘除运算量为乘除运算量为.2)1(nn51例例 用直接三角分解法解用直接三角分解法解 223291076824312321xxx)2()1()3()4()2()8()7()6()10(解解21 3224 326 4)1(22 2328 14)1(37 32)1(3310 LUA 30024031211301200152 22329113121321yyy 9149321yyy 9149324312321321yyyxxx 321321xxx 先求先求Ly=b 得得 y 再求再求 Ux=y 得得x53 由于在求出由于在求出
31、uij,lij 和和 yi 后后,aij和和bi就无需保留了就无需保留了,故故上机计算时上机计算时,可把可把 L,U和和y存在存在A,b所占单元所占单元,回代时回代时x 取代取代 y,整个计算过程中不需要增加新的存储单元整个计算过程中不需要增加新的存储单元.在求一系列在求一系列系数矩阵相同而右端项不同系数矩阵相同而右端项不同的线性的线性方程组方程组 Ax=b(k),(k=1,2,m)时时(如求逆矩阵如求逆矩阵),用用三角分解法更为简便三角分解法更为简便.每解一个方程组每解一个方程组Ax=b(k)仅需仅需要增加要增加 n2 次乘除法运算次乘除法运算.解线性方程组解线性方程组Ax=b 的的 Doo
32、little三角分解法的计三角分解法的计算量约为算量约为n3/3,与与Gauss消去法相同消去法相同.54 Crout分解分解定理定理 设设 A 为为n 阶方阵阶方阵,若若 A 的前的前(n 1)阶顺阶顺序主子式序主子式 Ai (i=1,2,n 1)均不为均不为0,则矩阵则矩阵A可以唯一分解为可以唯一分解为其中其中 L 为下三角阵为下三角阵,U 为为单位单位上三角阵上三角阵.LUA 55 列主元的三角列主元的三角分解分解定理定理 设设 A 为为n 阶非奇异方阵阶非奇异方阵,则存在排列则存在排列矩阵矩阵P,使得使得其中其中 L 为单位下三角阵且为单位下三角阵且|lij|1,U 为上三角为上三角阵
33、阵.(在同样条件下在同样条件下,也可推出也可推出PA有有Crout分解分解)LUPA 定义定义(排列阵排列阵)单位阵经若干次两行互单位阵经若干次两行互换后而得到的矩阵换后而得到的矩阵.56 nnnnnnnnnddddxxxxbacbacbacb12112111122211dAx 特殊的稀疏矩阵特殊的稀疏矩阵 解三对角方程组的追赶法解三对角方程组的追赶法解三对角方程组解三对角方程组57 追赶法是求三对角线性方程组的三角分解法追赶法是求三对角线性方程组的三角分解法.追赶法本质上是追赶法本质上是Gauss消元法消元法.二阶常微分方程边值问题的差分离散方程组二阶常微分方程边值问题的差分离散方程组,热热
34、传导方程以及船体数学中建立的三次样条函数的三传导方程以及船体数学中建立的三次样条函数的三转角或三弯矩方程组均为三对角方程组转角或三弯矩方程组均为三对角方程组.58 4433322211bacbacbacb122bal 443332)2(2110bacbacbcb122)2(2clbb )2(233bal 443)3(32)2(21100bacbcbcb233)3(3clbb )3(344bal )4(43)3(32)2(211000bcbcbcb344)4(4clbb 例例 4阶三对角矩阵的三角分解阶三对角矩阵的三角分解59 )4(43)3(32)2(211000bcbcbcb 1111432
35、lll 4433322211bacbacbacb例例 4阶三对角矩阵的三角分解阶三对角矩阵的三角分解单位下二对角阵单位下二对角阵上二对角阵上二对角阵60定理定理 设三对角矩阵设三对角矩阵A满足下列条件满足下列条件 0|)12(0|,|0|11nniiiiiabnicacabcb则它可以分解成则它可以分解成A=LU nnnnnbacbacbacb11122211 nnnucucuculll132211321111 三对角矩阵的三角分解三对角矩阵的三角分解对角占优阵对角占优阵单位下二对角阵单位下二对角阵 上二对角阵上二对角阵工程中得到的工程中得到的三对角阵多数三对角阵多数满足此条件满足此条件61
36、nnnnnbacbacbacb11122211 nnnucucuculll132211321111 三对角矩阵三角分解中三对角矩阵三角分解中LU的求解次序的求解次序nniiullululu 1322162 nnnnnbacbacbacb11122211 nnnucucuculll13221132111111bu 122ual 1222clbu 对对 k=3,n 1 kkkual1 kkkkclbu 三对角矩阵三角分解中三对角矩阵三角分解中LU的求解的求解 右端用矩阵乘法展开右端用矩阵乘法展开,比较两边元素得比较两边元素得乘除运算量乘除运算量 2n263 解解三对角方程组的追赶法三对角方程组的追
37、赶法dLUxdAx 先求先求Ly=d,得得 y;再求再求 Ux=y,得得x.y追：消元过程追：消元过程赶：回代过程赶：回代过程64 Ly=d nkydydykkkk,3,2,111l 乘除运算量为乘除运算量为 n1 nnnddddyyyylll32132132111165 Ux=y 1,2,1,1nkuxyxuyxkkkkknnnc 乘除运算量为乘除运算量为 2n1 nnnnyyyyxxxxucucucu321321132211追赶法总的乘除追赶法总的乘除运算量运算量 5n466 追赶法的实质就是追赶法的实质就是Gauss消元法消元法,只是由于系只是由于系数中出现了大量的零数中出现了大量的零,
38、在计算过程中将它们撇开在计算过程中将它们撇开,从而使计算公式大大简化从而使计算公式大大简化,也大大减少了计算量也大大减少了计算量.为节省计算机存储单元为节省计算机存储单元,计算得到的计算得到的 lk,uk 分别分别存放在存放在 ak,bk 的存储单元内的存储单元内,而而 yk,xk 存放在存放在 dk 的的存储单元内存储单元内.当系数矩阵为满足定理条件的当系数矩阵为满足定理条件的严格对角占优阵严格对角占优阵时时,追赶法具有良好的数值稳定性追赶法具有良好的数值稳定性.674 平方根法与改进的平方根法平方根法与改进的平方根法求解对称正定方程组求解对称正定方程组,即即其中其中A为对称正定阵为对称正定
39、阵.bAx 利用对称性利用对称性,平方根法的计算量是平方根法的计算量是Gauss消去消去法的一半法的一半.68 对称正定阵的对称正定阵的Cholesky分解分解定理定理 若若 A 是对称正定阵是对称正定阵,则存在唯一的非奇异则存在唯一的非奇异下三角阵下三角阵L,使得使得且且 L的对角元素皆为正数的对角元素皆为正数,即即 lii 0 (i=1,2,n).TLLA 证明证明 A可进行可进行Doolittle分解分解.A为对称正定阵为对称正定阵,它的它的n个顺序主子式均大于个顺序主子式均大于0,设设ULA 其中其中 为单位下三角阵为单位下三角阵,U为上三角矩阵为上三角矩阵.L用反证法用反证法69 n
40、nuuuU2211 nnuuu2211 111UD 对对U进一步分解进一步分解.iiijuu A对称对称.TLU A正定正定对角矩阵对角矩阵D的对角元均为正的对角元均为正.令令),(diag2211nnuuuD 则则TTTLLLDDLUDLULA 其中其中 为非奇异下三角阵为非奇异下三角阵,且对角元均为正且对角元均为正.DLL A=LLT.iju70 333222312111333231222111333231222111llllllllllllaaaaaaA对称对称00 A为为3阶对称正定阵阶对称正定阵,A=L LT,怎样求怎样求L?23323223132223121311122222121
41、11211llllllllllllll对称对称1131311121211111 ,lallalal 322232213122222221 ,allllall 33233232231 alll:column1st:column2nd:column3rd22213132322212222/)(,lllallal 2322313333llal 71例例 对下列矩阵进行对下列矩阵进行Cholesky分解分解,241111 al,1112121 lal,3113131 lal,1111al,112121lal,113131lal,2212222lal ,2231213232lllal .23223133
42、33llal 22565172624A,41172212222 lal,24/8)(2221313232 lllal349222322313333 llal 323041002LTLLA 72 A为为n阶对称正定阵阶对称正定阵,A=L LT,怎样求怎样求L?21l11l22l1nl2nl2kl1klkklnklnnl L L中元素的求解次序中元素的求解次序 依次求依次求L的第一列的第一列,第二列第二列,第第n列列.73)0,0,0,0,(111ll )0,0,0,(22212lll )0,0,(3332313llll ),(321nnnnnnlllll.nnnnnnnnnnnnllllllll
43、llllaaaaaaA122121112122211121222111对称对称)(),(1jillllajkjkikjiij A为为n阶对称正定阵阶对称正定阵,A=L LT,怎样求怎样求L?li 为为 L的第的第 i 个行向量个行向量74 ),(),(),(),(),(),(21221211nnnnTllllllllllllLLA)(),(1jillllajkjkikjiij 对称对称 A为为n阶对称正定阵阶对称正定阵,A=L LT,怎样求怎样求L?752111111),(llla 11111),(llllaiii 求求 L 的第一列的第一列2222212222),(lllla 2222112
44、2),(llllllaiiii 求求 L 的第二列的第二列 A为为n阶对称正定阵阶对称正定阵,A=L LT,怎样求怎样求L?1111al),3,2(/1111nilalii 2212222lal ),3()(2221122nilllaliii 76 A为为n阶对称正定阵阶对称正定阵,A=L LT,怎样求怎样求L?2212221),(kkkkkkkkkklllllla kkikkkkikikikiiklllllllllla 112211),(设已经求得设已经求得 L 的前的前k1列列,现求现求L 的第的第 k 列列(k=3,4,n)(212221 kkkkkkkklllal,ki 对对kkkkk
45、ikikiikiklllllllal)(112211 Matlab函数：函数：chol 77 6023021912295631269 Acolumn1st112121/lal 431231 /l13341 /l1111al column2nd22222221all 1)2(5222 l3222322131allll 11/)2)(4 9(32 l4222422141allll 01/)2)1)(2(42 l例例 求正定阵求正定阵 的的Cholesky分解分解.3911 l23621 /l113131/lal 114141/lal 解解 L1423 011 78column4th44244243
46、242241allll 1)2()0()1(6 22244 lcolumn3rd33233232231alll 2)1()4(212233 l43433342324131allllll 22/)1)(0()4)(1(0(43 l例例 求正定阵求正定阵 的的Cholesky分解分解.解解 6023021912295631269 A L1423 011 22179bxLLbAxT 先求先求 Ly=b,得得 y,再求再求 LTx=y,得得 x.解正定线性方程组的解正定线性方程组的平方根法平方根法或或Cholesky分解法分解法.y平方根法平方根法或或Cholesky分解法分解法 设设A为对称正定阵为
47、对称正定阵80定理定理 若若 A 是对称正定阵是对称正定阵,则存在唯一的则存在唯一的单位单位下三下三角阵角阵L和对角阵和对角阵D,使得使得且且 D的对角元素皆为正数的对角元素皆为正数.对称正定阵的对称正定阵的LDLT分解分解TLDLA LUA nnuuuU2211iju nnuuu2211 111iiijuuUD 证明证明 A对称对称.TLU A正定正定对角矩阵对角矩阵D的对角元均为正的对角元均为正.81 对称正定阵的对称正定阵的LDLT分解本质上是对分解本质上是对A作作Doolittle分解分解,即即LU分解分解.LDLT分解中的分解中的 D=LU分解中的分解中的U的对角部分的对角部分LDL
48、T分解中的分解中的 L=LU分解中的分解中的L82 对称正定矩阵对称正定矩阵A的的LU分解分解,计算量可以计算量可以节省一半节省一半),2,1(11njaujj ),2(1111niuulii 求求U的第的第1行行 求求L的第的第1列列 对称正定阵的对称正定阵的LDLT分解中分解中L,D的计算的计算 先对对称正定阵先对对称正定阵A作作LU分解分解83),1(nkiuulkkkiik )()()1()1(2211njkulululaujkkkjkjkkjkj 求求U的第的第k行行 (k=2,3,n)求求L的第的第k列列 (k=2,3,n)对称正定阵的对称正定阵的LDLT分解中分解中L,D的计算的
49、计算节省了计算量节省了计算量 nnuuuD2211 求求D84例例 求矩阵求矩阵 222322191631522624A的的 LDLT分解分解.解解42 642 46415 231 91919 2 42 413 3232 1614122 )4()2()6()5()1()19()2()3()2()22(4244 9385 131121012123001210001L 16944DTLDLA 86 解正定线性方程组的解正定线性方程组的改进改进平方根法平方根法或或LDLT分解法分解法.bxLDLbAxT 先求先求 Ly=b,得得 y,再求再求 LTx=D 1y,得得 x.y改进改进平方根法平方根法或
50、或LDLT分解法分解法 设设A为对称正定阵为对称正定阵87 平方根法与改进的平方根法的优点平方根法与改进的平方根法的优点 计算无须选主元计算无须选主元,由于正定性由于正定性,计算过程是计算过程是数数值稳定值稳定的的 计算量是计算量是Gauss消元法的一半消元法的一半 由于对称性由于对称性,实际计算可存储一半实际计算可存储一半 是求解中小型稠密正定线性方程组的好算法是求解中小型稠密正定线性方程组的好算法885 误差分析误差分析 用直接法解线性方程组，初始数据会有误差，用直接法解线性方程组，初始数据会有误差，计算过程同样会产生误差，这就需要对这些误计算过程同样会产生误差，这就需要对这些误差作一些分
