数值分析第四版第五章直接法

1、第五章. 解线性代数方程组的直接方法 1.引言,（一般有限步内得不到精确解）,解线性方程组的两类方法：,直接法: 经过有限次运算后可求得方程组精确解的方法,(不计舍入误差!),迭代法：从解的某个近似值出发，通过构造一个无穷,序列去逼近精确解的方法。,2. Gauss (高斯)消去法,相当于第i个方程-第一个方程数新的第i方程同解！第一方程不动！,Gauss 消去法计算过程,在A(1):b(1)中,红方框中的元素是要化为0的部分；A(2):b(2)中,红方框中的元素全部已发生变化,故上标由(1)改(2).,即对增广阵进行初等行变换:,(i=2,3,n),(i,j=2,3,n),(i=2,3,n)

2、,(i=2,3,n),第k次消元(1kn-1),设第k-1次消元已完成,且 0,此时增广矩阵如下:,本次消元的目的是对框内部分作类似第一次消元的处理,消掉第k+1个方程到第n个方程中的xk项,即把 到 化为零.计算公式如下:,(i=k+1,n),(i,j=k+1,n),(i=k+1,n),(i=k+1,n),只要 0,(k=1,2,n-1)消元过程就可以进行下去.当k=n-1时,消元过程完成,得:,它的方阵部分A(n)是一个上三角形矩阵,它对应的方程组是一个上三角形方程组,只要 0,就可以回代求解公式为:,(i=n-1,n-2,1),综合以上讨论,高斯消元法解线性方程的公式为: 1消元 令 (

3、i,j=1,2,n),(i=k+1,k+2,n),(i,j=k+1,k+2,n),(i=k+1,k+2,n),2回代,若 0,(i=n-1,n-2,1),对k=1到n-1,若 0 ,进行: (i=k+1,k+2,n),系数矩阵与常数项：,返回,(即为P169 定理6),消去第一列的 n-1 个系数要计算n*(n-1） 个乘法。,Gauss消去法乘法计算量,3. 高斯主元素消去法,3.1 列主元消元法 设已用列主元消元法完成Ax=b的第k-1次消元(1kn-1), 此时方程组Ax=bA(k)x=b(k),有如下形式:,若 =0,则必有方框内的元素全为零,此时易证,即方程组Ax=b无确定解,应给出

4、信息退出计算,(k j n),然后用高斯消元法进行消元. 这样从k=1做到k=n-1,就完成了消元过程,只要A0,列主元高斯消元法必可进行下去.,进行第k次消元前,先进行、两个步骤:,A=A(k)=0,若,0,则交换ik行和k行元素,即:,在方框内的一列内选出绝对值最大者,即确定ik,全主元消元法 在列主元消去法中,若每次选主元不局限所在列中,而在整个主子矩阵,中选取,便称为全主元高斯消元法,此时增加的步骤为: ,确定ik , jk, 若 =0, 给出=0的信息,退出计算.,作如下行交换和列交换 行交换:,(k j n),列交换:,(k i n),值得注意的是,在全主元的消元法中,由于进行了列

5、交换, x各分量的顺序已被打乱.因此必须在每次列交换的同时,让 机器“记住”作了一次怎样的交换,在回代解出后将x各分 量换回原来相应的位置,这样增加了程序设计的复杂性. 此外作步比较大小时,全主元消元法将耗用更多的机时. 但全主元消元法比列主元消元法数值稳定性更好一些.实际 应用中,这两种选主元技术都在使用.,选主元素的高斯消元法是一种实用的算法,它可以应用于任意的方程组Ax=b,只要A0.,高斯-若当(Gauss-Jordan)消元法 很多问题都需要计算方阵的逆阵.线性代数中有多种求逆方法,本节讨论求A-1的数值方法.A非奇异,求A-1,设X=A-1,AX=I,将X和I列分块,得n个线性方程

6、组: Axi=ei , i=1,2,n. 解这n个线性方程组就求得A-1,原则上,所有解的方程组的方法都能求A-1.实际计算中使用较多的是高斯-若当消元法.,高斯-若当消元法 高斯-若当消元法是高斯消元法的一种变形.高斯消元法是消去对角元下方的元素.若同时消去对角元上方和下方的元素,而且将对角元化为1,就是高斯-若当消元法.,设高斯-若当消元法已完成k-1步,于是Ax=b化为等价方程组A(k)x=b(k),增广阵,在第k步计算时,考虑将第k行上下的第k列元素都化为零,且akk化为1.对k=1,2,n 1按列选主元,确定ik使,2换行,交换增广阵第k行和第ik行,(j=k,k+1,n),3.计算

7、乘数 mik= - aikakk mkk=1akk.,4消元 aij=aij+mikakj(i=1,2,n,且ik,j=k+1,n) bi=bi+mikbk (i=1,2,n,且ik),5主行计算 akj=akjmkk (j=k,k+1,n) bk=bkmkk,当k=n时,显然xi=bi,i=1,2,n就是Ax=b的解.,(i=1,2,n, ik),高斯-若当消元法的消元过程比高斯消元法略复杂,但省去了回代过程,它的计算量约为n3/2,大于高斯消元法.也称为无回代的高斯消元法.,求方阵的逆 高斯-若当消元法解方程组并不比高斯消元法优越,但用于矩阵求逆是适宜的,实际上它是初等变换方法求逆的一种规

8、范化算法:,例3,所以,将高斯-若当消元法算法略加改造便成了求逆的算法. AX=I 对 k=1, 2, , n,1按列选主元, 确定ik,2换行, 交换第 k 行和 ik 行,(j=k,k+1,n),(j=1,2,n),3计算乘数 mik=-aikakk , (i=1,2, ik) mkk=1akk,5主行计算 akj = akj mkk ( j=k,k+1,n) xkj = xkj mkk ( j=1,2,n ) A-1=(xij)nn,应该注意到,在上述算法中都加进了选列主元的步骤, 它可以保证在方阵A非奇异的情况下都能求得A-1.求逆算 法中不宜使用全面选主元,否则会改变逆阵元素的排列,

9、 增加程序设计的复杂性.,4消元 aij = aij + mikakj,(i=1,2,n,且ik, j=k+1,n) xij = xij + mikxkj,(i=1,2,n, 且ik, j=1,2,n),Ax=b是线性方程组,A是nn方阵,并设A的各阶顺序主子式不为零。 令A(1)=A,当高斯消元法进行第一步后,相当于用一个初等矩阵左乘A(1).不难看出,这个初等矩阵为,其中mi1 (i=2,3,n) 由式(2.9)确定,即,5 矩阵的三角分解,同样第 k 步消元有,进行n-1步后,得到A(n) ,记U=A(n),显然U 的下三角部分全化 为零元素,它是一个上三角阵。整个消元过程可表达如下:,

10、其次指出,已知U是上三角矩阵,下面讨论L的性态,首先指出,L相当是由各L-1k (k=1,2,n-1)的所有左下元素拼凑 后加上对角元1而得.,当A进行LU分解后,Ax=b就容易解了. 即Ax=b等价于:,这样,Ax=b分解为解两个三角形方程组,三角形方程组是极易求解的.(*)即是高斯消元法的矩阵表述.,定理7 矩阵 Ann,只要A的各阶顺序主子式非零,则A可以分解为一个单位下三角阵L和一个上三角阵U的乘积,即A=LU,且这种分解是唯一的. (P173),L是下三角矩阵,且所有的对角元为1,故称为单位下三角阵.这样,A=LU称为A的LU分解,其中L是单位下三角阵,U是上三角阵.,例1,选主元素

11、的矩阵表示,直接LU分解法 A的LU分解可以用高斯消元法完成,但也可以用矩阵乘法原理推出另一种方法,结果是完全一致的,由矩阵乘法公式可推出直接LU分解算法,也称杜利特(Doolittle) 紧凑算法,现总结如下:,1.矩阵分解:A=LU，记,例如：,2.解 L y = b,3.解 U x = y,杜利特算法可用文字描述如下: 将A的元素划分为形如“”的n框, 如A的第一行和第一列元素为第一框,紧靠第一框内部的第二行和第二列元素为第二框,依此类推,ann一个元素为第n框.为便于描述我们称原A矩阵中元素为a元素,计算后每框中的每行元素为u元素(包含对角元素),每列元素为l元素(不包含对角元素).,

12、杜利特算法实际上就是高斯消元法的另一种形式.它的计算量与高斯消元法一样.但它不是逐次对A进行变换,而是一次性地算出L和U的元素.L和U的元素算出后,不必另辟存贮单元存放,可直接存放在A的对应元素的位置,节省存贮单元,因此也称为紧凑格式法.,也可以用增广矩阵A:b进行以上操作,这时每框要多算一个 u元素,结果矩阵中的最后一列就是 Ly=b 的解,回代时只要解 Ux=y 就行了.,当用手工计算小型线性方程组的解时,用此方法比较方便.,这样我们可以在A上直接修改计算,作如下的操作:,得到的矩阵由L和U的元素拼合而成,它的上三角部分(包含 对角元素)即是LU分解后的U矩阵;它的下三角部分(不包含 对角

13、元素)加上一个单位阵(即在对角元处全部写为1),就是L阵.,(1)第1框: u元素等于对应的a元素;l元素等于对应的a元素 除以本框对角元素;,(2)第2到n框: u元素等于对应的a元素减去已经算过的各框 中同行同列元素的乘积;l元素除进行以上操作外,还要除 以本框对角元素,方阵行列式求法 在实际问题中,有时会遇到求方阵的行列式.在线性代数中讲到的行列式定义算法和展开算法均不适用于阶数较高的行列式.而LU分解是计算行列式的十分方便和适用的算法. A=LU,即只要将U阵的对角元相乘就可得A的行列式. 为避免计算中断,还应加上选主元素过程,此时每做一次行交换(或列交换),行列式要改变一次符号,有:

14、,其中 p 是进行的行交换次数(对列主元消元法而言),或是行交换和列交换次数的总和(对全主元消元法而言).若选不出非零的主元素,则必有det A=0.,克劳特(Crout)分解 将LU分解换一个提法:要求L为一般下三角阵,U为单位上 三角阵,就是克劳特分解.,LT显然是单位上三角阵,记,这就是克劳特分解。,有,实际上将A转置为AT,进行LU分解 AT=LU 则 A=UTLT,显然当A的各阶顺序主子式非零时,它是存在且唯一的. 克劳特分解对应的解法称克劳特算法,也可用于解线性 方程组。它的特点是在回代时不做除法.它在下文的追 赶法中有应用.,UT显然是一般下三角阵,记,平方根法 实际问题中Ax=

15、b,A若是对称正定矩阵,则高斯消元法简化为 平方根法或改进的平方根法.,取,记,记,记,于是,矩阵的LDU分解,证: 由条件有A=LU,由分解过程知,定理 若A的各阶顺序主子式非零,则A可以分解为A=LDU, 其中L是单位下三角阵,U是单位上三角阵,D是对角阵,且这种 分解是唯一的.,定理 设 A 对称正定,则存在三角分解 A=LLT,L是非 奇异下三角矩阵,且当限定 L 的对角元为正时, 这种分解是唯一的.,证 因为 A 对称正定,所以A各阶顺序主子式为正.,所以有A=LDU,AT=UTDLT=A=LDU.,由LDU分解的唯一性知:UT=L, U=LT,所以 A=LDLT.,下面证D的对角元

16、为正数,因为 |L|=10,所以 LTyi=ei 必有非零解 yi 0,yiTAyi=yiTLDLTyi=(LTyi)TD(LTyi)=eiTDei=di,对称正定矩阵的乔累斯基(Cholesky)分解(P189),形如 A=LLT 的分解称为对称正定矩阵的乔累斯基分解。,这是一个二次型,由A 对称正定知 di 0 (i=1,2, ,n).,由证明过程容易看出,是唯一的,对称正定矩阵的A=LLT分解对应于解对称正定方程组Ax=b 的平方根法.,平方根,由矩阵乘法原理容易推出L的元素 lij 的算法:,对j=1,2,n,计算,Ax=b可化为,解法是:,这就是平方根法,它适用于系数阵是对称正定矩阵

17、的方程组.它的运算量以乘除法计是n3/6左右,只是高斯消元法的一半,显然这是因为只计算L不算U的缘故.,平方根法不用考虑选主元,这也是它的优点.它的缺点是要计算n次开平方,为避免开平方运算,发展了平方根法的改进形式,它对应于A=LDLT分解.(方法略),追赶法 (p193),在很多情况下,如三次样条插值,常微分方程的边值问题等 都归结为求解系数矩阵为对角占优的三对角方程组Ax=f, 即:,其中|i-j|1时,aij=0,且满足如下的对角占优条件: (1)|b1|c1|0,|bn|an|0 (2)|bi|ai|+|ci|, aici0, i=2,3,n-1.,用矩阵乘法比较和Crout分解可推导

18、总结算法步骤如下:,1.分解计算,2.解Ly=f,3.解Ux=y,实际计算中Ax=f的阶数往往很高,应注意A的存贮技术.已知数 据只用4个一维数组就可存完.即ai,bi,ci,fi各占一个 一维数组, i和i可存放在bi,ci的位置,yi 和xi 则可放在fi的位置,整个运算可在4个一维数组中运行.追赶 法的计算量很小,只是5(n-1)次乘除法.追赶法的计算也不要选 主元素.,5 向量和矩阵的范数,在分析方程组的解的误差及下章中迭代法的收敛时,常产 生一个问题,即如何判断向量x的“大小”,对矩阵也有类似的问 题.本节介绍n维向量和nn矩阵的范数.,向量范数,定义 x 和y 是 Rn 中的任意向

19、量, 向量范数是定义在Rn上的实值函数,它满足:,(1) x 0,并且,当且仅当x=0时, x =0;,(2) k x =|k| x , k是一个实数;,(3) x + y x + y ,容易看出, 实数的绝对值,复数的模, 三维向量的模都 满足以上三条, n维向量的范数概念是它们的自然推广.,常使用的向量范数有三种,设x=(x1,x2,xn)T,容易验证,它们都满足定义中的三个条件.,例 x=(1,0.5,0,-0.3)T, 求,解:,定义 设A是nn矩阵,定义,为矩阵A的范数.,矩阵范数,这样定义的范数有如下性质:,从向量范数出发,可以定义矩阵的范数.,(1) A0,并且,当且仅当A是零矩

20、阵时, A =0,(2) kA= |k|A,(3)两个同阶方阵A,B,有A+BA+B,(4)A是nn矩阵,x是n维向量,有AxAx,(5)A,B都是nn矩阵,有ABAB,矩阵范数最常用的有以下三种:,（1 是ATA的最大特征值）,它们分别与向量的三种范数对应,即用一种向量范数可定义.,相应的矩阵范数,定理15,即对任意给定的两种范数 有下列关系: c1 c2 其中的c1,c2是正的常数, 表示向量,也可以是 矩阵的范数.,Rn空间上的范数等价,从以上定理看出,当向量或矩阵的任一种范数趋于零时,其它,各种范数也趋于零.因此讨论向量和矩阵序列的收敛性时,可不指明使用的何种范数；证明时,也只要就某一

21、种范数,证明就行了.,序列的收敛,定义,如,称向量序列x (k)收敛于向量x. 记作:,有了向量和矩阵范数的概念,就可以定义向量和矩阵,定义,如,称常方阵序列A (k)收敛于方阵A 记作:,为任一种范数.,则称为A的谱半径.,定理19 方阵谱半径和方阵范数有如下关系:,定理16,其中,谱半径,定义 设n 阶方阵A的特征值为j (j=1,2,n),证:设i是A的任一特征值,xi为对应的特征向量,两边取范数,用矩阵性质(4)有:,|i|x iAx i,因为 x i 0,所以x i0,i=1,2,n,Axi= ixi,所以 |i |A,所以 (A)=max|i |A.,的充要条件是(A)1.,例,求

22、A1,A2 ,A, (A).,解: 显然A1 =4, A=4,57,定理 设A是nn阶矩阵,A的各次幂组成的矩阵序列 I,A,A2, ,Ak, 收敛于零,即,证明从略.,定理20 对于实对称矩阵A, 则,定理21 如果 ,则 为非奇异矩阵, 且,条件数及病态方程组 线性方程组 Ax=b 的解是由系数阵A及右端向量b决定的.由实际问题中得到的方程组中,A的元素和b的分量,总不 可避免地带有误差,因此也必然对解向量x产生影响.,这就提出一个问题:当A有误差A,b有误差b时,解向量x有多大误差?即当A和b有微小变化时,x的变化有多大? 若A和b的微小变化,也只导致x的微小变化,则称此问题是 “良态”

23、的；反之,若A和b的微小变化会导致 x 的很大变化,则称此问题为“病态”问题.,6 误差分析,在以下的讨论中, 设A非奇异, b0,所以|x|0.,1.设Ax=b中仅b向量有误差b ,对应的解x发生误差x , 即:,注意到Ax=b,所以Ax=b,若A非奇异,有 x=A-1b.,所以 xA-1 b .,又因 b=AxAx,所以,两式相除,有,即x 的相对误差小于等于b 的相对误差的AA-1倍.,2.A 有误差A ,b无误差,此时有类似的结论。,定义若nn方阵A非奇异,则称AA-1为A的 条件数,记为 Cond(A)=AA-1 由于选用的范数不同,条件数也不同,在有必要时, 可记为 Condp(A)=ApA-1p (p=1,2, ).(P207),3.由以