清华第五版数值分析第5章课件.ppt

1、 第五章第五章 线性代数方程组的数值解法线性代数方程组的数值解法 ,0,X,.11ijn nTAXbAAT nna,),)(x(bxbb 矩矩阵阵表表示示记记为为这这里里我我们们假假设设11112211211222221122 nnnnnnnnnnna xa xa xba xa xa xba xa xa xb 阶阶线线性性方方程程组组：引言引言 关于线性方程组的数值解法一般有两类。关于线性方程组的数值解法一般有两类。直接法直接法：经过有限步算术运算，可求得方程组：经过有限步算术运算，可求得方程组的精确解的方法（若在计算过程中没有舍入误的精确解的方法（若在计算过程中没有舍入误差）差）.代表性的算

2、法是高斯代表性的算法是高斯(Gauss)(Gauss)消去法。消去法。计计算代价高算代价高.迭代法：迭代法：用某种极限过程去逐步逼近线性方程用某种极限过程去逐步逼近线性方程组精确解的方法组精确解的方法.简单实用。简单实用。思思路路首先将首先将A通过通过初等行变换初等行变换化为化为上三角阵上三角阵，此此过程称为过程称为消去过程消去过程，再求解如下形状的方，再求解如下形状的方程组，此过程称为程组，此过程称为回代求解回代求解。5.1 高斯消去法高斯消去法=,.,n,iubxu.bxu.xubxu.xuxuiinnnnnnnn2102222211212111 其中，其中，12,.11,ni)/uxu(

3、bx/ubxiinijjijiinnnn 上三角方程组的回代上三角方程组的回代求解求解：二二 高斯消去法步骤如下：高斯消去法步骤如下：第一步：第一步：niiaaai,2,1,0)1()1()1(11111 行行第第行行第第当当 )1()1()1()1()1()1(2)1()1()1()1()1()1(2122221111211nnnnbaaabaaabaaannn )2()2()2(2)2(2)2(2)2(22)1()1()1()1(00111211nnnnnbaabaabaaan),.,2(),.,2;,.,2(1)1(1)1()2()1(11)1()2(nimbbbnjniamaaiiij

4、iijij niaamii,.,3,2,)1(11)1(11 令令则则设设 Ax=b.记记A(1)=A b(1)=b第二步：第二步：niiaaai,3,2,0)2(22)2(2)2(22 行行第第行行第第当当 )(n)(nn)(n)()(n)()()(n)()(nbaabaabaaabaaaa333333333332222223222(1)1(1)1(1)13(1)12(1)1100000 )(n)(nn)(n)()(n)(nbaabaabaaa22222222222(1)1(1)1(1)12(1)1100)3()33(2222322223n,.,ibmbbn,.,j;n,.,iamaa)(i

5、)(i)(i)(ji)(ij)(ij niaamii,.,3,2,)2(22)2(22 令令则则 重复上述过程重复上述过程,最后得最后得 )n(n)k(k)()()n(nn)k(kn)k(kk)(n)()(n)()()n()n(b.b.bba.a.a.a.aa.aabA22112222211112111)n()n(bx A再再解解原方程组的同解方程组原方程组的同解方程组 n,.,ki,aam)k(kk)k(ikikn,.,ki,bmbbn,.,kj;n,.,ki,amaan,.,k)k(kik)k(i)k(i)k(kjik)k(ij)k(ij111112111 其中其中对于对于Gauss Ga

6、uss 消去法消去法(1)(1)(1)(1)11112211(2)(2)(2)22222()()()01 2nnnnnnnnnniiia xa x.a xbax.axb.axba,i,.,n其中，()()()()()11,.2 1nnnnnnniiiiiijjiij ixb/ain,x(ba x)/a 上三角方程组的回代上三角方程组的回代求解求解：定理：定理：若若A的各阶的各阶顺序主子式不等于顺序主子式不等于0，则高斯消，则高斯消去法能顺序进行消元，得到唯一解。去法能顺序进行消元，得到唯一解。iiiiiaaaaA.)det(1111 高斯消去法条件(高斯消去法适用范围)证明：高斯消去法能进行下

7、去的充要条件是而消元过程不会改变顺序主子式的值。所以解的唯一性又要求(1)(2)()1122det()iiiiAa aaL()0iiia()1det()/det()iiiiiaAA(1,1)inL|0A 高斯消去法存在的问题：高斯消去法存在的问题：2.2.如果某个如果某个 ,但很小，会引入较大的误差。但很小，会引入较大的误差。1.1.在高斯消去法消去过程中可能出现在高斯消去法消去过程中可能出现 的情况，这时高斯消去法将无法进行的情况，这时高斯消去法将无法进行()0kkka()0kkka：消消元元因因子子)(k(kkkikikaam 三三 高斯列主元消去法高斯列主元消去法 nnnnnnn)()(

8、baaabaaabaaab,Ab,AbAX2122222111121111 的增广矩阵为的增广矩阵为设方程组设方程组具体步骤为：具体步骤为：然后进行消元，得然后进行消元，得行，行，行和第行和第交换第交换第第一步：选第一步：选1111,1,max1iaainii 基本思想基本思想：在每轮消元之前，选列主元素：在每轮消元之前，选列主元素（绝对值最大的元素）（绝对值最大的元素）)(n)(nn)(n)()(n)()()(n)()()()(baabaabaaab,A22222222222111111211122.,2,)3()3(2)2(22)2(2,max2bAiaainii然然后后进进行行消消元元，

9、得得行行，行行和和第第交交换换第第第第二二步步：选选 )2()2()2(2)2(2)2(2)2(2211121100nnnnnnbaabaabaaa以第二步为例以第二步为例:挑选从第二行开始挑选从第二行开始的第二列中的最大的第二列中的最大元素元素,交换行将其变交换行将其变为主元为主元.kik,aakk)k(knnik)k(k,imaxk次消元次消元然后进行第然后进行第行，行，行和第行和第交换第交换第步：选步：选第第 如此至多经过如此至多经过n-1n-1步步,就得到与之同解的上三角形方就得到与之同解的上三角形方 程组的增广矩阵程组的增广矩阵,再用回代过程即可得方程组的解再用回代过程即可得方程组的

10、解.12324 6349 2511 34xxx 例：用Gauss列主元消去法解方程组24 6 349 2 549 2 524 6 311 3 411 3 4 解：49 2511052255 110424 49 2555 110424110522 492555110424300 45 12rr2112rr3114rr23rr3225rr12313953,20220 xxx 回代求解得：高斯列主元消去法的适用范围 定理：系数矩阵非奇异，这高斯列主元消去法可行。高斯消元法每一消元步骤相当于对系数矩阵左乘一个行初等矩阵，具体为nnI 111nnmE 111212nnmE 1111313第一行乘第一行乘

11、 m21加到第二行加到第二行nnI 111第一行乘第一行乘 m31加到第三行加到第三行nnnnmE 1111nnI 111第一行乘第一行乘 mn1加到第加到第n行行记记 L1=EnE3E2 高斯消元法第一步：高斯消元法第一步：n,irmrii32)(11 )1()1()1()1()1()1(2)1()1()1()1()1()1(2122221111211nnnnbaaabaaabaaannn )2()2()2(2)2(2)2(2)2(22)1()1()1()1(00111211nnnnnbaabaabaaan等价于等价于L1A(1)=A(2),L1b(1)=b(2)A(1)b(1)A(2)b(

12、2)1001001121nmm)()(iiaam111111 令令高斯消元法的矩阵形式高斯消元法的矩阵形式 )3()3()3(3)3(3)3(3)3(33)2(2)2(2)2(23)2(221113121100000nnnnnnnbaabaabaaabaaaa )2()2()2(2)2(2)2(2)2(22)1()1()1()1(00111211nnnnnbaabaabaaanA(2)b(2)A(3)b(3)n,irmrii43)(22 )()(iiaam222222 令令第二步消元第二步消元等价于等价于L2A(2)=A(3),L2b(2)=b(3)100010001000012322.m.m

13、.Ln其中其中Lk A(k)=A(k+1),Lkb(k)=b(k+1)一般第一般第k步消元的矩阵表示为步消元的矩阵表示为Lk=1.11,1knkkmm 其中其中 1kL,1.11,1knkkmm Ua.a.aa.aaa.aaaALL.LL)n(nn)(n)()(n)()()(n)()()(nn 0000003333322223222111131121111221最最后后可可得得A(n)为上三角阵为上三角阵为单位下三角阵，为单位下三角阵，其中其中所以所以U1.111.).(1213231211112121111221LLUUmmmmmmULLLLULLLLAnnnnnnnn L 为单位下三角阵而

14、为单位下三角阵而 U 为为一般一般上三角阵的分解称上三角阵的分解称为为Doolittle 分解分解（2）L 为一般下三角阵而为一般下三角阵而 U 为为单位单位上三角阵的分解称为上三角阵的分解称为Crout 分解分解。基于这些想法我们得到下面的三角分解解线性方程组的方法基于这些想法我们得到下面的三角分解解线性方程组的方法 L5.2 矩阵的三角分解法矩阵的三角分解法一一 矩矩阵的三角分解阵的三角分解 若矩阵若矩阵A有分解：有分解：A=LU，其中，其中L为下三角阵，为下三角阵，U为上三角阵，则称该分解为为上三角阵，则称该分解为A的的LU分解。分解。若矩阵若矩阵A有分解有分解A=LU，则解线性方程组，

15、则解线性方程组 就等价于求解就等价于求解Ax b()LybL UUxyxb当系数矩阵当系数矩阵A相同，有很多不同的右端向量时，相同，有很多不同的右端向量时，用此方法很有效。如求用此方法很有效。如求即求方程即求方程的解，用此法很合适。的解，用此法很合适。下面看什么情况下，可对系数矩阵做三角分解。下面看什么情况下，可对系数矩阵做三角分解。-1A(1,2,)iAxe inL三角分解的存在唯一性定理 定理。分解，且分解是唯一的分解，且分解是唯一的必可作必可作则方阵则方阵零，即零，即的各阶顺序主子式不为的各阶顺序主子式不为阶方阵阶方阵如果如果LUAAaaaaa,0,0,0An2221121111 引理：

16、若 则()ijALUUu1122det()(1,)iiiAu uu inLL存在唯一性证明证明：存在性：高斯消去法能进行下去，则三角分解存在。唯一性：设则下三角阵=上三角阵，所以所以1122ALUL U111212L LU U111212L LU UI1212LLUU,4L Un 此此分分解解在在于于如如何何算算出出的的各各元元素素，以以为为例例 000000 1010010001 4434334423221413121143424132312144434241343332312423222114131211uuuuuuuuuullllllaaaaaaaaaaaaaaaauululululul

17、ululululuululuulululululuuluuluululuuuu4434432442144133432342134122421241114134243214313323321331223212311131241421231321221221112114131211 二二 矩矩阵的直接三角分解阵的直接三角分解直接计算直接计算 A 的的 LU 分解分解(例例)（续）（续）uulululululululululuululuulululululuuluuluululuuuu4434432442144133432342134122421241114134243214313323321331

18、223212311131241421231321221221112114131211 ;,;,114141113131112121114141313121211111ualualualauauauaulu 列列行行ninjualauiijj,2,1,111111 ;)-(,)-(;-,-,-2212414242221231323221421242413212323122122222uulaluulalulauulauulaulu 列列行行.,3,)-(,2,-221222212122n injuulalulauiiijjj 1)1)2)2)34432442144144444ulululauu

