第五章线性方程组的直接解法

1、1第五章第五章 线性方程组的直接解法线性方程组的直接解法 /* direct methord for solving linear systems */上一页上一页 下一页下一页 返回返回 第一节第一节 gaussgauss消去法消去法第二节第二节 直接三角分解方法直接三角分解方法第三节第三节 方程组的性态与误差估计方程组的性态与误差估计2求解求解bxa bax 为为线线性性方方程程组组的的矩矩阵阵形形式式)1(22112222212111212111 nnnnnnnnnnbxaxaxabxaxaxabxaxaxa nnnnnnbbbxxxaaaaa111111,., 0 解解唯唯一一若若 a

2、上一页上一页 下一页下一页 返回返回 31 gauss1 gauss消去法消去法一、一、 高斯顺序消去法高斯顺序消去法思思路路首先将首先将a化为上三角阵化为上三角阵 /* upper-triangular matrix */，再回代求解再回代求解 /* backward substitution */。= 是一种古老的求解线性方程组的方法是一种古老的求解线性方程组的方法, 按自然顺序进按自然顺序进行消元的方法行消元的方法. 上一页上一页 下一页下一页 返回返回 4例例1 解方程组解方程组 02115472321321321xxxxxxxxx解：解：step1 消元消元 012111154711

3、21312212rrrr 5 . 35 . 05 . 203330711223)35.2(rr 620033307112 6233372332321得得同同解解方方程程组组xxxxxx上一页上一页 下一页下一页 返回返回 5step2 对上三角形方程进行回代求解对上三角形方程进行回代求解, 得得 12/ )7(23/ )33(326321323xxxxxx 下面我们来一般性地讨论求解下面我们来一般性地讨论求解n阶线性方程组的高阶线性方程组的高斯顺序消去法斯顺序消去法.上一页上一页 下一页下一页 返回返回 6消元消元 ,)()0()0(nnijaaa 令令 )0()0(1)0(nbbbbstep

4、 1：设设 ，计算因子，计算因子0011)(a)., 2(/)0(11)0(11niaalii 将增广矩阵将增广矩阵/* augmented matrix */ 第第 i 行行 li1 第第1 1行行，得到，得到与与(1)(1)式等价的方程组式等价的方程组)0()0()1(bxa 式式变变为为则则 )0()0(1)0()0()0(2)0(1)0(1)0(12)0(11)0(,nnnnnnbbbbaaaaaaaa), 3 , 2,(,)0(11)0()1()0(11)0()1()1()1(njiblbbalaabxaiiijiijij 其其中中简简记记为为 )1()1(2)0(121)1()1(

5、2)1(2)1(22)0(1)0(12)0(1100nnnnnnnbbbxxxaaaaaaa上一页上一页 下一页下一页 返回返回 7.),(,)()(01011化化为为的的将将以以后后各各行行中中行行为为基基础础以以第第设设nkiakakikkkk)., 1(/)1()1(nkiaalkkkkikik 计计算算step 2：一般第一般第 k 次消元次消元 (1k n-1n-1) 具具有有形形式式其其中中程程组组式式等等价价的的方方即即已已计计算算出出与与步步计计算算已已完完成成设设第第)1()1()1(,)1(,1 kkkabxak )1()1()1()1()1(2)1(22)0(1)0(12

6、)0(11)1(knnknkkknkkknnkaaaaaaaaaa上一页上一页 下一页下一页 返回返回 8即即行行乘乘以以行行减减去去第第然然后后用用第第,iklki)()(,)(,)(,)()(,)()()()()(knnkknknkkkkkknkkkkkknkaaaaaaaaaaa1111111101012011上一页上一页 下一页下一页 返回返回 ),( )()()(nkkiblbbkkikkiki2111).,()()()(nkjalaakkjikkijkij111，9 )1()1(2)0(121)1()1(2)1(22)0(1)0(12)0(11.nnnnnnnnbbbxxxaaaa

7、aa)2()1()1( nnbxastep 3：继续上述过程继续上述过程, 且设且设 aii（i-1）0(i=1,2,n-1),直直到完成第到完成第 n-1 次消元次消元, 最后得到与最后得到与 a(0)x=b(0) 等价的三角等价的三角形方程组形方程组 a(n-1)x=b(n-1).将将(1)式化为式化为(2)式的过程称为式的过程称为消元过程消元过程.上一页上一页 下一页下一页 返回返回 10回代回代 )1 ,2,1()(),2()1(1)1()1()1()1( nnkaxabxabxkkknkjjkkjkkknnnnnn得得求求解解公公式式求求解解三三角角形形方方程程组组定理定理 若若a的

8、所有的所有顺序主子式顺序主子式 /* determinant of leading principal submatrices */ 均不为均不为0，则高斯消元无需换行即可进行到底，则高斯消元无需换行即可进行到底，得到唯一解。得到唯一解。iiiiiaaaaa.)det(1111 事实上，只要事实上，只要 a 非奇异，即非奇异，即 a 1 存在，则可通过逐次消存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出唯一解。元及行交换，将方程组化为三角形方程组，求出唯一解。上一页上一页 下一页下一页 返回返回 11高斯顺序消去法流程图高斯顺序消去法流程图1, 1, )1,(, dknnan

9、增增广广矩矩阵阵输输入入方方程程阶阶数数0 kkafkkikikaaanki/, 1 t停停输输出出失失败败信信息息,kjikijijaaaankj , 1 k=k+1kkadd 0 nnaf11111, ,/)(,nnkaaaaaaddkknkjjnkjnknknn停停输输出出失失败败信信息息,t结结束束系系数数行行列列式式即即输输出出,), 1)(1,dnkxaknk 消元过程消元过程回代过程回代过程上一页上一页 下一页下一页 返回返回 12.,)(步步的的主主元元为为第第系系数数其其首首项项个个方方程程为为主主方方程程步步消消元元时时保保留留的的第第称称第第定定义义kakkkkk1 顺序

10、消去法的缺点为当主元顺序消去法的缺点为当主元akk(k -1)=0时时, 消元过程消元过程不能继续进行不能继续进行. 或者当或者当akk (k -1) 0时时, 虽然消元过程可以虽然消元过程可以进行进行, 但若但若akk (k -1) 0时时, 时时, 会出现会出现很小的数作除数的现象很小的数作除数的现象,使舍入误差增大使舍入误差增大,导致解的严重导致解的严重失真失真.)1()1(/ kkkkikikaal计计算算上一页上一页 下一页下一页 返回返回 13例：例：解方程组解方程组 3200001. 100001. 02121xxxx/* 精确解为精确解为 */1, 121 xx用用gauss消

11、去法计算：消去法计算：1521024103000010100001020000101000010100001010000rr.1112xx,664102000010200000101000010100001010000.,0112xx回回代代求求解解15241050101000010100001010000103000010100001020000rr.101000010100000103000010100001020000.二、二、主元素消去法主元素消去法上一页上一页 下一页下一页 返回返回 两行互换得两行互换得，若将若将 2114 由上例可以看出由上例可以看出,为提高算法的数值稳定性为提高

12、算法的数值稳定性, 应选取绝对值应选取绝对值尽可能大的元素作主元尽可能大的元素作主元. 此此时时增增广广矩矩阵阵为为方方程程组组变变为为步步消消元元成成设设已已用用列列主主元元消消去去法法完完,)11(1)1()1( kkbxankk )1()1()1(2)0(1)1()1()1()1()1(2)1(22)0(1)0(12)0(11)1()1(),(knkkknnknkkknkkknnkkbbbbaaaaaaaaaba0|max|)1(,)1( knikkiikkkaa即即作作为为主主元元大大的的元元素素在在方方框框内内选选取取绝绝对对值值最最步步消消元元时时在在进进行行第第,)1(, kki

13、kak.,然然后后继继续续消消元元行行元元素素行行与与交交换换kik上一页上一页 下一页下一页 返回返回 按列部分选主元的消去法也称按列部分选主元的消去法也称列主元消去法列主元消去法。15 21 . 0301045132321321321xxxxxxxxx程程组组用用列列主主元元素素消消去去法法解解方方例例解：解： 211 .03010451321,ba 211 .0313210104521rr 255 .20112 .100104513125/35/1rrrr 112 .10255 .200104523rr 96.14 .100255 .2001045235 .2/2 .1rr 96.14

14、.1255 .201045332321xxxxxx 4 .122 .1321xxx得得解解上一页上一页 下一页下一页 返回返回 16 一些特殊情况一些特殊情况, 主元就在对角线上主元就在对角线上,不需选主元不需选主元. 元素满足如下条件的矩阵元素满足如下条件的矩阵),2, 1(0|1niaanijjijii 即对角线上每一元素的绝对值均大于同行其他各元素绝对即对角线上每一元素的绝对值均大于同行其他各元素绝对值之和值之和,这样的矩阵称为这样的矩阵称为按行按行严格对角占优矩阵严格对角占优矩阵,简称简称严格对角严格对角占优矩阵占优矩阵.例：例： 511283124a性质：严格对角占优矩阵必定非奇异性

15、质：严格对角占优矩阵必定非奇异.,)1(用用顺顺序序消消去去法法求求解解即即可可因因此此不不必必选选主主元元必必为为步步主主元元则则消消元元过过程程中中第第对对称称且且严严格格对对角角占占优优若若 kkkaka上一页上一页 下一页下一页 返回返回 17三、高斯三、高斯- -约当约当消去法消去法 （求矩阵逆）求矩阵逆）与与 gauss消去法消去法 的主要区别：的主要区别： 每一步不计算每一步不计算 lik ，而是先将当前主元，而是先将当前主元 akk(k-1) 变为变为 1; 把把 akk(k-1) 所在列的上、下元素全消为所在列的上、下元素全消为0;baxibxa1 这种方法是不是比这种方法是

16、不是比gauss消去法更好消去法更好? gauss消去法过程中，消去对角线下方和上方的元素，消去法过程中，消去对角线下方和上方的元素，称这种方法为称这种方法为高斯高斯- -约当消去法约当消去法.上一页上一页 下一页下一页 返回返回 这种方法不用回代过程了，这种方法不用回代过程了，看上去是要好些？看上去是要好些？18 运算量运算量 /* amount of computation */由于计算机中乘除由于计算机中乘除 /* multiplications / divisions */ 运算的时运算的时间远远超过加减间远远超过加减 /* additions / subtractions */ 运算

17、的时间，故运算的时间，故估计某种算法的运算量时，往往只估计估计某种算法的运算量时，往往只估计乘除的次数乘除的次数，而且通，而且通常以乘除次数常以乘除次数的的最高次幂最高次幂为运算量的为运算量的数量级数量级。 gauss消去法消去法:step k：设设 ，计算因子，计算因子且计算且计算0)1( kkka)., 1(/)1()1(nkiaalkkkkikik )., 1,()1()1()()1()1()(nkjiblbbalaakkikkikikkjikkijkij 共进行共进行n 1步步 ) 1() 1(2)0(121) 1() 1(2) 1(22)0(1)0(12)0(11.nnnnnnnnb

18、bbxxxaaaaaa)1()1(/ nnnnnnabx) 1., 1()1(1)1()1( niaxabxiiinijjiijiii(n k) 次次(n k)2 次次(n k) 次次(n k) (n k + 2) 次次nnnknknnk6523)2)(2311 消元时乘除次数：消元时乘除次数：1 次次(n i +1) 次次22)1(1211nninni 回代时乘除次数：回代时乘除次数：gauss消去法的总乘除次数为消去法的总乘除次数为 ，运算量为，运算量为 级。级。nnn31323 33n上一页上一页 下一页下一页 返回返回 19 gauss-jordan 消去法消去法: 运算量约为运算量约

19、为 。故通常只用于求逆矩阵，而不用于。故通常只用于求逆矩阵，而不用于解方程组。求逆矩阵即解方程组。求逆矩阵即 。23n 1| aiia上一页上一页 下一页下一页 返回返回 202 2 直接三角分解法直接三角分解法一、一、 lu分解法分解法 /* lu factorization. */就就 n=3的情况分析顺序消去法的消元过程的情况分析顺序消去法的消元过程. 333323122322211131211baaabaaabaaaba32011111111,irlriaalaiiii行行消消元元第第令令设设 )1(3)1(33)1(32)1(2)1(23)1(221131211)1()1(00baa

20、baabaaaba1011)1()1()1(3121)1(babamllm 记记)()()(122132321220aala令令设设23233rlr 行行消消元元第第 )2(3)2(33) 1(2) 1(23) 1(221131211000babaabaaayu111)1()1()2(32)2(yubamlm 记记上一页上一页 下一页下一页 返回返回 21可见可见, 消元过程相当于下述矩阵乘法运算：消元过程相当于下述矩阵乘法运算：)1()2(yubamm 由分块乘法可得：由分块乘法可得：(*)1()2()1()2(ybmmuamm 1211112)()()()()()()(mmmml令令直接计

21、算可得直接计算可得111101011011323121323121lllllll由由(*)式可得式可得blylua yuxblybluxbaxuxylua令令 )2(33)1(23)1(22131211000aaaaaau上一页上一页 下一页下一页 返回返回 22step 1:)0(/111111 aaalii记记m(1) =111121nll ，则，则 )0()0()1(bam)0(1)0(1)0(11.baan) 1(a) 1 (bstep n 1: ) 1() 1(2)0(1) 1() 1(2) 1(22)0(1)0(12)0(11) 1()2() 1(.nnnnnnnnnbbbaaaa

22、aabammm 其中其中 m (k)= 111, 1knkkll 以上分析推广到以上分析推广到n阶线性方程组可得阶线性方程组可得上一页上一页 下一页下一页 返回返回 23 1)(km111,1knkkll 1)1(1)2(1)1()(.)()(nmmm111jil,记为记为l单位下三角阵单位下三角阵/* unitary lower-triangular matrix */记记 u =) 1() 1(2) 1(2211211. nnnnnaaaaaaa 的的 lu 分解分解/* lu factorization */也称也称 a 的的doolittle分解分解 1112121nnlllllua

23、上一页上一页 下一页下一页 返回返回 24 道立特分解法道立特分解法 /* doolittle factorization */： lu 分解的紧凑格式分解的紧凑格式 /* compact form */ nnnnnnnnnuuuuuullaaaa2221121112111111.11.), 1(,11njaujj )., 32(1),21(1), 32(1),21(11111nillniuunianiaaijij，列列元元素素的的第第及及，行行元元素素的的第第可可以以求求出出矩矩阵阵，列列元元素素及及第第，行行元元素素的的第第因因此此由由lua 令令根据矩阵乘法法则根据矩阵乘法法则,先比较等

24、式两边第先比较等式两边第1行和第行和第1列元素有：列元素有：), 3 , 2(,1111111niuaaaliii 上一页上一页 下一页下一页 返回返回 25设已定出设已定出u 的第的第1行到第行到第k-1行的元素行的元素 nkjkkknjnjuuuuuuuuuu, 1, 11, 12222111211 l 的第的第1列到第列到第k-1列的元素列的元素 1,211,2121111knnnkkkklllllll).,21(), 1,(nkkilklnkkjukuikkj，元元素素列列的的第第与与行行元元素素的的第第现现在在确确定定 因因此此有有：的的内内积积列列向向量量的的第第与与行行向向量量的

25、的第第是是矩矩阵阵,)00,()00, 1,()(, 1211,21tjjkjjkjjkkkkkjuuuuujulllklkja 上一页上一页 下一页下一页 返回返回 26., 1,111nkkjuululakjkqqjkqnqqjkqkj )1(., 1,11nkkjulaukqqjkqkjkj 因此因此., 2, 1,111nkkiulululakkikkqqkiqnqqkiqik 同理有同理有)2(., 2, 1,11nkkiuulalkkkqqkiqikik 得得(1),(2)两式就是两式就是 lu分解的一般计算公式分解的一般计算公式, 其结果与高斯消去法所得其结果与高斯消去法所得结果

26、完全一样结果完全一样. 避免了中间过程的计算避免了中间过程的计算,故称为故称为a的直接分解公式的直接分解公式.,除除一一下下要要用用只只是是最最后后还还一一样样时时规规律律与与而而计计算算的的乘乘积积之之和和同同一一列列的的与与框框中中与与它它位位于于同同一一行行的的仅仅用用前前外外除除用用到到时时计计算算kkkjiksjkskjkjuululkau1上一页上一页 下一页下一页 返回返回 27框框第第框框第第框框第第框框第第nkulllllluuulluuuuluuuuunnnknnikiiknkjkkkknjknjk2121212122222211111211按上图逐框求出矩阵按上图逐框求出

27、矩阵a的的lu分解分解,这种计算方法称为这种计算方法称为紧凑格式法紧凑格式法。上一页上一页 下一页下一页 返回返回 .,nkkjulaukqqjkqkjkj111.,nkkiuulalkkkqqkiqikik211128定理定理 若矩阵若矩阵a非奇异非奇异, 则则a能分解为能分解为lu 的充分必要条件是的充分必要条件是a的的所有顺序主子式所有顺序主子式 均不为均不为0.定理定理 若非奇异矩阵若非奇异矩阵a有有lu 分解分解,则此分解是唯一的则此分解是唯一的.上一页上一页 下一页下一页 返回返回 29例：利用系数矩阵的例：利用系数矩阵的lu分解分解, 求解方程组求解方程组 45141267311

28、61234312111114321xxxx解：解：lu分解的紧凑格式为分解的紧凑格式为111u上一页上一页 下一页下一页 返回返回 113u112u114u11121l41431l61641l111222u011123u211324u1141332l51611142l214233u2325061343l12141134)(u2152315261744)(u 21512201111112356114111lu30 4514126123561141114321yyyy则求解原方程组等价于求解下面两个方程组则求解原方程组等价于求解下面两个方程组ly=b及及ux=y 43214321215122011

29、111yyyyxxxx15466y上一页上一页 下一页下一页 返回返回 2121x31二、二、三对角方程组三对角方程组追赶法追赶法 nnnnnnndddxxxbacbacbacb212111122211dax 即即若若a满足满足gauss消去法可行的条件，则可用消去法可行的条件，则可用lu分解法求解分解法求解111132nllllnnuccucuu12211niclbuniualbuiiiiiii, , , , 32321111其中：其中：上一页上一页 下一页下一页 返回返回 32解方程组解方程组ax=d分为两步，即求解分为两步，即求解ly=d和和ux=y，计算公式如下：，计算公式如下：niy

30、ldydyiiii, , 32111上述方法为求解上述方法为求解三对角方程组的追赶法三对角方程组的追赶法，也称，也称thomas算法算法.上一页上一页 下一页下一页 返回返回 12111, , nnixcyuxuyxiiiiinnn 追赶法追赶法公式简单，计算量和存储量都小，整个求解过程公式简单，计算量和存储量都小，整个求解过程只需要只需要5n- -4次乘除运算。次乘除运算。33上一页上一页 下一页下一页 返回返回 三、三、平方根法平方根法 对称对称 正定正定矩阵的分解法矩阵的分解法将将对称对称 正定阵正定阵 a 做做 lu 分解分解记为记为) 3(dm 1112221111112221122

31、211211uuuuuuuuuuuuuuuunnnnnnnnldmlua 于于是是)(tttdlmaaa 的的对对称称性性有有由由有有分分解解的的唯唯一一性性可可知知由由分分解解的的因因此此上上式式也也是是为为上上三三角角阵阵为为单单位位下下三三角角阵阵其其中中,luluadlmtt.,ttlmlm 即即34定理定理 设设n阶对称正定矩阵阶对称正定矩阵a，则存在唯一的单位下三角阵，则存在唯一的单位下三角阵l及对及对角阵角阵d 使得使得 。tdlla 称为矩阵称为矩阵a 的的乔累斯基分解乔累斯基分解上一页上一页 下一页下一页 返回返回 定理定理 设矩阵设矩阵a对称正定，则存在唯一的对角元为正的下

32、三角阵对称正定，则存在唯一的对角元为正的下三角阵 l，使得，使得 。tlla 称为对称正定矩称为对称正定矩阵阵a 的的乔累斯基分解乔累斯基分解 利用利用乔累斯基（乔累斯基（cholesky）分解）分解式来求解式来求解ax=b的方法的方法也称也称cholesky方法或方法或平方根法平方根法353 方程组的性态与误差估计方程组的性态与误差估计上一页上一页 下一页下一页 返回返回 一、矩阵的条件数一、矩阵的条件数 例，考查以下三个方程组及其准确解例，考查以下三个方程组及其准确解 121 1.00012.0001112xx121 1.00012.0002112yy121 1.00052.0001112

33、zz1211xx 1202yy 121.80.2zz其准确解其准确解其准确解其准确解其准确解其准确解 可以看到，后两个方程组与第一个方程组相比，系数矩阵可以看到，后两个方程组与第一个方程组相比，系数矩阵或右端向量仅有或右端向量仅有0.0005以下的误差，但准确解却相差很大。对以下的误差，但准确解却相差很大。对这样的方程组，无论用多么稳定的算法求解，一旦计算中产生这样的方程组，无论用多么稳定的算法求解，一旦计算中产生误差就使解面目全非，所以该方程组的误差就使解面目全非，所以该方程组的性态很差性态很差。36上一页上一页 下一页下一页 返回返回 定义：定义：若方程组若方程组ax= =b的系数矩阵的系

34、数矩阵a与右端向量与右端向量b的微小变化的微小变化（小扰动）（小扰动）, ,将引起解向量将引起解向量x产生巨大变化，则称此方程组为产生巨大变化，则称此方程组为病病态方程组态方程组，其系数矩阵，其系数矩阵a称为称为病态矩阵病态矩阵，否则称，否则称ax=b为为良态方良态方程组程组，称，称a为为良态矩阵良态矩阵 . . 方程组的病态程度与方程组的病态程度与ax=b对对a和和b的扰动的的扰动的敏感程度敏感程度有关。有关。37求解求解 时，时，a 和和 的误差对解的误差对解 有何影响？有何影响？bxa bx 设设 a 精确，精确， 有误差有误差 ，得到的解为，得到的解为 ，即，即bb xx bbxxa

35、)(bax 1 |1bax 绝对误差放大因子绝对误差放大因子|xaaxb 又又|1bax |1bbaaxx 相对误差放大因子相对误差放大因子上一页上一页 下一页下一页 返回返回 38 设设 精确，精确，a有误差有误差 ，得到的解为，得到的解为 ，即，即ba xx bxxaa )( bxxaxxa )()( )(1xxaax |11aaaaaaxxx bxaaxaa )()(axxaa )(axxaaia )(1axaaaix 111)( (只要只要 a充分小，使得充分小，使得)1|11 aaaa |1|1|1111aaaaaaaaaaaaxx 是关键是关键的误差放大因子，称为的误差放大因子，称

36、为a的的条件数条件数，记为，记为cond (a) ,越越 则则 a 越病态，越病态，难得准确解。难得准确解。|1 aa大大上一页上一页 下一页下一页 返回返回 39注注: cond (a) 的具体大小与的具体大小与 | | 的取法有关，但相的取法有关，但相对大小一致。对大小一致。 cond (a) 取决于取决于a，与解题方法无关。，与解题方法无关。常用条件数有：常用条件数有：cond 1(a)cond (a)cond2 (a)(/ )(minmaxaaaatt 特别地，若特别地，若 a 对称，则对称，则|min|max)( acond2条件数的性质：条件数的性质： a可逆，则可逆，则 cond p (a) 1； a可逆，可逆， r 则则 cond ( a) = cond (a) ； a正交，则正交，则 cond 2 (a) =1； a可逆，可逆，r正交，则正交，则 cond 2 (ra) = cond 2 (ar) = cond2 (a) 。上一页上一页 下一页下一页 返回返回 40精确解精确解为为.11 x例例 97.199.1,98.099.099.01ba计算计算cond 2(a) 。 10000990099009800a 1 = 解：解：考察考察 a 的特征根的特征根 0)det(ai 000