线性方程组的数值解法LU分解法

3.5 LU分解法n我们知道对矩阵进行一次初等变换，就相当于用相应的初等矩阵去左乘原来的矩阵。因此我们从这个观点来考察Gauss消元法并用矩阵乘法来表示，即可得到求解线性方程组的另一种直接法：矩阵的三角分解。高斯消元过程的矩阵表示)2()1(1)2()2(2)2()2()1()1()1()1()1()1()1()1()1()1()1()1(131211)1(11)1(1111131121)1(11.1.00.1011.,32i)()(1),.,0:122211211212222111211AALaaaaaaaaaaaaaaaalllni-laalaaaannnnnnnnnnnniiin，其矩阵形式为，行行令消零时，将步等价于第则）（）（）（高斯消元过程的矩阵表示)3()3()3(3)3(3)3(33)2(2)2(23)2(22)1(1)1(13)1(12)1(11222)2(22)2(222322)2(22.00.00.0.:),.,4,3(1.00.0.100.0100.00102AaaaaaaaaaaaALAniaalllLannnnnn)()(iin，即有左乘时，用矩阵步等价于：若同理第高斯消元过程的矩阵表示1.1111.11.000.00.0.2321212111)()3(3)3(33)2(2)2(23)2(22)1(1)1(13)1(12)1(111221nnnnnnnnnnllLllLUaaaaaaaaaaALLLL因为以此类推可得高斯消元过程的矩阵表示为上三角阵为单位下三角阵，其中所以U1.111.).(1213231211112121111221LLUUllllllULLLLULLLLAnnnnnnnnLU分解法()U.L UALUxybALULyxx bb由此解线性方程组就等价于解两个三角方程：因此，关键问题在于能否对矩阵 直接进行 分解.其中 为单位下三角阵，=为上三角阵A矩阵分解理论),.2,1(0.1111nkaaaaAkkkkk矩阵分解理论 矩阵分解理论 3.2.2 Doolittle分解1112131112132122232122233132333132331111212111 212111313111 313111,31111(1,2,3)jjjjL Unaaauuuaaaluuaaallukauuajaau lluaau llu此分解在于如何算出的各元素，以为例时：由得；由得)(322332133133333323321331332212313232233212313213212323231321231221222222122122ululauuululakuulalululaulauuulaulauuulak得时：由得由；得由；得时：nA的各阶顺序主子式均不为零，即),.2,1(0.1111nkaaaaAkkkkkDoolittle分解11121111212122221222n1n2n12.1.1.1,nnnnnnnnnaaauuuaaaluuaaalluL U令用比较等式两边元素的方法逐行逐列求解各元素Doolittle分解。得再由；得由时：。得再由；得由时),.,4,3(),.,3,2(12),.,3,2(),.2,1(1:12212122222121212222121211111111 i1111niuulalululanjulauuulakniualluanjauuakiiiiiijijjjjjiiijjjjDoolittle分解1111211211,.1,000.10.,.,kttjktkjkjktkjtjktjjjjkkkkkjknkkkknkkjulauuuluuulllakjuuuk）（有步时：计算第Doolittle分解11111111,.1/)(000.0,1,.,.,ktkktkitikikktkkiktkitkkkikiiknkkknkiuulalululuullakill）（得，于是由由于计算Doolittle分解nnnnnnkkkttkitikikkttjktkjkjulluuluuunkuulalnkinkjulau.A,.,2,1/)(),.,1;,.,(2122221112111111的各位各元素在计算机内存于即Doolittle分解。可获解得及再由TniinijjijiiijjijiixxxnniuxuyxniylbyULxyxby),.,(1,.1,/)(,.,2,121111例题30191063619134652.D.1321xxxoolittle分解求解方程组试用例例题。，；，令、分解解：326246521101001636191346521311121131211332322131211323121luluuukuuuuuulllALU例题LUAululaukuulalulauulauk473652143121434/)(7)6(21935213223321331333322123132321321232312212222所以时：，时：例题TyyyyyyybLy)4,1,10(43034,12019,1030191014312112321321即得）解（、解方程组例题。所以方程组的解为解得：解TxxxxxxxyUx)1,2,3(3,2,14110473652)2(123321Doolittle分解).,.,()(),.,(/)(;/)2),.,)(;)1)(),.,(/)(),.,(,.,)2),.,(/)(),.,(),.,()1(nixniuxuyxayxniylbybyyUxbLynkiuulalankkjulauankniaalaLUAnibnjiaiiinijjijiinmnnijjijiikkkttkitikikikkttjktkjkjkjiiiiij2141212311323212212111111111111111输出：方程组的解；和解；做对；分解输入：3.6 平方根法(Cholesky分解法)n在应用数学中，线性方程组大多数的系数矩阵为对称正定这一性质，因此利用对称正定矩阵的三角分解式求解对称正定方程组的一种有效方法，且分解过程无需选主元，有良好的数值稳定性。对称矩阵的Cholesky分解n A对称：AT=A A正定：A的各阶顺序主子式均大于零。即 ),.2,1(0.1111nkaaaaAkkkkk对称矩阵的Cholesky分解n由Doolittle分解，A有唯一分解 。，也就是所以，有即TTTTTTTTTLLAULULLULULULULUALU)(A对称矩阵的Cholesky分解n定理 设A为对称正定矩阵，则存在唯一分解A=LDLT，其中L为单位下三角阵，D=diag(d1,d2,dn)且di0(i=1,n)对称矩阵的Cholesky分解n证明：ndddUULLUADoolittle21111,A令非奇异的上三角阵。为为单位下三角阵，其中分解为分解可唯一是对称正定，由因为对称矩阵的Cholesky分解11121221212A0000.;.00,.,nnnnAddd ddd dddd dd 又因为 是正定的，则 的顺序主子式均大于零故有得；由得；由得。即均大于零.对称矩阵的Cholesky分解。所以为单位上三角阵。为对角阵其中故LDUAUDDUdddddddUn,1.1*.*1*.*12211211对称矩阵的Cholesky分解。，所以故有对称，即又因为TTTTTTLDLAULADLULDUAA)(推论：设A为对称正定矩阵，则存在唯一分解 其中L为具有主对角元素为正数的下三角矩阵。TLLA 对称矩阵的Cholesky分解n证明：0,.,0,0)(4.2.3),.,(2121212121nTTTndddLLLLDLDLDLAddddiagD的主对角元素为其中则由定理令Cholesky分解的求法332322131211333231222111333231232221131211212221113?.,llllllllllllaaaaaaaaanlllllllLLLAAijnnnnT为例。以如何求令则对称正定设Cholesky分解的求法。，得由；，得时：由。；同理得，得由；，得时：由2221313232223221313222122222222212211313111212111212111112111121lllalllllalalllaklallalllaallakCholesky分解的求法njnjilllallalnlallllakjjjkjkikijijjkjkjjjjii,.,2,1,.,1/)()(,311211122123333323323223133有阶行列式推广到，得时：由Cholesky分解法TTLLAyxLbLybAXcholesky其中分解法解线性方程组用Cholesky分解法缺点及优点 优点：可以减少存储单元。缺点：存在开方运算，可能会出现根号下负数。改进Cholesky分解法n改进的cholesky分解A=LDLTnnnnnnnnTnnnnnndlddldlddldldlddllllllDLLAdddDllllllL.1.111)(1.11133322322211311211112132312121121323121由改进的cholesky分解1121111211,.,2,1)1,.,2,1(/)(),.,2,1()1,.,2,1(jkkikiiijjkjkkikijijjkikikiijijjkjkkikijnidladijdldlalniddlaijdlldlaji由此可得有逐行相乘，并注意到改进的cholesky分解)1,.,2,1,.,3,2(,1111ijnilcaddcllcacdcldlcikikikiiijijijjkjkikijijjijijjijij，成所以可将上述公式改写，则可令为减少计算量TnnnTdydydyyDdddDDLLACholeskyyxbybx),.,(111221112111故即等价于求分解法解线性方程组用改进的cholesky分解算法；做对做对分解，输入1111)2(1,.,2,1)1(,.,3,2.2);,.2,1(),.,2,1(.1ikikikiiiiijijijijjkikikijijTiijlcadadclalcacijniLDLAnibni,ja改进的cholesky分解算法。输出：；解；解解),.,2,1()3()1,.,1()2(),.,2()1(.3111111nixnixldyxdyxyDxLniylbybybLyAinikkkiiiinnnTikkikiibx例题123411512231236CholeskeyxxxADoolitte 例试用改进的分解算法解方程组解：对系数矩阵 做分解例题TLDLA11125.025.0111.7541125.025.01175.11.751141125.025.011125.075.175.125.01143125.075.175.125.01143225.02225.0114321221114例题TybyyyyyyyL)3,47,5(3,75.1,5635110.25125.01321321即得得由例题111233215(,1,3)410.250.251.25111133,2,1(1,2,3)TTTDLDxxxxxxyxyx而，由得得所以方程的解：nA=LDLT分解，既适合于解对称正定方程组，也适合求解A为对称，而各阶顺序主子式不为零的方程组n而对A=LLT只适合于对称正定方程组3.7 三对角方程组求解的追赶法11222111|i(1),|0ijn nnnnijnnAaacbacAbacbjaa设三对角矩阵，对一切有即三对角方程组求解的追赶法11112222223111111111,1111(2,3,.,)nnnnnnnnniiiiii iADoolittleLUacqcbacpqcpbacqcbapqqabpinqqapc 如果 满足分解则有A=其中三对角方程组求解的追赶法),.,2(1111),.,(1113213213221niypfyfyffffyyyypppfffULAiiiinnnTnfyxfyfx解得，故有其中等价于求则求三对角方程组求解的追赶法组的追赶法。以上称为解三对角方程解得再由)1,.,1(1121121112211niqxcyxqyxyyyyxxxxqcqcqcqiiiiinnnnnnnnnn三对角方程组求解的追赶法n其计算工作量为5n-4次乘除法。工作量小，其实现的条件为qi不为零。有以下定理可得证三对角矩阵求解的充分性条件。且追赶法可以实现。非奇异则矩阵，且满足形如设三对角矩阵定理,|,|)3()1,.,3,2(|)2(),.2,1(0)1(,11).2.3(2.2.311AabcbnicabnicAnniiii解三对角矩阵线性方程组的追赶法程序框图补充1:矩阵求逆-分块乘法 1112121212,(,)(,)(,)(,)1,2,nnnnjjAAEEAxxxEe eeA xxxe eeAxejn-1将A按列分块则补充2:初等变换法求矩阵的逆12212|GaussJordanA EXXAnnnAAn 初等变换由顺序消去法知，求解上式即为（）（E）则。算法上增加了个存储单元即可实现，但当 很大时，这个单元还是很困难的。如果将的单元，最后仍用来存放则可节省个单元，这就是原地求逆。补充3:矩阵原地求逆法11()1()()()()()()111:.1.()1()1.nnkkkikkkkkkKikkkkkkknknnGaussJordanL LL AElaikaLllikalAL LL实现 由消去法矩阵形式其中故n为使求逆过程不断提高求解精度，因此增加选主元工作，最常用的是选列主元求逆。因此增加一个数组Z(n)，记录选主元的交换号，最后在消元工作完成后，根据Z(n)对A中的元素进行相应的列交换，得到A-1GaussJordan原地求逆法算法（原地求逆法）);,.2,1()3(;,0)2(;)(|max|)1(),.2,1.3);,.,2,1()(.2),.2,1,(.1njaathenklifDifaDilikzaankniiiznjialjkjlkkkiknikkiijk停机输入奇异标志；选主元：做对；输入：结束。输出：；做对),.2,1,(.5);,.,2,1(thenif)(1,.,1.4);,.,2,1()7();,.,2,1,()6();,.,2,1()5(;1)4(njianiaaktkZtnkkiniacakjkinjiaaaakjnjacaacaijitikikkikkjikijijkjkkjkkkkk例题1122111221AA例设矩阵求逆矩阵。102011001520310221100010001122111221消去法解法一：JordanGauss例题原地求逆法解法二：所以得JordanGaussA1203514611203514611000100011200110211003104011例题11251212102121121251121012112122111121121221111122213)1(1)1(AAczk例题212210212512131121021251212112113)2(2)2(AAczk，例题32320151361420125133142012512131123)3(3)3(AAczk，例题1203514611203514610211531643)1(3)2(13AAzz所以，由与第三列交换第一列与第三列交换第二列