19、行行.,4 )-(,3,-332321313313133niulnjuluulalulauiiii2j32jjj ；列列行行 3)/uulul(a l lulululauululauu3323421341434334332432143134342332133133333;,3)3)一般计算公式一般计算公式n ,2,i ,n,1,j ,111111 ualauiijj 计算量与计算量与 Gauss Gauss 消去法同消去法同.1111 2,3,()/1,kkjkjkqqjqkikikiqqkkkqknjkniknual ulal uu 对对计计算算(1)LU计计算算矩矩阵阵，；(2)11-1

20、2,3,1 Lybiiniijijjybyybl 解解方方程程组组：LU 分解求解线性方程组分解求解线性方程组 (1,1)13 nnnn innijjiiiij iUxyxyuxyu xu （）解解方方程程组组:a11 a12 a13 a1n u11 u12 u13 u1n a21 a22 a23 a2n l21 u22 u23 u2n a31 a32 a33 a3n l31 l32 u33 u3n an1 an2 an3 ann ln1 ln2 ln3 unn.(1)(3)(5)(2n-1)(2)(4)(6)(2n)紧凑格式紧凑格式旧元素减去左边行与旧元素减去左边行与顶上列向量的点积顶上列向

21、量的点积计算行不用除法计算行不用除法计算列要除主对角元计算列要除主对角元注注：111111-11-11,ijjikkjkjkqqjqkikiqqkqikkkaualuualul ual u 例例 求矩阵的求矩阵的Doolittle分解分解 11242142612332442A 1122214161232/32442 112221013632/32442 1519220501363232442/111111-11-11,ijjikkjkjkqqjqkikiqqkqikkkaualuualul ual u 9000050036302442U,/L,LUA 1519220101001230001其中

