730 likes | 1.06k Views
第五章 解线性方程组的直接方法. §5.1 引 言. 线性方程组 :. (5.1). 1. 结束. 常记为矩阵形式 Ax = b (5.2). 在自然科学和工程技术中 , 很多问题归结为 解线性方程组 . 有的问题的数学模型中虽不直接表现为含线性方程组 , 但它的数值解法中将问题“离散化”或“线性化”为线性方程组 . 因此线性方程组的求解是数值分析课程中最基本的内容之一. 此时 A 是一个 n × n 方阵 , x 和 b 是 n 维列向量 .
E N D
第五章 解线性方程组的直接方法 §5.1 引 言 线性方程组: (5.1) 1 结束
常记为矩阵形式 Ax=b (5.2) 在自然科学和工程技术中,很多问题归结为解线性方程组.有的问题的数学模型中虽不直接表现为含线性方程组,但它的数值解法中将问题“离散化”或“线性化”为线性方程组.因此线性方程组的求解是数值分析课程中最基本的内容之一. 此时A是一个n×n方阵,x 和 b是n维列向量. 根据线性代数知识若|A | ≠0,(5.2)的解存在且唯一. 2 结束
关于线性方程组的解法一般分为两大类,一类是直接法,即经过有限次的算术运算,可以求得(5.1)的精确解(假定计算过程没有舍入误差).如线性代数课程中提到的克莱姆算法就是一种直接法.但该法对高阶方程组计算量太大,不是一种实用的算法〔1〕.实用的直接法中具有代表性的算法是高斯消元法,其它算法都是它的变形和应用.关于线性方程组的解法一般分为两大类,一类是直接法,即经过有限次的算术运算,可以求得(5.1)的精确解(假定计算过程没有舍入误差).如线性代数课程中提到的克莱姆算法就是一种直接法.但该法对高阶方程组计算量太大,不是一种实用的算法〔1〕.实用的直接法中具有代表性的算法是高斯消元法,其它算法都是它的变形和应用. 另一类是迭代法,它将(5.1)变形为某种迭代公式,给出初 始解 x0 ,用迭代公式得到近似解的序列{xk},k=0,1,2, ,在一定的条件下 xk→x* (精确解). 3 结束
迭代法显然有一个收敛条件和收敛速度问题. 这两种解法都有广泛的应用,我们将分别讨论,本章介绍 直接法. §5.2 高斯(Gauss)消元法 高斯消元法是一种古老的方法.我们在中学学过消元法,高斯消元法就是它的标准化的、适合在计算机上自动计算的一种方法。 4 结束
5.2.1 高斯消元法的基本思想 例1解方程组 (5.3) (5.4) (5.5) 第一步,将(5.3)乘-2加到(5.4);(5.3)乘-1加到(5.5), 得到 (5.3) (5.6) (5.7) 5 结束
第二步,将(5.6)乘-2/3加到(5.7),得到 (5.3) (5.6) (5.8) 回代:解(5.8)得x3,将x3代入(5.6)得x2,将x2, x3代入(5.3)得x1,得到解 x*=(2,1,-1)T 容易看出第一步和第二步相当于增广矩阵[A:b]在作行变换,用ri表示增广阵[A:b]的第i行: 6 结束
7 结束
由此看出上述过程是逐次消去未知数的系数,将Ax=b化为等价的三角形方程组,然后回代解之,这就是高斯消元法.由此看出上述过程是逐次消去未知数的系数,将Ax=b化为等价的三角形方程组,然后回代解之,这就是高斯消元法. 5.2.2 高斯消元法公式 记 Ax=b 为 A(1)x=b(1) , A(1) 和 b(1) 的元素记为 和 , i , j=1,2,,n .第一次消元,目的是消掉第二个方程到第n个方程中的 x1 项,得到 A(2)x=b(2),这个过程须假定 ≠0. 8 结束
在[A(1):b(1)]中,红方框中的元素是要化为0的部分;[A(2):b(2)]中,红方框中的元素全部已发生变化,故上标由(1)改(2),计算公式为:在[A(1):b(1)]中,红方框中的元素是要化为0的部分;[A(2):b(2)]中,红方框中的元素全部已发生变化,故上标由(1)改(2),计算公式为: 9 结束
(i=2,3,,n) (i,j=2,3,,n) (i=2,3,,n) (i=2,3,,n) 第 k 次消元 (1≤k≤n-1) 设第k-1次消元已完成,且 ≠0,此时增广矩阵如下: 10 结束
本次消元的目的是对框内部分作类似第一次消元的处理,消掉第k+1个方程到第n个方程中的xk项,即把 到 化为零.计算公式如下: 11 结束
(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时,消元过程完成,得: 12 结束
它的方阵部分A(n)是一个上三角形矩阵,它对应的方程组是一个上三角形方程组,只要 ≠0,就可以回代求解,公式为 (i=n-1,n-2,,1) 综合以上讨论,高斯消元法解线性方程的公式为: 13 结束
1消元 ①令 (i,j=1,2,…,n) ②对k=1到n-1,若 ≠0 ,进行: (i=k+1,k+2,,n) (i=k+1,k+2,,n) (5.9) (i,j=k+1,k+2,,n) (i=k+1,k+2,,n) 14 结束
2回代,若 ≠0 (i=n-1,n-2,,1) (5.10) 5.2.3 高斯消元法的条件 以上过程中,消元过程要求 ≠0 (i=1,2,,n-1),回代过程则进一步要求 ≠0,但就方程组Ax=b讲, 是否等于0是无法事先看出的. 15 结束
注意A的顺序主子式Di(i=1,2,,n)在消元过程中不变.这是因为消元所作的变换是“将某行的若干倍加到另一行”注意A的顺序主子式Di(i=1,2,,n)在消元过程中不变.这是因为消元所作的变换是“将某行的若干倍加到另一行” 上,据线性代数知识,此类变换不改变行列式的值.若高斯消元过程已进行了k-1步(此时当然应有 ≠0,i≤k-1),这时计算 A(k) 的顺序主子式: 16 结束
有递推公式 (i=2,3,,k) 显然, , 可知,消元过程能进行到底的充要 条件是Di≠0 ,(i=1,2,,n-1),若要回代过程也能完成,还应加 上Dn=|A|≠0,综合上述有: 定理5.1高斯消元法消元过程能进行到底的充要条件是系数阵A的1到n-1阶顺序主子式不为零;Ax=b能用高斯消元法解的充要条件是A的各阶顺序主子式不为零. 17 结束
5.2.4 高斯消元法的计算量估计 消元过程的工作量,参看公式(5.9),k是消元次数,k=1,2,,n-1,第k步消元时,计算 lik(i=k+1,,n)需要 n-k次除法;计算 (i,j=k+1, ,n)需要 (n-k)2 次乘法及 (n-k)2次加减法;计算 需要 n-k次乘法及 n-k次减法,合计: 乘除法次数: 18 结束
加减法的次数: 回代过程的工作量,参见公式(5.10),求 xk需 n-k 次加减法, n-k次乘法和1次除法,合计为 乘除法次数: 加减法次数: 19 结束
总的运算次数为 乘除法 (当n较大时) 加减法 (当n较大时) 一般讲乘除法的运算比加减法占用机时多得多,往往只统计乘除法次数而称高斯消元法的运算量为 次. 20 结束
§5.3 高斯主元素消元法 在上节的算法中,消元时可能出现 =0的情况,高斯消元法将无法继续;即使 ≠0,但 <<1,此时用它作除数,也会导致其它元素数量级严重增加,带来舍入误差的扩散,使解严重失真. 例2线性方程组 准确解是(1,1)T. 21 结束
现设我们的计算机为四位浮点数,方程组输入计算机后成为现设我们的计算机为四位浮点数,方程组输入计算机后成为 (5.11) 用高斯消元法: l12=0.2×106 , r2=r2-l12×r1 回代:x2=0.100 0×10=1, x1=0, 解严重失真. 22 结束
若将 r1和 r2交换. 消元,l12=0.5×10-5,r2=r2-l12×r1 回代 , x1=0.1000×10 , x2=0.1000×10 ,得到准确解. 23 结束
从上例中可以看出,对方程组作简单的行交换有时会显著从上例中可以看出,对方程组作简单的行交换有时会显著 改善解的精度.在实际使用高斯消元法时,常结合使用 “选 主元”技术以避免零主元或“小主元”出现,从而保证高斯消 元法能进行或保证解的数值稳定性. 5.3.1 列主元消元法 设已用列主元消元法完成Ax=b的第k-1次消元(1≤k≤n-1), 此时方程组Ax=b→A(k)x=b(k),有如下形式: 24 结束
(5.12) 进行第k次消元前,先进行①、②两个步骤: ①在方框内的一列内选出绝对值最大者,即 ,确定ik.若 =0,则必有方框内的元素 25 结束
全为零,此时易证|A|=|A(k)|=0,即方程组Ax=b无确定解,应给出信息退出计算.全为零,此时易证|A|=|A(k)|=0,即方程组Ax=b无确定解,应给出信息退出计算. ②若 ≠0,则交换ik行和k行元素即 (k j n) 然后用高斯消元法进行消元. 这样从k=1做到k=n-1,就完成了消元过程,只要|A|≠0,列主元高斯消元法必可进行下去. 26 结束
5.3.2 全主元消元法 在(5.12)中,若每次选主元不局限在方框列中,而在整个主子矩阵 中选取,便称为全主元高斯消元法, 此时增加的步骤为: ① ,确定ik, jk,若 =0,给出 |A|=0的信息,退出计算. 27 结束
②作如下行交换和列交换 行交换: (k j n) 列交换: (k i n) 值得注意的是,在全主元的消元法中,由于进行了列交换 ,x 各分量的顺序已被打乱.因此必须在每次列交换的同时, 让机器“记住”作了一次怎样的交换,在回代解出后将x各 分量换回原来相应的位置,这样增加了程序设计的复杂性. 此外作①步比较大小时,全主元消元法将耗用更多的机时. 28 结束
但全主元消元法比列主元消元法数值稳定性更好一些.实际应用中,这两种选主元技术都在使用.但全主元消元法比列主元消元法数值稳定性更好一些.实际应用中,这两种选主元技术都在使用. 选主元素的高斯消元法是一种实用的算法,它可以应用于任意的方程组Ax=b,只要|A|≠0. §5.4 高斯-若当(Gauss-Jordan)消元法 很多问题都需要计算方阵的逆阵.线性代数中有多种求逆方法,本节讨论求A-1的数值方法.A非奇异,求A-1 29 结束
设X=A-1,AX=I,将X和I列分块,得n个线性方程组: Axi=ei , i=1,2,,n. 解这n个线性方程组就求得A-1,原则上,所有解的方程组的方法都能求A-1.实际计算中使用较多的是高斯-若当消元法. 5.4.1 高斯-若当消元法 高斯-若当消元法是高斯消元法的一种变形.高斯消元法是消去对角元下方的元素.若同时消去对角元上方和下方的元素,而且将对角元化为1,就是高斯-若当消元法. 设高斯-若当消元法已完成k-1步,于是Ax=b化为等价方程组A(k)x=b(k),增广阵 30 结束
(5.13) 在第k步计算时,考虑将第k行上下的第k列元素都化为零,且akk化为1.对k=1,2,,n 1按列选主元,确定ik使 31 结束
2换行,交换增广阵第k行和第ik行 (j=k,k+1,,n), 3计算乘数 mik= - aik/akk(i=1,2,,n, i≠k) mkk=1/akk. (5.14) 4消元 aij=aij+mikakj(i=1,2,,n,且i≠k,j=k+1,,n) bi=bi+mikbk(i=1,2,,n,且i≠k) (5.15) 32 结束
5主行计算 akj=akj×mkk(j=k,k+1,,n) bk=bk×mkk (5.16) 当k=n时, 显然xi=bi,i=1,2,,n 就是 Ax=b的解. 33 结束
高斯-若当消元法的消元过程比高斯消元法略复杂,但省去了回代过程,它的计算量约为n3/2,大于高斯消元法.也称为无回代的高斯消元法.高斯-若当消元法的消元过程比高斯消元法略复杂,但省去了回代过程,它的计算量约为n3/2,大于高斯消元法.也称为无回代的高斯消元法. 5.4.2 求方阵的逆 高斯-若当消元法解方程组并不比高斯消元法优越,但用于矩阵求逆是适宜的,实际上它是初等变换方法求逆的一种规范化算法: 例3 34 结束
35 结束
所以 将高斯-若当消元法算法略加改造便成了求逆的算法. AX=I 对k=1,2,,n 1按列选主元,确定ik 2换行,交换第k行和ik行 (j=k,k+1,,n) (j=1,2,,n) 3计算乘数 mik=-aik/akk, (i=1,2,, i≠k) mkk=1/akk, (5.16) 36 结束
4消元 aij=aij+mik akj,(i=1,2,…,n,且i≠k,j=k+1,…,n) xij=xij+mik xkj,(i=1,2,…,n,且i≠k,j=1,2,…,n) (5.18) 5主行计算 akj=akj×mkk(j=k,k+1,…,n) xkj=xkj×mkk(j=1,2,…,n) (5.19) A-1=(xij)n×n 应该注意到,在上述算法中都加进了选列主元的步骤,它可以保证在方阵A非奇异的情况下都能求得A-1.求逆算法中不宜使用全面选主元,否则会改变逆阵元素的排列,增加 程序设计的复杂性. 37 结束
§5.5 矩阵的 LU分解 Ax=b是线性方程组,A是n×n方阵,并设A的各阶顺序主 子式不为零。令 A(1)=A,当高斯消元法进行第一步后,相当于 用一个初等矩阵左乘A(1).不难看出,这个初等矩阵为 其中li1 (i=2,3,…,n) 由式(5.9)确定,即 38 结束
同样第k步消元有 进行n-1步后,得到A(n) ,记U=A(n),显然U的下三角部分全化为零元素,它是一个上三角阵。整个消元过程可表达如下: 已知U是上三角矩阵,下面讨论L的性态,首先指出 39 结束
其次指出 证明见习题。 40 结束
L相当是由各L-1k (k=1,2,…,n-1)的所有左下元素拼凑后加上对角元1而得. L是下三角矩阵,且所有的对角元为1,故称为单位下三角阵.这样,A=LU称为A的LU分解,其中L是单位下三角阵,U是上三角阵. 定理5.2矩阵 An×n,只要A的各阶顺序主子式非零,则A可以分解为一个单位下三角阵L和一个上三角阵U的乘积,即A=LU,且这种分解是唯一的. 当A进行LU分解后,Ax=b就容易解了. 即Ax=b等价于: 这样,Ax=b分解为解两个三角形方程组,三角形方程组是极易求解的.(5.24)即是高斯消元法的矩阵表述. 41 结束
5.5.2 直接LU分解 A的LU分解可以用高斯消元法完成,但也可以用矩阵乘法原理推出另一种方法,结果是完全一致的.设: 由矩阵乘法公式可推出直接LU分解算法,也称杜利特(Doolittle) 算法,现总结如下: 1.矩阵分解:A=LU 对k=1,2,…,n 42 结束
2.解 L y = b 3.解 U x = y 杜利特算法实际上就是高斯消元法的另一种形式.它的计算量与高斯消元法一样.但它不是逐次对A进行变换,而是一次性地算出L和U的元素.L和U的元素算出后,不必另辟存贮单元存放,可直接存放在A的对应元素的位置,节省存贮单元,因此也称为紧凑格式法. 43 结束
杜利特算法可用文字描述如下: 将A的元素划分为形如“┏”的n框, 如A的第一行和第一列元素为第一框,紧靠第一框内部的第二行和第二列元素为第二框,依此类推,ann一个元素为第n框.为便于描述我们称原A矩阵中元素为a元素,计算后每框中的每行元素为u元素(包含对角元素),每列元素为l元素(不包含对角元素). 这样我们可以在A上直接修改计算,作如下的操作: (1)第1框:u元素等于对应的a元素;l元素等于对应的a元素除以本框对角元素; (2)第2到n框:u元素等于对应的a元素减去已经算过的各框中同行同列元素的乘积;l元素除进行以上操作外,还要除以本框对角元素 44 结束
得到的矩阵由L和U的元素拼合而成,它的上三角部分(包含得到的矩阵由L和U的元素拼合而成,它的上三角部分(包含 对角元素)即是LU分解后的U矩阵;它的下三角部分(不包含对角 元素)加上一个单位阵(即在对角元处全部写为1),就是L矩阵. 也可以用增广矩阵[A:b]进行以上操作,这时每框要多算一个 u 元 素,结果矩阵中的最后一列就是 Ly=b的解,回代时只要解 Ux=y就 行了. 当用手工计算小型线性方程组的解时,用此方法比较方便. 45 结束
5.5.3 方阵行列式求法 在实际问题中,有时会遇到求方阵的行列式.在线性代数中讲到的行列式定义算法和展开算法均不适用于阶数较高的行列式.而LU分解是计算行列式的十分方便和适用的算法. A=LU 即只要将U阵的对角元相乘就可得A的行列式. 为避免计算中断,还应加上选主元素过程,此时每做一次行交换(或列交换),行列式要改变一次符号,有: 46 结束
其中 p是进行的行交换次数(对列主元消元法而言),或是行交换和列交换次数的总和(对全主元消元法而言).若选不出非零的主元素,则必有 det A=0. 5.5.4 克劳特(Crout)分解 将LU分解换一个提法:要求L为一般下三角阵,U为单位上三角阵,就是克劳特分解. 实际上将A转置为AT,进行LU分解 AT=LU则 A=UTLT 47 结束
UT显然是一般下三角阵,记 有 LT显然是单位上三角阵,记 这就是克劳特分解。 显然当A的各阶顺序主子式非零时,它是存在且唯一的.克劳特分解对应的解法称克劳特算法,也可用于解线性方程组。它的特点是在回代时不做除法.它在下文的追赶法中有应用.限于篇幅,公式从略。 48 结束
§5.6 平方根法 实际问题中Ax=b,A若是对称正定矩阵,则高斯消元法简化为 平方根法或改进的平方根法. 5.6.1 矩阵的LDU分解 定理5.3若A的各阶顺序主子式非零,则A可以分解为A=LDU,其中 L 是单位下三角阵, U是单位上三角阵, D是对角阵, 且这 种分解是唯一的. 证: 由条件有A=LU,由分解过程知 取 49 结束
记 记 记 于是 5.6.2 对称正定矩阵的乔累斯基(Cholesky)分解 定理5.4设 A 对称正定,则存在三角分解 A=LLT,L是非奇异下三角矩阵 ,且当限定 L 的对角元为正时,这种分解是唯一的. 50 结束