22、其中 9519220501363232442/12325610413191963630 xxx 例例:用矩阵的直接三角分解法解方程组用矩阵的直接三角分解法解方程组或或 用用 Doolittle 分解法分解法111213212223313233111213312121311111(1)256100413191063611256462322ALUuuuluullukuuuaalluu 解解：分分解解，令令，；，。时时，LUAululaukuulalulauulauk 473652143121434/)(7)6(21935213223321331333322123132321321232312212

23、222所以所以时：时：，时：时：123123211021193413010,19201,34304(10,1,4)TLybyyyyyyy （）解解得得即即123321(3)25610371441,2,3(3,2,1)TUxyxxxxxxx 解解解解得得：所所以以方方程程组组的的解解为为。三三 解三对角线性方程组的解三对角线性方程组的11112222211111nnnnnnnnnbcxdabcxdAabcxdabxd 可以证明可以证明A的的Crout分解分解形式为：形式为：,111c ,222c ,2212ba ,ckkk ,1kkkkba ,11b LU11111,/bcb 2221222,/

24、bac 1,(2,3,)/,(2,3,1)kkkkkkkbaknckn 11111211122111122211nnnnnnnnnnaaabacbacbacb111111,/,2,3,/,2,3,1iiiiiiibcbbaincin 1111/()/,2,3,iiiiiydbyda yin 1,1,2,1nniiiixyxyxinn 追的过程追的过程赶的过程赶的过程 定理定理 设设A为对称正定矩阵，则存在唯一分解为对称正定矩阵，则存在唯一分解A=LDLT，其中，其中L为单位下三角阵，为单位下三角阵，D=diag(d1,d2,dn)且且di0(i=1,n)四四 系数矩阵为对称正定矩阵方程组的平方

25、根法系数矩阵为对称正定矩阵方程组的平方根法 A对称：对称：AT=A A正定：正定：A的各阶顺序主子式均大于零。的各阶顺序主子式均大于零。即即 ),.,2,1(0.1111nkaaaaAkkkkk 证明：证明：ndddUU,LLUAA,Doolittle21A令令上三角阵。上三角阵。为为为单位下三角阵为单位下三角阵，其中，其中分解为分解为可唯一可唯一分解分解是对称正定，由是对称正定，由因为因为均均大大于于零零即即。得得由由；得得；由由得得故故有有的的顺顺序序主主子子式式均均大大于于零零是是正正定定的的，则则又又因因为为nnnnddddddddddddA,.,00.;0000A2121221211

26、1 。所所以以为为单单位位上上三三角角阵阵。为为对对角角阵阵其其中中故故0002211211111LDUAU,DDU.d*.d*d*.d*dddUn DU0 即为即为L、U分解中的分解中的UDU0。，所以，所以故有故有对称，即对称，即又因为又因为TTTTTTLDLAULADLU)LDU(AA 000推论：推论：设设A为对称正定矩阵，则存在唯一分解为对称正定矩阵，则存在唯一分解 其中其中L为具有主对角元素为正数的下三角矩为具有主对角元素为正数的下三角矩阵。阵。TLLA Cholesky分解U从而：从而：其中其中 为下三角阵。为下三角阵。12nddDd 111/21/222nnddddDDdd 1

27、/21/21/21/211)(TTTTALDLLD D LLDLDLL （）1/21LLDCholesky分解的求法分解的求法 nnnnTl.ll.lllLLLA,A21222111令令则则对称正定对称正定设设121/2111(),(1,2,)1,2,()()1,iiiiiikkjijijikjkjjklalinjnlal llijn 计算公式计算公式:求解对称正定方程组求解对称正定方程组Ax=b的的平方根法平方根法(计算公式计算公式)：1 1分解计算分解计算 （对（对L L按列计算）按列计算）：TLLA 21112)()1(/ikikiiiilal ,l)lla(ljjjkjkikijij

28、11)2()121(n,ji;n,j 求求解解计计算算.2，),2,1()()3(11nilylbyiiikkikii 。)1,2,1,()()4(1 nnilxlyxiinikkkiii xyxLybLyT求求，求求，),2,1(ni 平方根法平方根法缺点及优点缺点及优点 优点：可以减少存储单元。优点：可以减少存储单元。缺点：存在开方运算，可能会出现根号下缺点：存在开方运算，可能会出现根号下负数。负数。改进平方根法改进平方根法 分解分解A=LDLT1212313212111 211 3111212232223132333121111.1()1.1.1.1,nnnnnnTnnnnnnndldl

29、lLDllldAL DLdd ld ld lldd ld llldd llll 由由nd ,.TTijii jiALULDLUDLud l改进平方根法11(,1,)(1,2,)kkjkjkppjpkiikkkual ujk knulikknu LL令令DLTx=y，1)解下三角形方程组解下三角形方程组Ly=b得得)n,i(ylbykikikii2111 2)解上三角形方程组解上三角形方程组DLTx=y得得 1(,1,2,1)niiki kk iiiyxl xin nu L123411512231236xxxADoolitte 例例试试用用改改进进的的平平方方根根法法解解方方程程组组解解：对对系

30、系数数矩矩阵阵做做分分解解4114114111220.25220.251.751.751230.25230.251341114110.251.751.750.2511.751.750.25110.25111A =LU11(,1,)(1,2,)LLkkjkjkppjpkiikkkual ujk knulikknu 123123150.25130.251165,1.75,37(5,3)4TLyyyyyyyby 由由得得得得即即123321541171.751.754133,2,1(1,2,3)Txxxxxxx 得得所所以以方方程程的的解解：求解求解 时，时，A 和和 的误差对解的误差对解 有何影响

31、？有何影响？bxA bx 设设 A 精确，精确，有误差有误差 ，得到的解为，得到的解为 ，即，即bb xx bbxxA )(bAx 1|1bAx 绝对误差放大因子绝对误差放大因子|xAxAb 又又|1bAx|1bbAAxx 相对误差放大因子相对误差放大因子5 线性方程组的性态和解的误差分析线性方程组的性态和解的误差分析6 Error Analysis for .bxA 设设 精确，精确，A有误差有误差 ，得到的解为，得到的解为 ，即，即bA xx bxxAA )(bxxAxxA )()()(1xxAAx|11AAAAAAxxx bxAAxAA )()(xAxAA )(xAxAAIA )(1xA

32、AAAIx 111)(Wait a minute Who said that(I+A 1 A)is invertible?(只要只要 A充分小，使得充分小，使得)1|11 AAAA|1|1|1111AAAAAAAAAAAAxx 是关键是关键的误差放大因子，称为的误差放大因子，称为A的的条件数条件数，记为，记为cond(A),越越 则则 A 越病态，越病态，难得准确解。难得准确解。|1 AA大大例：例：Hilbert 阵阵 12111131211211nnnnnHcond(H2)=27cond(H3)748cond(H6)=2.9 106cond(Hn)as n 注：注：现在用现在用Matlab

33、Matlab数学软件可以很方便求矩阵的状态数数学软件可以很方便求矩阵的状态数!定义定义2:2:设线性方程组的系数矩阵是非奇异的设线性方程组的系数矩阵是非奇异的,如果如果condcond(A A)越大越大,就称这个方程组就称这个方程组越病态越病态.反之反之,condcond(A A)越小越小,就称这个方就称这个方程组程组越良态越良态.一般判断矩阵是否病态，并不计算一般判断矩阵是否病态，并不计算A A 1 1，而由经验得出。，而由经验得出。行列式很大或很小（如某些行、列近似相关）；行列式很大或很小（如某些行、列近似相关）；元素间相差大数量级，且无规则；元素间相差大数量级，且无规则；主元消去过程中出

34、现小主元；主元消去过程中出现小主元；特征值相差大数量级。特征值相差大数量级。近似解的误差估计及改善：近似解的误差估计及改善：设设 的近似解为的近似解为，则一般有，则一般有bxA*x0*xAbr|*|xxrxbcond(A)误差上限误差上限 改善方法改善方法(1)：Step 1:近似解近似解 bxA;1xStep 2:;11xAbr Step 3:;111drdA Step 4:;112dxx 若若 可被精确解出，则有可被精确解出，则有 就是精确解了。就是精确解了。1dbAxAbAxx11112)(2x经验表明经验表明：若：若 A 不是非常病态（例如：不是非常病态（例如：），），则如此迭代可达到

35、机器精度；但若则如此迭代可达到机器精度；但若 A 病态，则此算法也病态，则此算法也不能改进。不能改进。1)(Acond 改善方法改善方法(2)(2)对方程组进行预处理，即适当选择非奇异对角阵对方程组进行预处理，即适当选择非奇异对角阵D D，C,C,使求解使求解 Ax=b Ax=b 的问题转化为求解等价方程组的问题转化为求解等价方程组 DACCDACC-1-1x=Db,x=Db,且使且使DACDAC 的条件数得到的条件数得到改善。改善。用双精度进行计算，以便改善和减轻病态矩阵的影响用双精度进行计算，以便改善和减轻病态矩阵的影响。)()()()()(1111CcondAcondDcondCCAADDDACDACDACcond