河南农业大学《数值计算方法》教案 汪松玉 - 1 - 第五章 线性代数方程组的直接解法 教学目的 1. 掌握解线性方程组的高斯消去法、 高斯选主元素消去 法; 2. 掌握用直接三角分解法解线性方程组的方法; 3. 了解解对称正定矩阵线性方程组的平方根法与解三 对角线方程组的追赶法; 4. 掌握向量,矩阵范数,矩阵的条件数等概念及方程组 的扰动分析. 教学重点及难点 重点:1. 解线性方程组的高斯消去法、高斯选主元素 消去法; 河南农业大学《数值计算方法》教案 汪松玉 - 2 - 2. 直接三角分解法解线性方程组的方法; 3. 向量,矩阵范数,矩阵的条件数等概念及方 程组的扰动分析; 难点:方程组的扰动分析. 河南农业大学《数值计算方法》教案 汪松玉 - 3 - 线性代数是数值计算方法的基础,解线性代数方程组的 有效方法在计算数学和科学计算中具有特殊的地位和作用. 如弹性力学、电路分析、热传导和振动、以及社会科学及定 量分析商业经济中的各种问题等许多实际问题的求解最终 转化为代数问题来处理,即归结为解线性方程组,如插值公 式的建立,微分方程离散格式的构造等.n 阶线性方程组的一 般形式是 11 1 12 2 1 1 21 1 22 2 2 2 1 1 2 2 n n n n n n nn n n a x a x a x b a x a x a x b a x a x a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??? ? (5.0.1) 其矩阵形式表示为 ? Ax b 河南农业大学《数值计算方法》教案 汪松玉 - 4 - 其中 11 12 1 1 1 21 22 2 2 2 1 2 , , n n n n nn n n a a a x b a a a x b a a a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A x b 一般 b≠0.当系数矩阵 A 非奇异(即det A ≠0)时, 方程组(5.0.1) 有唯一解 x . 求解上述线性代数方程组的方法主要有两类: (a) 直接法:对于给定的方程组,在没有舍入误差的假 设下, 能在预定的运算次数内求得精确解.最基本的直接法是 Gauss 消去法, 重要的直接法全都受到 Gauss 消去法的启发. 缺点:计算代价高. 河南农业大学《数值计算方法》教案 汪松玉 - 5 - (b) 迭代法:基于一定的递推格式,产生逼近方程组精 确解的近似序列, 用某种极限过程去逐步逼近精确解的方法. 收敛性是其为迭代法的前提,此外,存在收敛速度与误差估 计问题. 本章介绍解线性方程组的直接法(主要是高斯消去法及 其各种变形), 即如果计算过程中没有舍入误差, 经过有限步 算术运算可求得方程组的精确解,但在实际计算中由于舍入 误差的存在和影响,这类方法也只能求得方程组的近似解. 其中选主元消去法,三角分解法是解低阶稠密矩阵方程组的 有效方法. 河南农业大学《数值计算方法》教案 汪松玉 - 6 - 5.1 高斯消去法 5.1.1 高斯消去法的基本方法 高斯消去法是一种古老的方法,其基本思想是通过对方 程组(5.0.1)作下列运算(变换): (1) 交换方程组中任意两个方程的次序; (2) 方程组中任何一个方程乘以某一个非零常数; (3) 方程组中任何一个方程减去另一个方程的常数倍. 在系数矩阵 A 非奇异的假设下,可以得到一个与(5.0.1)等价 (或者说是具有相同解)的三角形方程组 河南农业大学《数值计算方法》教案 汪松玉 - 7 - (1) (1) (1) (1) (1) 11 1 12 2 1, 1 2 1, 1 (2) (2) (2) (2) 22 2 2, 1 2 2, 2 ( 1) ( 1) ( 1) 1, 1 2 1, 1 ( ) , n n n n n n n n n n n n n n n n n n n a x a x a x a x b a x a x a x b a x a x b a x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ( ) n n b ? ? ? ? ? ? ? ? ? ? (5.1.1) 这个过程称为消元. 显然,只要全部 ( ) , 0 k k k a ? (k=n, n-1,?,1),就可依次求出 1 1 , , , x x x n n ? ? , 这种按未知量编号从大到小的次序求解上三角 形方程组的过程称为回代. 这种将 n 阶线性方程组(5.0.1)化简为三角形方程组 (5.1.1)之后再回代求解的方法称之为高斯(Gauss)消去法. 河南农业大学《数值计算方法》教案 汪松玉 - 8 - 下面以(5.0.1)为例说明高斯消去法的步骤: 第一步 若11 0 a ? ,令1111 , ( 2,3, , ) i i a l i n a ? ? ? ,用1il?乘第 一个方程分别加到第i ( 2,3, , ) i n ? ? 个方程上去, 从而消去第 i ( 2,3, , ) i n ? ? 个方程的首项,方程组(5.0.1)变为下列形式: (1) (1) (1) (1) 11 1 12 2 1 1 (2) (2) (2) 22 2 2 2 (2) (2) (2) 2 2 n n n n n nn n n a x a x a x b a x a x b a x a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ? (5.1.2) 其中, 1 (1) 1 , 1,2, , j j a a j n ? ? ? , (1) 1 1 b b ? (2) 1 1 , , 2,3, , ij ij i j a a l a i j n ? ? ? ? , 河南农业大学《数值计算方法》教案 汪松玉 - 9 - (2) (1) 1 1 , 2,3, , i i i b b l b i n ? ? ? ? 第二步 若(2) 22 0 a ? ,令(2) 2 2 (2) 22 , ( 3,4, , ) i i a l i n a ? ? ? ,用1il?乘第二个方程分别加到第i ( 3,4, , ) i n ? ? 个方程上去, 从而消去 (2) 2 , ( 3,4, , ) i a i n ? ? . 一般地,设第 1 k ? 步后方程组化为如下形式 河南农业大学《数值计算方法》教案 汪松玉 - 10 - (1) (1) (1) (1) 11 1 12 2 1 1 (2) (2) (2) 22 2 2 2 n n n n k k k kk k kn n k k k k nk k nn n n a x a x a x b a x a x b a x a x b a x a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ? ???? ? (5.1.3) 若()0kkk a ? ,令()(),(1, 2, , ) k ik ik k kk a l i k k n a ? ? ? ? ? ,用1il?乘第k 个 方程分别加到第i 个(1, 2 , ) i k k n ? ? ? ? 方程上去,得到如下 方程组 河南农业大学《数值计算方法》教案 汪松玉 - 11 - (1) (1) (1) (1) 11 1 12 2 1 1 (2) (2) (2) 22 2 2 2 ( 1) ( 1) ( 1) 1, 1 1 1, 1 ( 1) ( 1) ( 1) , 1 1 n n n n k k k k k k k n n k k k k n k k nn n n a x a x a x b a x a x b a x a x b a x a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ? ???? ? (5.1.4) 其中 ( 1) , ( , 1, 2, , ) k k k ij ij ik kj a a l a i j k k n ? ? ? ? ? ? ? ( 1) , ( 1, 2, , ) k k k i i ik k b b l b i k k n ? ? ? ? ? ? ? 经过 1 n ? 步消元后,最后从(5.1.4)可得到一个上三角形方程 组 河南农业大学《数值计算方法》教案 汪松玉 - 12 - (1) (1) (1) (1) 11 1 12 2 1 1 (2) (2) (2) 22 2 2 2 ( ) ( ) n n n n k k k kk k kn n k n n nn n n a x a x a x b a x a x b a x a x b a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ? ???? (5.1.5) 现考虑回代过程的一般计算公式. 当()0nnn a ? 时,逐步回代得原方程组的解 河南农业大学《数值计算方法》教案 汪松玉 - 13 - ( ) ( ) 1 / , ( )/ , ( 1, 2, ,1) n nn k kj kk n n n n k k k k j j k x b a x b a x a k n n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.1.6) 因为高斯消去法的消元过程,实际上是对方程组的系数 进行运算,与未知量毫无关系,在实质上是对方程组(5.0.1) 的增广矩阵[A, b]进行初等行变换, 直至将系数矩阵化为上三 角形. 为了将高斯消去法编制成实用的计算机程序,就必须对 消元过程的具体计算步骤作出一些规定, 使之有统一的规律. 由不同的规定原则便产生了不同形式的算法.例如顺序高斯 消去法、列主元高斯消去法以及全主元高斯消去法,这些都 河南农业大学《数值计算方法》教案 汪松玉 - 14 - 是高斯消去法中最典型且常用的算法.其共同特点为: 它们的 主要步骤都是由消元和回代两大过程组成,主要区别在于消 元过程的计算步骤有各自的规范,而这些规范是与程序的简 明通用以及计算的稳定性有关的. 5.1.2. 顺序高斯消去法 顺序高斯消去法是按方程和未知量的自然顺序进行消元 的,它是各种高斯消去法的基础,下面举例说明具体计算步 骤.为方便起见我们约定用记号②- * ① k 表示第②行减去第 ①行的 k 倍,用①? ②表示交换第①与第②行,其余类推. 例5.1 用顺序高斯消去法解方程组 河南农业大学《数值计算方法》教案 汪松玉 - 15 - 1 2 3 1 2 3 1 2 3 2 2 3 3, 4 7 7 1, 2 4 5 7, x x x x x x x x x ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? (5.1.7) 解 先作消元过程.增广矩阵为 [A,b]= 2 2 3 3 4 7 7 1 2 4 5 7 ? ? ? ? ? ? ? ? ? ? ? ? 1 2 ? ? ???? ? 第 列消元 ② ① ③+① 2 2 3 3 0 3 1 5 0 6 8 4 ? ? ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 16 - 2 2 ? ? ???? ? 第 列消元 ③ ② 2 2 3 3 0 3 1 5 0 0 6 6 ? ? ? ? ? ? ? ? ? ? ? 由此得上三角形方程组 1 2 3 2 3 3 2 2 3 3, 3 5, 6 6, x x x x x x ? ? ? ? ? ? ? ? ? ? ? ? (5.1.8) 再作回代过程.由(5.1.8)第三个方程求得 3 1 x ? ,代入第二个方 程得 2 3 ( 5 )/3 2 x x ? ? ? ? ? .将23,xx代入第一个方 程可得 1 2 3 (3 2 2 )/ 2 2 x x x ? ? ? ? ,所以方程组的解为 1 2 3 2, 2, 1 x x x ? ? ? ? . 上例的求解用的是最基本的顺序高斯消去法,它的消元过程 河南农业大学《数值计算方法》教案 汪松玉 - 17 - 遵守如下规则: (1) 依从左到右、自上而下的次序将主对角元下方的元 素化为零; (2) 不作行交换,也不作将某行乘以一个非零数的变换; (3) 当进行第 k 列消元时(即使主对角元 kk a 下方化零), 将kk a 下面各行分别减去第 k 行的适当倍数,不作其他变换. 例5.1 主要是用手算表示出指定算法的计算步骤,对数 据形式未作具体要求,其实在计算机里上述整数或分数均为 浮点数,计算繁简程度相同. 现在,我们将例 5.1 中使用的基本高斯消去法推广到解 一般的 n n ? 阶线性方程组(5.0.1).设方程组的系数矩阵非奇 异,其增广矩阵记为 河南农业大学《数值计算方法》教案 汪松玉 - 18 - 11 12 1 21 22 2 1 2 (1) (1) (1) (1) 1 (1) (1) (1) (1) (1) (1) 2 (1) (1) (1) (1) [ , ] n n n n nn n a a a b a a a b b a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b (5.1.9) 其中, 1 (1) 1 , j j a a ? (1) , i i b b ? ( , 1,2, , ) i j n ? ? . 高斯消元法的消元过程由 1 n ? 步组成: 第一步 若(1) 11 0 a ? ,令(1) 1 1 (1) 11 , ( 2,3, , ) i i a l i n a ? ? ? ,用1il?乘(5.1.9)第一行分别加到第 i ( 2,3, , ) i n ? ? 行上去,从而消去 (1) (1) (1) 21 31 1 , , , n a a a ? ,得到 河南农业大学《数值计算方法》教案 汪松玉 - 19 - 1 11 12 1 2 22 2 (1) (1) (2) (1) (2) (2) (2) 2 (2) (2) (2) (2) (2) 0 [ , ] 0 n n nn n n b a a a b a a A b a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.1.10) 其中, (2) (1) (1) 1 1 , ( , 2,3, , ) ij ij i j a a l a i j n ? ? ? ? , (2) (1) (1) 1 1 ,( 2,3, , ) i i i b b l b i n ? ? ? ? 第二步 若(2) 22 0 a ? ,令(2) 2 2 (2) 22 , ( 3,4, , ) i i a l i n a ? ? ? ,用1il?乘(5.1.10)第二行并分别加到第i ( 3,4, , ) i n ? ? 行上去,从而消 去42 2 (2) (2) (2) 32 , , , n a a a ? .依次下去,假设进行了 1 k ? 步,得到 河南农业大学《数值计算方法》教案 汪松玉 - 20 - 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 , 1 1, 1, 1 1, 1 , 1 [ , ] k n k k k n k k k k k k kk k k kn k k k k k k k k k k n k k k k k nk n k nn n a a a a a b a a a a b a a a b a a a b a a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b ? ? ? ? ? (5.1.11) 第k步假设 ( ) 0 k kk a ? , ( ) ( ) , ( 1, 2 , ) k ik ik k kk a l i k k n a ? ? ? ? ? , 按照上述 方法消去(5.1.11)式第 k 列的 ( ) ( ) 1 , , k k k k nk a a ? ? ,得到 河南农业大学《数值计算方法》教案 汪松玉 - 21 - 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 ( 1) ( 1) , 1 ( 1) ( 1) ( 1) 1, 1 1, 1 ( 1) ( 1) ( 1) , 1 [ , ] k n k k k n k k k k k k kk k k kn k k k k k k k n k k k k n k nn n a a a a a b a a a a b a a a b a a b a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b ? ? ? ? ? ? (5.1.12) 其中 ( 1) ( 1) , , 1, 2, , , 1, 2, , k k k ij ij ik kj k k k i i ik k a a l a i j k k n b b l b i k k n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.1.13) (5.1.13)式是消元过程的一般计算公式, 式中作为分母的元素 ( ) k kk a 称为第 k 步的主元素(简称主元), ( 1, 2 , ) ik l i k k n ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 22 - 称为乘子.若()0kkk a ? ,则()()1,,kkkknk a a ? ? 中至少有一个元素,例如()krk a 不为零, (否则, 方程组(5.0.1)的系数矩阵 A 奇异).这样, 我们可取 ( ) k rk a 为主元,然后交换矩阵 ( ) ( ) , k k ? ? ? ? A b 的第k 行与第 r 行,把它交换到( , ) k k 的位置上. 经过 1 n ? 步消元后, 最后从(5.1.12)可得到一个上梯形矩 阵 河南农业大学《数值计算方法》教案 汪松玉 - 23 - 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 , 1 ( 1) ( 1) ( 1) 1, 1 1, 1 ( ) ( ) [ , ] k n k k k n n n k k k k kk k k kn k k k k k k k n k n n nn n a a a a a b a a a a b a a a b a a b a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b (5.1.14) 这里,我们假设整个消元过程中没有进行矩阵的行交换. ( ) n A 是一个上三角矩阵,与()()[,]nnAb相应的上三角方程组(5.1.5) 和(5.0.1)同解. 使用高斯消去法求解n 阶线性方程组总共需要乘除法运 算次数为 3 2 1 ( 3 - ) 3 n n n ? 次, 减法运算次数为 3 2 1 (2 3 -5 ) 6 n n n ? 次. 河南农业大学《数值计算方法》教案 汪松玉 - 24 - 由于一般计算机中乘除运算时间远远超过加减运算时间,所 以在估计一种计算方法的运算量时, 往往只估计乘除的次数. 又因为n 很大时, 1 p p n n ? ,所以在比较计算量的大小时, 常常只比较n 的最高次项.因此,我们可以说高斯消去法的乘 除次数为 3 1 3 n 数量级. 从上面讨论不难看出, 高斯消去法的消去过程把方程组 化成上三角形方程组, 必须通过回代才能得到解.如果在消元 的过程稍作改变,即将对角线上方的元素也化为零,则可使 方程(5.0.1)化简成为每个方程只含有一个未知量的方程组, 因此求解时无需回代的步骤.这种消去法称为高斯-约当 (Gauss-Jordan)消去法,也称为无回代的高斯消去法. 但是, 经过计算我们发现,这种方法所需要的计算量太大,其中乘 河南农业大学《数值计算方法》教案 汪松玉 - 25 - 除法次数为 3 2 1 1 2 2 n n n ? ? ,即为 3 1 2 n 数量级,超过高斯消去法, 所以一般来说, 求解线性方程组时采用高斯-约当消去法是不 明智的.但是,用它来求解矩阵的逆矩阵时却很方便,有兴趣 的读者请参阅其他相关文献. 5.1.3 选主元高斯消去法 上一小节我们介绍的是顺序高斯消去法,其消元过程能 进行到底的充要条件是主元素 ( ) 0 k kk a ? ,( 1,2, , ) k n ? ? (其中 (1) 11 11 a a ? ), 因为它们在计算公式中是被作为除数的.但即使在 消元过程中主元素 ( ) k kk a ( 1,2, , ) k n ? ? 全不为 0, 也不能保证能 够避免在计算过程中产生巨大的误差.因为在计算过程中, 若主元素()kkk a 的绝对值相对其它元素(如与同列的元素河南农业大学《数值计算方法》教案 汪松玉 - 26 - ( ) ( ) 1 , , k k k k nk a a ? ? 相比)过小, 在这种情况下, 将导致舍入误差的扩 散,最后使计算结果不可靠,甚至消元过程无法进行,下面 举例说明. 例5.2 考察方程组-5 1 2 1 2 10 1 2 x x x x ? ? ? ? ? ? ? , 其精确解为1100000/99999 x ? , 2 99998/99999. x ? 若使用四位浮点十进 制运算进行消元,得到 -5 1 2 -5 -5 2 10 + =1 10 =10 x x x ? ? ? ? ? 利用回代过程得方程组的解为: 1 2 0, 1 x x ? ? .若用消去法求 河南农业大学《数值计算方法》教案 汪松玉 - 27 - 解12-5 1 2 2 10 1 x x x x ? ? ? ? ? ? ? 则可得 1 1 x ? , 2 1 x ? . 显然,第一个结果是严重失真的,而第二个结果则比较 精确.究其原因,就是小主元的出现.用它来作除数带入了大 的舍入误差,再经传播,误差变得更大,造成计算结果与真 解相差甚远,因此第一种算法是不稳定的. 在顺序高斯消去法中,出现零主元将导致求解失败;而 出现小主元又会导致计算数值解精度的严重失真.为避免出 现不稳定的算法,我们在使用消去法求解方程组(5.0.1)时, 在消元过程中引入逐次选主元的思想,其基本思路是每次消 元之前在系数矩阵中按一定的范围选取绝对值最大的元素 作为主元素,以便减少舍入误差的影响,这种方法称为选主 元高斯消去法, 是对顺序高斯消去法的一种改进.常用的选主 河南农业大学《数值计算方法》教案 汪松玉 - 28 - 元高斯消去法是列主元高斯消去法、行主元消去法和全主元 高斯消去法. 设方程组(5.0.1)的系数矩阵A非奇异,进行第 k-1 次消 元后,增广矩阵为 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 , 1 1, 1, 1 1, 1 , 1 [ , ] k n k k k n n n k k k k kk k k kn k k k k k k k k k k n k k k k k nk n k nn n a a a a a b a a a a b a a a b a a a b a a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A ? ? ? ? ? ? ? ? ? ? ? b ? ? ? ? ? (5.1.15) 列主元高斯消去法:是在第 k 次(k=1,2,… ,n-1)消元时, 从系数矩阵第 k 列中 ( 1) 1, k k k a ? ? 以下的各元素中寻找绝对值最大 河南农业大学《数值计算方法》教案 汪松玉 - 29 - 的元素作为主元素(称为列主元素),增广矩阵中主元所在的 行称为主行, 主元所在的列称为主列.我们的工作就是要确定 r,使()krk a 为()kkk a 1 , , k k k k nk a a ? ? 中之最大者. 若r>k 则交换增广矩阵的第 k 行和第 r 行,使()krk a 代替 ( ) k kk a 成为主元素,然后进行消元.其余步骤与顺序高斯消去法 相同.由(5.1.15)可以看出,只要 det A≠0,就必存在 ( , r r k ? 1, , ) k n ? ? ,使得 ( ) k r k a ≠0,因此整个计算过程一定能进行到 底. 行主元高斯消去法:是在第 k 次(k=1,2,… ,n-1)消元时, 若从 1 2 , , , k k k kk kk kn a a a ? ? ? 中选取绝对值最大的元素作为主元素 (称为行主元素),即若 河南农业大学《数值计算方法》教案 汪松玉 - 30 - ( ) ( ) max k k k kj kj k j n a a ? ? ? 则选取 ( ) k k kj a 作主元,并且在进行第 k 次消元前交换增广矩阵 的第 k 列和 k j 列(必须记录这种交换信息,以便整理方程组 的解时用).经这样修改的高斯消去法称为行主元高斯消去法. 全主元高斯消去法也简称全主元消去法,它与列主元消 去法的不同之处,是在第 k 次(k=1,2, ?,n-1)消元时,从系 数矩阵(见(5.1.15))的右下角(n-k+1)阶子矩阵 ( ) ( ) ( ) ( ) k k kk kn k k nk nn a a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? 中,选取绝对值最大的元素作为主元素,如果它位于第 r 行 河南农业大学《数值计算方法》教案 汪松玉 - 31 - 第s列, 则通过交换 k, r 两行及交换 k, s 两列, 可以使主元素 位于 ( ) k kk a 的位置, 然后进行消元计算. 由于作列的交换改变 了方程中未知量的次序, 因此回代过程实际上是按未知量调 换后的编号顺序求解的. 下面举例说明计算步骤. 例5.3 用列主元高斯消去法求解下面的方程组(用三位 有效数字计算) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 9 5 . 0 4 5 5 2 5 . 0 5 5 . 0 2 01 . 0 3 2 1 3 2 1 3 2 1 x x x x x x x x x 解 用方程组的增广矩阵进行消元,每次选出的主元用 小方框标出 河南农业大学《数值计算方法》教案 汪松玉 - 32 - 0.01 2 0.5 5 [ , ] 1 0.5 2 5 9 5 4 0.5 A b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ③ ① ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 5 5 9 5 . 0 2 01 . 0 2 5 . 0 1 5 . 0 4 5 (第1列消元, 在0.01,-1,5 中取绝对值最大者 5 作主元, 交换①、③行使 5 换到 11 a 的位置, 再作消元计算.) 河南农业大学《数值计算方法》教案 汪松玉 - 33 - ? ? ? ? ? ? ? ? ? ? ? 500 1 5 1 ① ③ ① ② 5 4 0.5 9 0 1.30 2.10 6.80 5.02 0 2.01 0.501 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?③ ② ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 80 . 6 02 . 5 9 10 . 2 30 . 1 0 501 . 0 01 . 2 0 5 . 0 4 5 (第2列消元,在-1.30 和2.01 中取绝对值最大者 2.01 作 主元,交换②、③行使 2.01 换到 22 a 的位置,再作消元计算.) 河南农业大学《数值计算方法》教案 汪松玉 - 34 - ③ ② ? ? ? ? ?????? ? 130 201 . . 5 4 0.5 9 0 2.01 0.501 5.02 3.55 0 0 1.78 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 以上完成了消元过程, 得到一个上三角形方程组.再作回 代可求得 3 2 1 1.99, 2.00, 0.00100 x x x ? ? ? ? 事实上这个方程组的准确解是 2 , 2 , 0 3 2 1 ? ? ? ? x x x ,可见 用列主元消去法求得的近似解是令人满意的. 作为比较,用顺序高斯消去法重解例 5.3,对方程组的 增广矩阵作顺序消元得 河南农业大学《数值计算方法》教案 汪松玉 - 35 - 0.01 2 0.5 5 1 0.5 2 5 5 4 0.5 9 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 100 500 0.01 2 0.5 5 0 200 48.0 495 0 1000 251 2510 ? ? ? ? ? ? ? ? ? ? ???? ? ? ? ? ? ? ? ? ? ? ② ① ③ ① 5.00 0.01 2 0.5 5 0 200 48.0 495 0 0 11.0 30.0 ? ? ? ? ? ? ? ? ????? ? ? ? ? ? ? ? ? ? ③ ② 再作回代求得 3 2 1 2.73, 1.82 , 1.00 x x x ? ? ? ? ,这样的结果误 差太大,所以是无效的. 河南农业大学《数值计算方法》教案 汪松玉 - 36 - 例5.4 用全主元高斯消去法求解下面的方程组(用四位 有效数字计算) 1 2 3 1 2 3 1 2 3 12 3 3 15, 18 3 15, 6, x x x x x x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 它的准确解为 1 2 3 1, 2, 3 x x x ? ? ? . 解 用方程组的增广矩阵进行消元,每次选出的主元用 小方框标出 12 3 3 15 [ , ] 18 3 1 15 1 1 1 6 A b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 37 - ? ???? ? ① ② 18 3 1 15 12 3 3 15 1 1 1 6 ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 ?????? ? 第 列消元 18 3 1 15 0 1 2.333 5.000 0 1.167 0.944 5.167 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?????? ? 交换②、③列18 1 3 15 0 2.333 1 5.000 0 0.944 1.167 5.167 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 38 - ?????? ? 第2列消元 18 1 3 15 0 2.333 1 5.000 0 0 1.572 3.144 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 消元完成.由于交换过②、 ③两列, 所以改变了方程中 2 x 和3x的次序,最后得到与原方程组等价的上三角形方程组 1 3 2 3 2 2 18 3 15 2.333 5.000 1.572 3.144 x x x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? 回代即可逐步解出: 2 3 1 2.000, 3.000, 1.000. x x x ? ? ? 实际计算和理论分析表明, 列主元高斯消去法具有较好 的数值稳定性,是目前求解中小型稠密线性方程组最常用的 河南农业大学《数值计算方法》教案 汪松玉 - 39 - 方法之一. 在实际计算时,有些具有特殊系数矩阵的方程组,按顺 序消去法,若能保证主元的绝对值最大,或者设舍入误差能 得到控制,则不必选主元.这类矩阵有对称的严格对角占优 阵,对称正定阵等. 定义 5.1 所谓严格对角占优阵,是指 A 的元素满足如下 条件的矩阵: 1 1,2 , . n ii ij j j i a a i n ? ? ? ? ? ? 可以证明,严格对角占优矩阵必定是非奇异阵,而且经 过一步高斯消去法后仍然保持对角占优的优势. 河南农业大学《数值计算方法》教案 汪松玉 - 40 - 5.2 高斯消去法的矩阵形式 和矩阵的 LU 分解 顺序高斯消去法是对方程组(5.0.1)的增广矩阵 11 12 1 21 22 2 1 2 (1) (1) (1) (1) 1 (1) (1) (1) (1) (1) (1) 2 (1) (1) (1) (1) [ , ] n n n n nn n a a a b a a a b b a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b 进行一系列初等变换来实现的,根据线性代数知识,初等变 换与初等矩阵有密切关系的.所以,顺序消元法有如下形式: 第一步 若11 0 a ? ,令1111 ,( 2,3,4 , ) i i a l i n a ? ? ? ,组成初 河南农业大学《数值计算方法》教案 汪松玉 - 41 - 等下三角阵 21 1 1 1 0 1 0 1 n l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L (5.2.1) 则111 12 1 2 22 2 (1) (1) (2) (1) (2) (2) (2) 2 (1) (1) (2) (2) 1 (2) (2) (2) 0 0 n n nn n n b a a a b a a a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L A b A b 即用 1 i l ? 乘(5.1.9)的第一行分别加到第 ( 2,3, , ) i i n ? ? 行上去, 相当于用 1 L 左乘增广矩阵. 河南农业大学《数值计算方法》教案 汪松玉 - 42 - 第二步 若(2) 22 0 a ? ,令(2) 2 2 (2) 22 , ( 3,4 , ) i i a l i n a ? ? ? ,组成 2 32 2 1 0 0 1 0 1 0 1 n l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L (5.2.2) 用1il?乘(5.1.10)第二行并分别加到第 i( 3,4 , ) i n ? ? 行上去, 相当于用 2 L 左乘 (1) (1) 1[ , ] L A b . 整个消元过程是从左到右、自上而下地将主对角元下 方的 21 31 1 32 42 2 1 n n n n a a a a a a a ? ? ? ? 化为 0,这就相当于 河南农业大学《数值计算方法》教案 汪松玉 - 43 - 用形如(5.2.1),(5.2.2)的初等矩阵 1 2 3 1 , , , , n? ? L L L L 依次左乘 (1) (1) [ , ] A b .经过 1 n ? 步消元后,最后得到 (1) (1) 1 2 2 1[ , ] n n ? ? ? L L L L A b 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 , 1 ( 1) ( 1) ( 1) 1, 1 1, 1 ( ) ( ) k n k k k n k k k k kk k k kn k k k k k k k n k n n nn n a a a a a b a a a a b a a a b a a b a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ( ) ( ) [ , ] n n ? A b 河南农业大学《数值计算方法》教案 汪松玉 - 44 - 于是有 (1) ( ) 1 2 2 1 (1) ( ) 1 2 2 1 , . n n n n n n ? ? ? ? ? ? ? ? L L L L A A L L L L b b (5.2.3) 注意到 ( 1,2, , 1) i i n ? ? ? L 是可逆的且为主对角元全为 1 的下 三角阵,而有限个单位下三角阵的乘积仍是单位下三角阵, 故可令 1 1 2 1 ( ) n ? ? ? ? L L L L , 11 12 1 (1) (1) (1) (1) 1 (2) (2) (2) 2 22 2 ( ) ( ) , n n n n nn n b a a a b a a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? U y 则得: 1 (1) (1) ? ? L A b U y ,或1(1) 1 (1) , ? ? ? ? L A U L b y .高斯消 河南农业大学《数值计算方法》教案 汪松玉 - 45 - 去法的回代过程就是求解 ? Ux y . 又由 1 (1) ? ? L A U 可得: (1) ? A LU ,这说明,在()0kkk a ? 的 条件下,通过高斯消去法可以把方程组(5.0.1)的系数矩阵 A 分解为单位下三角阵和上三角阵的乘积, 称为 A 的LU 分解. 一旦实现了矩阵 A 的LU 分解,那么求解方程组 A ? x b 的问题就等价于求解两个三角形方程组: (1) L ? y b ,求y;(2) U ? x y ,求x.定义 5.2 设A为n阶矩阵,若A有分解式 ? A LU ,其中L为下三角矩阵,U 为上三角矩阵,则称之为 A 的LU 分 解或三角分解.特别若 L 是单位下三角矩阵, 则该分解式称为 杜利特尔(Doolittle)分解,如 河南农业大学《数值计算方法》教案 汪松玉 - 46 - 11 12 1 21 22 2 1 2 1 1 1 n n n n nn u u u l u u l l u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A 若U 是单位上三角矩阵,则称之为克劳特(Crout)分解,如A=????????????????????????1112112 2 1 22 21 11 ? ? ? ? ? ? ? n n nn n n u u u l l l l l l 的分解式就是克劳特分解式. 一个矩阵不一定存在三角分解式,当存在时一般不唯 河南农业大学《数值计算方法》教案 汪松玉 - 47 - 一.如2445?????????A,0112???????B,则B不存在三角分解(请自 行用反证法证明), A 的三角分解式不唯一,容易验证 2 0 1 2 4 1 0 3 ? ? ? ? ? ? ? ? ? ? ? ? ? ? A 是一个三角分解式,任取可逆的对角矩 阵2005???????D,可得新的三角分解式 1 2 0 1 2 4 1 0 3 ? ? ?? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ?? ? A D D 2 0 2 0 1/ 2 0 1 2 4 1 0 5 0 1/5 0 3 ? ?? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? 河南农业大学《数值计算方法》教案 汪松玉 - 48 - ? ? ? ? ? ? ? ? ? ? ? ? ? ? 5 / 3 0 1 2 / 1 5 8 0 4 那么,一个矩阵在什么条件下存在唯一的三角分解式呢? 定理 5.1 n 阶矩阵 A 存在唯一的杜利特尔分解和克劳特 分解的充要条件是 A 的顺序主子矩阵 k A (k=1,2,? ,n)非奇异. 推论 若定理 5.1 的条件满足,则A存在唯一的 LDU 分解式 ? A LDU ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 1 1 1 1 1 2 1 12 2 1 2 1 21 ? ? ? ? ? ? ? n n n n n u u u d d d l l l (5.2.4) 其中 L 为单位下三角阵,U 为单位上三角阵, D 为对角阵. 河南农业大学《数值计算方法》教案 汪松玉 - 49 - 反之也真. 设方程组(5.0.1)的系数矩阵 A 非奇异,并有三角分解 ? A LU .求解线性方程组的三角分解法的基本思想是:可以 先把 A 作LU 分解,存储 L 和U ,方程组就化为 x ? LU b 从而使问题转化为求解两个简单的三角形方程组 ? Ly b, ? Ux y (5.2.5) 再根据各个方程组不同的右端向量分别求解(5.2.5)的两个三 角形方程组即可. 下面首先介绍杜利特尔分解法.其步骤为: 第一步 三角分解 对k=1,2,…,n,计算U 和L的元素.为了在计算机上求解, 需找出一般计算公式.设A的杜利特尔分解式 LU A ? 为 河南农业大学《数值计算方法》教案 汪松玉 - 50 - ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? nn n n n n nn n n n n u u u u u u l l l a a a a a a a a a ? ? ? ? ? ? ? ? ? ? ? ? ? 2 22 1 12 11 2 1 21 2 1 2 22 21 1 12 11 1 1 1 (5.2.6) 由矩阵乘法知 A 的每个元素等于 L 的某一行乘U 的某一列, 若能按一定次序每次取出一个这样的等式,使其中只含一个 未求的元素,便可一步步直接算出 L 和U 的元素. 首先容易看出 A 和U 的第一行元素对应相等,即11,(1,2, , ) j j a u j n ? ? ? (5.2.7) 现在 11 u 已求, 故由 A 的第一列元素 11 1 1 u l a i i ? 可求出 L 的第一 列元素 1 1 11 / , ( 2,3, , ) i i l a u i n ? ? ? (5.2.8) 河南农业大学《数值计算方法》教案 汪松玉 - 51 - 是否可以如此继续,依次交替求出U 的第二行、 L 的第 二列, U 的第三行、L 的第三列…?不妨设U 的前 1 ? k 行和 L 的前 1 ? k 列已求, 欲求U 的第 k 行元素 ) ( k j ukj ? , 为此考虑 A 中对应元素 kj a , kj a 为L的第 k 行与U 的第 j 列对应位置元素 乘积之和,即kj a = 1 , 1 [ 1 0 0 ] k k k l l ? ? ? 1 1, j k j kj nj u u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 52 - = 1 1 , 1 1, k j k k k j kj l u l u u ? ? ? ? ? ? 因11,,?kkkll?及jkjuu11,,??分别为 L 中前 1 ? k 列、 U 中前 1 ? k 行元素,故上式唯一的未知数便可求出,即11,(,1, , ) k k j k j k s s j s u a l u j k k n ? ? ? ? ? ? ? ? (5.2.9) 再求 L 的第 k 列元素 ) ( k i lik ? ,考虑 A 中对应元素 ) ( k i aik ? ,与上面类似地可得 1 1 ( ) , ( 1, 2, , ) k ik is sk s ik k k a l u l i k k n u ? ? ? ? ? ? ? ? ? (5.2.10) 编程时对 n k , , 2 , 1 ? ? 循环使用(5.2.9)、(5.2.10)便可依次交替 计算U 的第 k 行、 L 的第 k ( n k , , 2 , 1 ? ? )列,直至全部算出. 河南农业大学《数值计算方法》教案 汪松玉 - 53 - 计算公式为 1 1 1 1 , , 1, , ( ) / , 1, 2, , k k j k j k s s j s k ik ik i s sk k k s u a l u j k k n l a l u u i k k n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.2.11) 第二步 求解过程 (1) 解?Ly b.对k=1,2,…,n 计算 1 1 k k k k j j j y b l y ? ? ? ?? , (5.2.12) (2) 解?Ux y .对k=n,n-1,…,1 计算 河南农业大学《数值计算方法》教案 汪松玉 - 54 - 1 ( )/ n k k k j j k k j k x y u x u ? ? ? ? ? . (5.2.13) 为了便于记忆公式并说明计算U 的行和 L 的列的次序,列 出下面的紧凑格式表: 11 11) ( u a 12 12 ) ( u a 13 13 ) ( u a 14 14 ) ( u a ? ① 1 1) ( y b 21 21) ( l a 22 22 ) ( u a 23 23 ) ( u a 24 24 ) ( u a ? ③ 2 2 ) ( y b 31 31) ( l a 32 32 ) ( l a 33 33 ) ( u a 34 34 ) ( u a ? ⑤ 3 3 ) ( y b 41 41) ( l a 42 42 ) ( l a 43 43 ) ( l a 44 44 ) ( u a ? ⑦ 4 4 ) ( y b ? ? ? ? ? ② ④ ⑥ ⑧ 其中①,②,?表示计算U 的行和 L 的列的次序,计算过程是 河南农业大学《数值计算方法》教案 汪松玉 - 55 - 按照一行一列、二行二列…的次序进行. 由公式(5.2.7)–(5.2.10)可知,上表中U 的第一行等于 A 的第一行,其余元素 j k u 等于 A 的对应元素 j k a 减去左边 L 的 同一行元素与上边U 的同一列元素对应乘积之和;计算 L 的 元素与计算U 的元素方法类似, 只是最后要除以同一列U 的 对角元.例如, 21 21 11 24 24 21 14 43 43 41 13 42 23 33 / , , ( )/ l a u u a l u l a l u l u u ? ? ? ? ? ? 等等.解方程组时,右端常数列b 作为矩阵的最后一列,计算 y 与计算U 的方法类似. 例 用三角分解法解线性方程组 河南农业大学《数值计算方法》教案 汪松玉 - 56 - . 7 17 3 5 3 0 1 0 3 4 2 1 1 0 1 0 0 2 0 1 4 3 2 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? x x x x 解 列表 1 11 ? u 0 12 ? u 2 13 ? u 0 14 ? u 5 1 ? y ① 0 21 ? l 1 22 ? u 0 23 ? u 1 24 ? u 3 2 ? y ③ 1 31 ? l 2 32 ? l 2 33 ? u 1 34 ? u 6 3 ? y ⑤ 0 41 ? l 1 42 ? l 0 43 ? l 2 44 ? u 4 4 ? y ⑦ ② ④ ⑥ 11 11 12 12 13 13 14 14 1 1 1, 0, 2, 0, 5; u a u a u a u a y b ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 57 - ; 0 1 0 , 1 1 1 , 0 1 0 11 41 41 11 31 31 11 21 21 ? ? ? ? ? ? ? ? ? a a l a a l a a l , 1 0 0 1 12 21 22 22 ? ? ? ? ? ? ? u l a u , 0 2 0 0 13 21 23 23 ? ? ? ? ? ? ? u l a u , 1 0 0 1 14 21 24 24 ? ? ? ? ? ? ? u l a u ; 3 5 0 3 1 21 2 2 ? ? ? ? ? ? ? b l b y 32 31 12 42 41 12 32 42 22 22 2 1 0 1 0 0 2, 1; 1 1 a l u a l u l l u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? , 2 0 2 2 1 4 23 32 13 31 33 33 ? ? ? ? ? ? ? ? ? ? ? u l u l a u , 1 1 2 0 1 3 24 32 14 31 34 34 ? ? ? ? ? ? ? ? ? ? ? u l u l a u ; 6 3 2 5 1 17 2 32 1 31 3 3 ? ? ? ? ? ? ? ? ? ? ? y l y l b y 河南农业大学《数值计算方法》教案 汪松玉 - 58 - 43 41 13 42 23 43 33 0 0 2 1 0 ; 2 a l u l u l u ? ? ? ? ? ? ? ? ? ? 44 44 41 14 42 24 43 34 3 0 0 1 1 0 1 2, u a l u l u l u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 4 4 41 1 42 2 43 3 7 0 5 1 3 0 6 4. y b l y l y l y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 由表可知 1 0 0 0 1 0 2 0 5 0 1 0 0 0 1 0 1 3 , , 1 2 1 0 0 0 2 1 6 0 1 0 1 0 0 0 2 4 L U ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? y 由U ? x y ,即 河南农业大学《数值计算方法》教案 汪松玉 - 59 - 1 2 3 4 1 0 2 0 5 0 1 0 1 3 0 0 2 1 6 0 0 0 2 4 x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 可得 4 2 4 x ? ,于是 4 3 2 1 2, 6 2 2 2, 3 1 2 1, 5 2 2 1 x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? . 例5.5 用杜利特尔分解法求解方程组 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 2 3 4 2 4 9 16 10 8 27 64 44 16 81 256 190 x x x x x x x x x x x x x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? . 河南农业大学《数值计算方法》教案 汪松玉 - 60 - 解设?ALU ,按紧凑格式表计算 L ,U 以及方程组 b = y L 的解 y ,得(1)1 (2)2 (3)3 (4)4 (2)2 (1)1 (4)2 (9)6 (16)12 (10)8 (1)1 (8)3 (27)6 (64)24 (44)18 (1)1 (16)7 (81)6 (256)24 (190)24 即1123421102612 8 , , 1 3 1 6 24 18 1 7 6 1 0 24 24 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L U y . 再求解 ? Ux y 求得原方程组的解为 河南农业大学《数值计算方法》教案 汪松玉 - 61 - 4 3 2 1 1 , 1, 1, 1. x x x x ? ? ? ? ? ? 接下来我们介绍矩阵的 Crout 分解,并导出计算公式. 设A可分解为下三角阵和单位上三角阵的乘积,即11 12 1 21 22 2 1 2 n n n n nn a a a a a a a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A 11 12 1 21 22 2 1 2 1 1 1 n n n n nn l u u l l u l l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? , 河南农业大学《数值计算方法》教案 汪松玉 - 62 - 则对 1,2, , i n ? ? 有1111,(,1, , ) 1, 2, , ) i ji ji jk ki k i i j ik ik kk ii k l a l u j i i n u a l u u j i i n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.2.14) 仿照杜利特尔分解的记忆方法,我们给出下面的紧凑格式 表: 11 11 ( ) a l 12 12 ) ( u a 13 13 ) ( u a 14 14 ) ( u a ? ② 1 1) ( y b 21 21) ( l a 22 22 ( ) a l 23 23 ) ( u a 24 24 ) ( u a ? ④ 2 2 ) ( y b 31 31) ( l a 32 32 ) ( l a 33 33 ( ) a l 34 34 ) ( u a ? ⑥ 3 3 ) ( y b 41 41) ( l a 42 42 ) ( l a 43 43 ) ( l a 44 44 ( ) a l ? ⑧ 4 4 ) ( y b 河南农业大学《数值计算方法》教案 汪松玉 - 63 - ? ? ? ? ? ① ③ ⑤ ⑦ 其中①,②, ?表示计算 L 的列和U 的行的次序, 计算过程是 按照一列一行、二列二行?的次序进行. 由公式(5.2.14)可知,上表中 L 的第一列等于 A 的第一 列,其余元素 j k u 等于 A 的对应元素 j k a 减去左边 L 的同一行 元素与上边U 的同一列元素对应乘积之和; 计算U 的元素与 计算 L 的元素方法类似, 只是最后要除以同一行 L 的对角元. 解方程组时,右端常数列b 作为矩阵的最后一列,计算 y 与 计算U 的方法类似. 例5.6 用Crout 分解法求解例 5.5 的线性方程组. 解 按紧凑格式,并对 b 作统一处理,有(1)1 (2)2 (3)3 (4)4 (2)2 河南农业大学《数值计算方法》教案 汪松玉 - 64 - (1)1 (4)2 (9)3 (16)6 (10)4 (1)1 (8)6 (27)6 (64)4 (44)3 (1)1 (16)14 (81)36 (256)24 (190)1 即1123421201364,,166143114 36 24 0 1 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L U y 再求解 ? Ux y 求得原方程组的解为 4 3 2 1 1 , 1, 1, 1. x x x x ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 65 - 5.3 解三对角线性方程组的追赶法 利用三角分解,易于导出一些特殊方程组的解法. 但在 一些实际问题中,例如常微分方程数值解、求热传导方程及 三次样条插值函数等问题中都会遇到系数矩阵是三对角矩 阵的方程组,并且系数矩阵大都有对角占优性质,因此不必 选主元,就能保证算法的稳定性. 设有方程组 ? Ax d ,其中 A 为三对角矩阵,即1112222333111,nnnnnnbcdabcdabcabcdab?????????????????????????????????????????Ad(5.3.1) 河南农业大学《数值计算方法》教案 汪松玉 - 66 - 关于三对角矩阵 A ,有如下性质: 定理 5.2 当三对角矩阵(5.3.1)满足条件 1 1 0 0 , 0,( 2,3, , 1) 0 k k k k k n n b c b a c a c k n b a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.3.2) 时,其行列式det ( ) 0 ? A ,即A是非奇异的. 由于方程组中零系数占多数,因此在进行顺序消元过程 时有必要修改并简化消元计算的公式,免去大量不必要的计 算.作A的Crout 分解,并把d 组作为 A 的最后一列统一处 理,可得 河南农业大学《数值计算方法》教案 汪松玉 - 67 - 1 2 2 3 , n n l a l a a l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L 1 2 1 1 1 , 1 n u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? U 1 2 n y y y ? ? ? ? ? ? ? ? ? ? ? ? ? ? y 其中 1 1 0 0 , 1,2, , ) / , ( 1,2, , 1) 0. i i i i i i i i i i i i l b a u y d y a l i n u c l i n u y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.3.3) 解方程 ? Ux y ,得: 河南农业大学《数值计算方法》教案 汪松玉 - 68 - 1, 1, 2, 2,1. n n i i i i x y x y u x i n n ? ? ? ? ? ? ? ? ? ? ? (5.3.4) 上述计算过程也称为求解三对角线性方程组的追赶法, 它由两组递推公式组成:其消元过程是按 i 从小到大的次序 来计算 , , i i i l y u ,(5.3.3)称为追的过程;回代过程求 i x 的次序 正相反,(5.3.4)称为赶的过程. 容易算出追赶法的运算量仅 为8n . 解三对角方程组(5.3.1)的追赶法,它分为"追"和"赶" 两个过程: (1) 追的过程(消元过程):按(5.3.3)的顺序计算系数: 1 2 n l l l ? ? ? ? 、 1 2 n y y y ? ? ? ? 及121nuuu?????河南农业大学《数值计算方法》教案 汪松玉 - 69 - (2) 赶的过程(回代过程):按5.3.4)的逆序求出解: 1 1 n n x x x ? ? ? ? ? . 注 追赶法具有计算量少,方法简单,算法稳定等优点, 在计算机上实现时,只需用三个一维数组分别存放方程组 (5.3.1)的系数矩阵 A 的三条对角线上的元素{ i i i a b c , 两组工作单元保存{ }, { } i i u y . 例5.7 用追赶法求解三对角方程组 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 10 4 2 4 1 1 4 1 1 4 3 2 1 x x x . 解23123,1, n a a c c ? ? ? ? ? ? 10 , 4 , 2 , 4 3 2 1 3 2 1 ? ? ? ? ? ? d d d b b b . 河南农业大学《数值计算方法》教案 汪松玉 - 70 - (追) 4 15 1 , 4 56 0 1 15 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L 1 1 - 0 4 4 1 - , 15 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? U 1 2 6 5 3 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? y . (赶) 3 3 2 2 2 3 3, 6 5 4 15 3 2 x y x y u x ? ? ? ? ? ? ? ? , 1 1 1 2 1 2 1 4 2 1 x y u x ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 71 - 5.4 对称正定矩阵的平方根法和 LDLT 分解 在很多实际问题中,归结为求解方程组 b Ax ? ,其中 n n R A ? ? 恰好局有特殊性质,为对称正定阵,即A.计算机 上求解这种类型方程组的有效方法之一就是平方根法,也就 是利用对称定矩阵三角分解的一种分解法. 定义 设A为n阶实方阵,如果 A 满足条件 (1) A AT ? ; (2) 对任意非零向量 n R x ? ,则有 0 T x Ax ? ; 则称 A 为对称正定矩阵. 关于对称正定矩阵有以下性质. 性质 1 如果 , n n R A ? ? 为对称正定阵,则(1) A 非奇异矩阵,且1?A亦是对称正定阵; 河南农业大学《数值计算方法》教案 汪松玉 - 72 - (2) 记kA为A的顺序主子阵,kA亦是对称正定阵),,2,1(nk??;(3) A 的特征值 ) , , 2 , 1 ( 0 ) ( n i A i ? ? ? ? ; (4) A 的顺序主子式都大于零,即). , , 2 , 1 ( 0 ) det( n k Ak ? ? ? 性质 2 (1)设nnRA??为对称矩阵,且A特征值 ( ) 0 i A ? ? ( 1, , i ? ? ) n ,则A为对称正定矩阵. (2) 设nnRA??为对称矩阵且 det( k A ) ) , , 2 , 1 ( 0 n k ? ? ? 则A为对称正定矩阵. 当A为n阶对称正定矩阵时,则A可分解成一个下三角 阵及其转置矩阵的乘积,即T?ALL (5.4.1) 河南农业大学《数值计算方法》教案 汪松玉 - 73 - 其中 L 是下三角阵. 事实上,将A作LU 分解, 11 12 1 21 22 2 1 2 1 1 1 n n n n nn u u u l u u l l u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A LU , 并写出其分块形式,有(1,2, , ) k k k k n ? ? ? A L U ,其中 11 12 1 21 22 2 1 2 , k k k k k kk a a a a a a a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A 21 1 2 1 1 , 1 k k k l l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L 河南农业大学《数值计算方法》教案 汪松玉 - 74 - 11 12 1 22 2 k k k kk u u u u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? U 由于 A 的各阶主子行列式都大于零,所以有 1 det det det 0, ( 1,2, , ). k k k k ii i u k n ? ? ? ? ? ? ? ? A L U (5.4.2) 故有 0, ( 1,2, , ) ii u i n ? ? ? . 进一步将U 分解为 DR ,其中 河南农业大学《数值计算方法》教案 汪松玉 - 75 - 11 22 0 , 0 nn u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? D 1 12 11 11 2 22 1, 1, 1 1 1 1 n n n n n n u u u u u u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? R . 可以证明上述分解是唯一的.再由 T ? A A,可得 T T R DL = LDR , 故有 , T T ? ? R L L R,即T?ALDL (5.4.3) 其中 L 为单位下三角矩阵,D 为对角矩阵.称(5.4.3)为对称正 河南农业大学《数值计算方法》教案 汪松玉 - 76 - 定矩阵 A 的TLDL 分解. 若A是对称正定矩阵, 0, ( 1,2, , ) ii u i n ? ? ? ,将D进一 步分解为: ? D DD ,其中 11 22 nn u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? D 则(5.4.3)又可以表示成 ( )( ) T T T ? ? ? A LDDL LD LD LL 记?LLD ,得分解式 T ? A LL (5.4.4) 河南农业大学《数值计算方法》教案 汪松玉 - 77 - 这里 L 是下三角矩阵,并且对角元全是正数. 通常称(5.4.4)为对称正定矩阵 A 的乔累斯基(Cholesky) 分解,或LLT 分解. 当方程组(5.0.1)的系数矩阵 A 为对称正定矩阵时,可以 用乔累斯基分解法求解,即利用 A 的TLL 分解或 T LDL 分解 式来求解,分别又称为平方根法和 T LDL 分解法.用平方根法 解方程组的计算过程分两个步骤: (1) 计算分解式 T ? A LL .设11 11 21 1 21 22 22 2 1 2 n n T n n nn nn l l l l l l l l l l l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A LL 河南农业大学《数值计算方法》教案 汪松玉 - 78 - 根据矩阵乘法有 1 2 2 2 1 1 ,( 1,2, , ) k k k k ks ks kk s s a l l l k n ? ? ? ? ? ? ? ? ? ? 以及 1 1 1 , k k jk js ks j s k s jk k k s s a l l l l l l ? ? ? ? ? ? ? ? ( 1, 2, , ) j i i n ? ? ? ? 于是有 1 2 1 1 1 , ( ) , ( 1, 2, , ) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? k kk kk ks s k jk js ks s jk kk l a l a l l l j k k n l (5.4.5) 河南农业大学《数值计算方法》教案 汪松玉 - 79 - 为了便于记忆,我们仿照杜利特尔分解的记忆方法,把元素 设想成如下排列 11 11 ( ) a l 12 12 ( ) a l 13 13 ( ) a l 14 14 ( ) a l ? 21 21 ( ) a l 22 22 ( ) a l 23 23 ( ) a l 24 24 ( ) a l ? 31 31 ( ) a l 32 32 ( ) a l 33 33 ( ) a l 34 34 ( ) a l ? 41 41 ( ) a l 42 42 ( ) a l 43 43 ( ) a l 44 44 ( ) a l ? ? ? ? ? ① ② ③ ④ 其中 , ( , 1,2, , ) ij ji ij ji a a l l i j n ? ? ? ? ,由(5.4.5)式可知, L 的 对角元 ii l 等于 A 的对应元素 ii a 减去一个内积,然后开平方. 内积的每一项仍然是同行元素 ik l 与同列元素 ki l 的乘积,但需 注意的是 ik ki l l ? . L 的非对角元素 ( ) ji l j i ? 的计算完全与 LU 河南农业大学《数值计算方法》教案 汪松玉 - 80 - 分解中 L 的元的计算相同,只是这里的对角元的名称是 ii l . 运算次序是按列进行. (2) 依次求解两个三角形方程组 ? Ly b及T?Lxy,即 ① b Ly ? 求y;②yxLT ? 求x;③11()/ ,( 1,2, , ) i i i ik k ii k y b l y l i n ? ? ? ? ? ? ? ; ④ 1 1, ,2,1) n i i ki k ii k i x y l x l i n n ? ? ? ? ? ? ? ? Cholesky 分解算法 设有 n n R A ? ? 为对称正定阵, 下述算 法计算下三角阵 L 使TLL A ? 且计算 ) ( j i lij ? 覆盖 ij a . 对于 n i , , 2 , 1 ? ? (行格式) 河南农业大学《数值计算方法》教案 汪松玉 - 81 - (1) 对于 1 , , 1 ? ? i j ? , jj j k jk ik ij ij ij a a a a l a / ) ( 1 1 ? ? ? ? ? ? (2) 2 / 1 1 1 2 ) ( ? ? ? ? ? ? i k ik ii ii ii a a l a 本算法大约需要 6 / 2 n 次乘除法. 由对称正定阵 A 分解计算公式(2)可得 ? ? ? ? i k ik ii n i l a 1 2 ) , , 2 , 1 ( , ? 于是, ii n i ii ik a a l ? ? ? ? 1 2 max ) , , 2 , 1 ; , , 2 , 1 ( i k n i ? ? ? ? 上式说明在不选主元的平方根法中,矩阵 L 元素有界, 或者说在分解过程中产生 L 元素 ik l 的数量级不会增长,且 河南农业大学《数值计算方法》教案 汪松玉 - 82 - ) , , 1 ( 0 n i lii ? ? ? ,由此可以说明,平方根算法是一个数值稳 定的方法. 对于正定矩阵的三角分解 T LL A ? 计算最约 6 / 3 n 次乘除 法,大约为 LU 分解计算的一半. 此外,计算过程式稳定的. 事实上,由关系式 2 1 ? ? ? k k k ks s a l 得||? ij ii l a 说明 Cholesky 分解中的量 ij l 能够得到控制, 因此其计算过程 是稳定的. 例5.8 用平方根法求解方程组 河南农业大学《数值计算方法》教案 汪松玉 - 83 - ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 7 8 5 11 0 2 0 2 1 2 1 1 3 2 1 x x x 解设(1)1 (1)1 (2)2 (1)1 (2)1 (0)-2 (2)2 (0)-2 (11) 3 于是,有100110223????????????L.由 ? Ly b 解得 1 2 3 5, 3, 3 y y y ? ? ? , 河南农业大学《数值计算方法》教案 汪松玉 - 84 - 由T?Lxy解得 1 , 5 , 2 3 2 1 ? ? ? ? x x x . 为了避免开方运算,有时直接使用对称矩阵 A 的TLDL 分解来计算,令T?ALDL 1 21 1 21 2 2 1 2 1 1 1 1 1 1 n n n n n d l l l d l l l d ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 根据矩阵乘法可以求出 L 和D的元素,然后将方程组(5.0.1) 即T?LDL x b 转化为两个三角形方程组 ? Ly b, 1 T ? ? L x D y , 由前一方程解出 y ,代入后一方程便可解出 x .具体计算步骤 可参看有关书籍.
下载该文档
文档格式：PDF
更新时间：2014-08-08
下载次数：0
点击次数：1
河南农业大学《数值计算方法》教案 汪松玉 - 1 - 第五章 线性代数方程组的直接解法 教学目的 1. 掌握解线性方程组的高斯消去法、 高斯选主元素消去 法; 2. 掌握用直接三角分解法解线性方程组的方法; 3. 了解解对称正定矩阵线性方程组的平方根法与解三 对角线方程组的追赶法; 4. 掌握向量,矩阵范数,矩阵的条件数等概念及方程组 的扰动分析. 教学重点及难点 重点:1. 解线性方程组的高斯消去法、高斯选主元素 消去法; 河南农业大学《数值计算方法》教案 汪松玉 - 2 - 2. 直接三角分解法解线性方程组的方法; 3. 向量,矩阵范数,矩阵的条件数等概念及方 程组的扰动分析; 难点:方程组的扰动分析. 河南农业大学《数值计算方法》教案 汪松玉 - 3 - 线性代数是数值计算方法的基础,解线性代数方程组的 有效方法在计算数学和科学计算中具有特殊的地位和作用. 如弹性力学、电路分析、热传导和振动、以及社会科学及定 量分析商业经济中的各种问题等许多实际问题的求解最终 转化为代数问题来处理,即归结为解线性方程组,如插值公 式的建立,微分方程离散格式的构造等.n 阶线性方程组的一 般形式是 11 1 12 2 1 1 21 1 22 2 2 2 1 1 2 2 n n n n n n nn n n a x a x a x b a x a x a x b a x a x a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??? ? (5.0.1) 其矩阵形式表示为 ? Ax b 河南农业大学《数值计算方法》教案 汪松玉 - 4 - 其中 11 12 1 1 1 21 22 2 2 2 1 2 , , n n n n nn n n a a a x b a a a x b a a a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A x b 一般 b≠0.当系数矩阵 A 非奇异(即det A ≠0)时, 方程组(5.0.1) 有唯一解 x . 求解上述线性代数方程组的方法主要有两类: (a) 直接法:对于给定的方程组,在没有舍入误差的假 设下, 能在预定的运算次数内求得精确解.最基本的直接法是 Gauss 消去法, 重要的直接法全都受到 Gauss 消去法的启发. 缺点:计算代价高. 河南农业大学《数值计算方法》教案 汪松玉 - 5 - (b) 迭代法:基于一定的递推格式,产生逼近方程组精 确解的近似序列, 用某种极限过程去逐步逼近精确解的方法. 收敛性是其为迭代法的前提,此外,存在收敛速度与误差估 计问题. 本章介绍解线性方程组的直接法(主要是高斯消去法及 其各种变形), 即如果计算过程中没有舍入误差, 经过有限步 算术运算可求得方程组的精确解,但在实际计算中由于舍入 误差的存在和影响,这类方法也只能求得方程组的近似解. 其中选主元消去法,三角分解法是解低阶稠密矩阵方程组的 有效方法. 河南农业大学《数值计算方法》教案 汪松玉 - 6 - 5.1 高斯消去法 5.1.1 高斯消去法的基本方法 高斯消去法是一种古老的方法,其基本思想是通过对方 程组(5.0.1)作下列运算(变换): (1) 交换方程组中任意两个方程的次序; (2) 方程组中任何一个方程乘以某一个非零常数; (3) 方程组中任何一个方程减去另一个方程的常数倍. 在系数矩阵 A 非奇异的假设下,可以得到一个与(5.0.1)等价 (或者说是具有相同解)的三角形方程组 河南农业大学《数值计算方法》教案 汪松玉 - 7 - (1) (1) (1) (1) (1) 11 1 12 2 1, 1 2 1, 1 (2) (2) (2) (2) 22 2 2, 1 2 2, 2 ( 1) ( 1) ( 1) 1, 1 2 1, 1 ( ) , n n n n n n n n n n n n n n n n n n n a x a x a x a x b a x a x a x b a x a x b a x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ( ) n n b ? ? ? ? ? ? ? ? ? ? (5.1.1) 这个过程称为消元. 显然,只要全部 ( ) , 0 k k k a ? (k=n, n-1,?,1),就可依次求出 1 1 , , , x x x n n ? ? , 这种按未知量编号从大到小的次序求解上三角 形方程组的过程称为回代. 这种将 n 阶线性方程组(5.0.1)化简为三角形方程组 (5.1.1)之后再回代求解的方法称之为高斯(Gauss)消去法. 河南农业大学《数值计算方法》教案 汪松玉 - 8 - 下面以(5.0.1)为例说明高斯消去法的步骤: 第一步 若11 0 a ? ,令1111 , ( 2,3, , ) i i a l i n a ? ? ? ,用1il?乘第 一个方程分别加到第i ( 2,3, , ) i n ? ? 个方程上去, 从而消去第 i ( 2,3, , ) i n ? ? 个方程的首项,方程组(5.0.1)变为下列形式: (1) (1) (1) (1) 11 1 12 2 1 1 (2) (2) (2) 22 2 2 2 (2) (2) (2) 2 2 n n n n n nn n n a x a x a x b a x a x b a x a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ? (5.1.2) 其中, 1 (1) 1 , 1,2, , j j a a j n ? ? ? , (1) 1 1 b b ? (2) 1 1 , , 2,3, , ij ij i j a a l a i j n ? ? ? ? , 河南农业大学《数值计算方法》教案 汪松玉 - 9 - (2) (1) 1 1 , 2,3, , i i i b b l b i n ? ? ? ? 第二步 若(2) 22 0 a ? ,令(2) 2 2 (2) 22 , ( 3,4, , ) i i a l i n a ? ? ? ,用1il?乘第二个方程分别加到第i ( 3,4, , ) i n ? ? 个方程上去, 从而消去 (2) 2 , ( 3,4, , ) i a i n ? ? . 一般地,设第 1 k ? 步后方程组化为如下形式 河南农业大学《数值计算方法》教案 汪松玉 - 10 - (1) (1) (1) (1) 11 1 12 2 1 1 (2) (2) (2) 22 2 2 2 n n n n k k k kk k kn n k k k k nk k nn n n a x a x a x b a x a x b a x a x b a x a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ? ???? ? (5.1.3) 若()0kkk a ? ,令()(),(1, 2, , ) k ik ik k kk a l i k k n a ? ? ? ? ? ,用1il?乘第k 个 方程分别加到第i 个(1, 2 , ) i k k n ? ? ? ? 方程上去,得到如下 方程组 河南农业大学《数值计算方法》教案 汪松玉 - 11 - (1) (1) (1) (1) 11 1 12 2 1 1 (2) (2) (2) 22 2 2 2 ( 1) ( 1) ( 1) 1, 1 1 1, 1 ( 1) ( 1) ( 1) , 1 1 n n n n k k k k k k k n n k k k k n k k nn n n a x a x a x b a x a x b a x a x b a x a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ? ???? ? (5.1.4) 其中 ( 1) , ( , 1, 2, , ) k k k ij ij ik kj a a l a i j k k n ? ? ? ? ? ? ? ( 1) , ( 1, 2, , ) k k k i i ik k b b l b i k k n ? ? ? ? ? ? ? 经过 1 n ? 步消元后,最后从(5.1.4)可得到一个上三角形方程 组 河南农业大学《数值计算方法》教案 汪松玉 - 12 - (1) (1) (1) (1) 11 1 12 2 1 1 (2) (2) (2) 22 2 2 2 ( ) ( ) n n n n k k k kk k kn n k n n nn n n a x a x a x b a x a x b a x a x b a x b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ? ???? (5.1.5) 现考虑回代过程的一般计算公式. 当()0nnn a ? 时,逐步回代得原方程组的解 河南农业大学《数值计算方法》教案 汪松玉 - 13 - ( ) ( ) 1 / , ( )/ , ( 1, 2, ,1) n nn k kj kk n n n n k k k k j j k x b a x b a x a k n n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.1.6) 因为高斯消去法的消元过程,实际上是对方程组的系数 进行运算,与未知量毫无关系,在实质上是对方程组(5.0.1) 的增广矩阵[A, b]进行初等行变换, 直至将系数矩阵化为上三 角形. 为了将高斯消去法编制成实用的计算机程序,就必须对 消元过程的具体计算步骤作出一些规定, 使之有统一的规律. 由不同的规定原则便产生了不同形式的算法.例如顺序高斯 消去法、列主元高斯消去法以及全主元高斯消去法,这些都 河南农业大学《数值计算方法》教案 汪松玉 - 14 - 是高斯消去法中最典型且常用的算法.其共同特点为: 它们的 主要步骤都是由消元和回代两大过程组成,主要区别在于消 元过程的计算步骤有各自的规范,而这些规范是与程序的简 明通用以及计算的稳定性有关的. 5.1.2. 顺序高斯消去法 顺序高斯消去法是按方程和未知量的自然顺序进行消元 的,它是各种高斯消去法的基础,下面举例说明具体计算步 骤.为方便起见我们约定用记号②- * ① k 表示第②行减去第 ①行的 k 倍,用①? ②表示交换第①与第②行,其余类推. 例5.1 用顺序高斯消去法解方程组 河南农业大学《数值计算方法》教案 汪松玉 - 15 - 1 2 3 1 2 3 1 2 3 2 2 3 3, 4 7 7 1, 2 4 5 7, x x x x x x x x x ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? (5.1.7) 解 先作消元过程.增广矩阵为 [A,b]= 2 2 3 3 4 7 7 1 2 4 5 7 ? ? ? ? ? ? ? ? ? ? ? ? 1 2 ? ? ???? ? 第 列消元 ② ① ③+① 2 2 3 3 0 3 1 5 0 6 8 4 ? ? ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 16 - 2 2 ? ? ???? ? 第 列消元 ③ ② 2 2 3 3 0 3 1 5 0 0 6 6 ? ? ? ? ? ? ? ? ? ? ? 由此得上三角形方程组 1 2 3 2 3 3 2 2 3 3, 3 5, 6 6, x x x x x x ? ? ? ? ? ? ? ? ? ? ? ? (5.1.8) 再作回代过程.由(5.1.8)第三个方程求得 3 1 x ? ,代入第二个方 程得 2 3 ( 5 )/3 2 x x ? ? ? ? ? .将23,xx代入第一个方 程可得 1 2 3 (3 2 2 )/ 2 2 x x x ? ? ? ? ,所以方程组的解为 1 2 3 2, 2, 1 x x x ? ? ? ? . 上例的求解用的是最基本的顺序高斯消去法,它的消元过程 河南农业大学《数值计算方法》教案 汪松玉 - 17 - 遵守如下规则: (1) 依从左到右、自上而下的次序将主对角元下方的元 素化为零; (2) 不作行交换,也不作将某行乘以一个非零数的变换; (3) 当进行第 k 列消元时(即使主对角元 kk a 下方化零), 将kk a 下面各行分别减去第 k 行的适当倍数,不作其他变换. 例5.1 主要是用手算表示出指定算法的计算步骤,对数 据形式未作具体要求,其实在计算机里上述整数或分数均为 浮点数,计算繁简程度相同. 现在,我们将例 5.1 中使用的基本高斯消去法推广到解 一般的 n n ? 阶线性方程组(5.0.1).设方程组的系数矩阵非奇 异,其增广矩阵记为 河南农业大学《数值计算方法》教案 汪松玉 - 18 - 11 12 1 21 22 2 1 2 (1) (1) (1) (1) 1 (1) (1) (1) (1) (1) (1) 2 (1) (1) (1) (1) [ , ] n n n n nn n a a a b a a a b b a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b (5.1.9) 其中, 1 (1) 1 , j j a a ? (1) , i i b b ? ( , 1,2, , ) i j n ? ? . 高斯消元法的消元过程由 1 n ? 步组成: 第一步 若(1) 11 0 a ? ,令(1) 1 1 (1) 11 , ( 2,3, , ) i i a l i n a ? ? ? ,用1il?乘(5.1.9)第一行分别加到第 i ( 2,3, , ) i n ? ? 行上去,从而消去 (1) (1) (1) 21 31 1 , , , n a a a ? ,得到 河南农业大学《数值计算方法》教案 汪松玉 - 19 - 1 11 12 1 2 22 2 (1) (1) (2) (1) (2) (2) (2) 2 (2) (2) (2) (2) (2) 0 [ , ] 0 n n nn n n b a a a b a a A b a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.1.10) 其中, (2) (1) (1) 1 1 , ( , 2,3, , ) ij ij i j a a l a i j n ? ? ? ? , (2) (1) (1) 1 1 ,( 2,3, , ) i i i b b l b i n ? ? ? ? 第二步 若(2) 22 0 a ? ,令(2) 2 2 (2) 22 , ( 3,4, , ) i i a l i n a ? ? ? ,用1il?乘(5.1.10)第二行并分别加到第i ( 3,4, , ) i n ? ? 行上去,从而消 去42 2 (2) (2) (2) 32 , , , n a a a ? .依次下去,假设进行了 1 k ? 步,得到 河南农业大学《数值计算方法》教案 汪松玉 - 20 - 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 , 1 1, 1, 1 1, 1 , 1 [ , ] k n k k k n k k k k k k kk k k kn k k k k k k k k k k n k k k k k nk n k nn n a a a a a b a a a a b a a a b a a a b a a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b ? ? ? ? ? (5.1.11) 第k步假设 ( ) 0 k kk a ? , ( ) ( ) , ( 1, 2 , ) k ik ik k kk a l i k k n a ? ? ? ? ? , 按照上述 方法消去(5.1.11)式第 k 列的 ( ) ( ) 1 , , k k k k nk a a ? ? ,得到 河南农业大学《数值计算方法》教案 汪松玉 - 21 - 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 ( 1) ( 1) , 1 ( 1) ( 1) ( 1) 1, 1 1, 1 ( 1) ( 1) ( 1) , 1 [ , ] k n k k k n k k k k k k kk k k kn k k k k k k k n k k k k n k nn n a a a a a b a a a a b a a a b a a b a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b ? ? ? ? ? ? (5.1.12) 其中 ( 1) ( 1) , , 1, 2, , , 1, 2, , k k k ij ij ik kj k k k i i ik k a a l a i j k k n b b l b i k k n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.1.13) (5.1.13)式是消元过程的一般计算公式, 式中作为分母的元素 ( ) k kk a 称为第 k 步的主元素(简称主元), ( 1, 2 , ) ik l i k k n ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 22 - 称为乘子.若()0kkk a ? ,则()()1,,kkkknk a a ? ? 中至少有一个元素,例如()krk a 不为零, (否则, 方程组(5.0.1)的系数矩阵 A 奇异).这样, 我们可取 ( ) k rk a 为主元,然后交换矩阵 ( ) ( ) , k k ? ? ? ? A b 的第k 行与第 r 行,把它交换到( , ) k k 的位置上. 经过 1 n ? 步消元后, 最后从(5.1.12)可得到一个上梯形矩 阵 河南农业大学《数值计算方法》教案 汪松玉 - 23 - 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 , 1 ( 1) ( 1) ( 1) 1, 1 1, 1 ( ) ( ) [ , ] k n k k k n n n k k k k kk k k kn k k k k k k k n k n n nn n a a a a a b a a a a b a a a b a a b a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b (5.1.14) 这里,我们假设整个消元过程中没有进行矩阵的行交换. ( ) n A 是一个上三角矩阵,与()()[,]nnAb相应的上三角方程组(5.1.5) 和(5.0.1)同解. 使用高斯消去法求解n 阶线性方程组总共需要乘除法运 算次数为 3 2 1 ( 3 - ) 3 n n n ? 次, 减法运算次数为 3 2 1 (2 3 -5 ) 6 n n n ? 次. 河南农业大学《数值计算方法》教案 汪松玉 - 24 - 由于一般计算机中乘除运算时间远远超过加减运算时间,所 以在估计一种计算方法的运算量时, 往往只估计乘除的次数. 又因为n 很大时, 1 p p n n ? ,所以在比较计算量的大小时, 常常只比较n 的最高次项.因此,我们可以说高斯消去法的乘 除次数为 3 1 3 n 数量级. 从上面讨论不难看出, 高斯消去法的消去过程把方程组 化成上三角形方程组, 必须通过回代才能得到解.如果在消元 的过程稍作改变,即将对角线上方的元素也化为零,则可使 方程(5.0.1)化简成为每个方程只含有一个未知量的方程组, 因此求解时无需回代的步骤.这种消去法称为高斯-约当 (Gauss-Jordan)消去法,也称为无回代的高斯消去法. 但是, 经过计算我们发现,这种方法所需要的计算量太大,其中乘 河南农业大学《数值计算方法》教案 汪松玉 - 25 - 除法次数为 3 2 1 1 2 2 n n n ? ? ,即为 3 1 2 n 数量级,超过高斯消去法, 所以一般来说, 求解线性方程组时采用高斯-约当消去法是不 明智的.但是,用它来求解矩阵的逆矩阵时却很方便,有兴趣 的读者请参阅其他相关文献. 5.1.3 选主元高斯消去法 上一小节我们介绍的是顺序高斯消去法,其消元过程能 进行到底的充要条件是主元素 ( ) 0 k kk a ? ,( 1,2, , ) k n ? ? (其中 (1) 11 11 a a ? ), 因为它们在计算公式中是被作为除数的.但即使在 消元过程中主元素 ( ) k kk a ( 1,2, , ) k n ? ? 全不为 0, 也不能保证能 够避免在计算过程中产生巨大的误差.因为在计算过程中, 若主元素()kkk a 的绝对值相对其它元素(如与同列的元素河南农业大学《数值计算方法》教案 汪松玉 - 26 - ( ) ( ) 1 , , k k k k nk a a ? ? 相比)过小, 在这种情况下, 将导致舍入误差的扩 散,最后使计算结果不可靠,甚至消元过程无法进行,下面 举例说明. 例5.2 考察方程组-5 1 2 1 2 10 1 2 x x x x ? ? ? ? ? ? ? , 其精确解为1100000/99999 x ? , 2 99998/99999. x ? 若使用四位浮点十进 制运算进行消元,得到 -5 1 2 -5 -5 2 10 + =1 10 =10 x x x ? ? ? ? ? 利用回代过程得方程组的解为: 1 2 0, 1 x x ? ? .若用消去法求 河南农业大学《数值计算方法》教案 汪松玉 - 27 - 解12-5 1 2 2 10 1 x x x x ? ? ? ? ? ? ? 则可得 1 1 x ? , 2 1 x ? . 显然,第一个结果是严重失真的,而第二个结果则比较 精确.究其原因,就是小主元的出现.用它来作除数带入了大 的舍入误差,再经传播,误差变得更大,造成计算结果与真 解相差甚远,因此第一种算法是不稳定的. 在顺序高斯消去法中,出现零主元将导致求解失败;而 出现小主元又会导致计算数值解精度的严重失真.为避免出 现不稳定的算法,我们在使用消去法求解方程组(5.0.1)时, 在消元过程中引入逐次选主元的思想,其基本思路是每次消 元之前在系数矩阵中按一定的范围选取绝对值最大的元素 作为主元素,以便减少舍入误差的影响,这种方法称为选主 元高斯消去法, 是对顺序高斯消去法的一种改进.常用的选主 河南农业大学《数值计算方法》教案 汪松玉 - 28 - 元高斯消去法是列主元高斯消去法、行主元消去法和全主元 高斯消去法. 设方程组(5.0.1)的系数矩阵A非奇异,进行第 k-1 次消 元后,增广矩阵为 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 , 1 1, 1, 1 1, 1 , 1 [ , ] k n k k k n n n k k k k kk k k kn k k k k k k k k k k n k k k k k nk n k nn n a a a a a b a a a a b a a a b a a a b a a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A ? ? ? ? ? ? ? ? ? ? ? b ? ? ? ? ? (5.1.15) 列主元高斯消去法:是在第 k 次(k=1,2,… ,n-1)消元时, 从系数矩阵第 k 列中 ( 1) 1, k k k a ? ? 以下的各元素中寻找绝对值最大 河南农业大学《数值计算方法》教案 汪松玉 - 29 - 的元素作为主元素(称为列主元素),增广矩阵中主元所在的 行称为主行, 主元所在的列称为主列.我们的工作就是要确定 r,使()krk a 为()kkk a 1 , , k k k k nk a a ? ? 中之最大者. 若r>k 则交换增广矩阵的第 k 行和第 r 行,使()krk a 代替 ( ) k kk a 成为主元素,然后进行消元.其余步骤与顺序高斯消去法 相同.由(5.1.15)可以看出,只要 det A≠0,就必存在 ( , r r k ? 1, , ) k n ? ? ,使得 ( ) k r k a ≠0,因此整个计算过程一定能进行到 底. 行主元高斯消去法:是在第 k 次(k=1,2,… ,n-1)消元时, 若从 1 2 , , , k k k kk kk kn a a a ? ? ? 中选取绝对值最大的元素作为主元素 (称为行主元素),即若 河南农业大学《数值计算方法》教案 汪松玉 - 30 - ( ) ( ) max k k k kj kj k j n a a ? ? ? 则选取 ( ) k k kj a 作主元,并且在进行第 k 次消元前交换增广矩阵 的第 k 列和 k j 列(必须记录这种交换信息,以便整理方程组 的解时用).经这样修改的高斯消去法称为行主元高斯消去法. 全主元高斯消去法也简称全主元消去法,它与列主元消 去法的不同之处,是在第 k 次(k=1,2, ?,n-1)消元时,从系 数矩阵(见(5.1.15))的右下角(n-k+1)阶子矩阵 ( ) ( ) ( ) ( ) k k kk kn k k nk nn a a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? 中,选取绝对值最大的元素作为主元素,如果它位于第 r 行 河南农业大学《数值计算方法》教案 汪松玉 - 31 - 第s列, 则通过交换 k, r 两行及交换 k, s 两列, 可以使主元素 位于 ( ) k kk a 的位置, 然后进行消元计算. 由于作列的交换改变 了方程中未知量的次序, 因此回代过程实际上是按未知量调 换后的编号顺序求解的. 下面举例说明计算步骤. 例5.3 用列主元高斯消去法求解下面的方程组(用三位 有效数字计算) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 9 5 . 0 4 5 5 2 5 . 0 5 5 . 0 2 01 . 0 3 2 1 3 2 1 3 2 1 x x x x x x x x x 解 用方程组的增广矩阵进行消元,每次选出的主元用 小方框标出 河南农业大学《数值计算方法》教案 汪松玉 - 32 - 0.01 2 0.5 5 [ , ] 1 0.5 2 5 9 5 4 0.5 A b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ③ ① ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 5 5 9 5 . 0 2 01 . 0 2 5 . 0 1 5 . 0 4 5 (第1列消元, 在0.01,-1,5 中取绝对值最大者 5 作主元, 交换①、③行使 5 换到 11 a 的位置, 再作消元计算.) 河南农业大学《数值计算方法》教案 汪松玉 - 33 - ? ? ? ? ? ? ? ? ? ? ? 500 1 5 1 ① ③ ① ② 5 4 0.5 9 0 1.30 2.10 6.80 5.02 0 2.01 0.501 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?③ ② ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 80 . 6 02 . 5 9 10 . 2 30 . 1 0 501 . 0 01 . 2 0 5 . 0 4 5 (第2列消元,在-1.30 和2.01 中取绝对值最大者 2.01 作 主元,交换②、③行使 2.01 换到 22 a 的位置,再作消元计算.) 河南农业大学《数值计算方法》教案 汪松玉 - 34 - ③ ② ? ? ? ? ?????? ? 130 201 . . 5 4 0.5 9 0 2.01 0.501 5.02 3.55 0 0 1.78 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 以上完成了消元过程, 得到一个上三角形方程组.再作回 代可求得 3 2 1 1.99, 2.00, 0.00100 x x x ? ? ? ? 事实上这个方程组的准确解是 2 , 2 , 0 3 2 1 ? ? ? ? x x x ,可见 用列主元消去法求得的近似解是令人满意的. 作为比较,用顺序高斯消去法重解例 5.3,对方程组的 增广矩阵作顺序消元得 河南农业大学《数值计算方法》教案 汪松玉 - 35 - 0.01 2 0.5 5 1 0.5 2 5 5 4 0.5 9 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 100 500 0.01 2 0.5 5 0 200 48.0 495 0 1000 251 2510 ? ? ? ? ? ? ? ? ? ? ???? ? ? ? ? ? ? ? ? ? ? ② ① ③ ① 5.00 0.01 2 0.5 5 0 200 48.0 495 0 0 11.0 30.0 ? ? ? ? ? ? ? ? ????? ? ? ? ? ? ? ? ? ? ③ ② 再作回代求得 3 2 1 2.73, 1.82 , 1.00 x x x ? ? ? ? ,这样的结果误 差太大,所以是无效的. 河南农业大学《数值计算方法》教案 汪松玉 - 36 - 例5.4 用全主元高斯消去法求解下面的方程组(用四位 有效数字计算) 1 2 3 1 2 3 1 2 3 12 3 3 15, 18 3 15, 6, x x x x x x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 它的准确解为 1 2 3 1, 2, 3 x x x ? ? ? . 解 用方程组的增广矩阵进行消元,每次选出的主元用 小方框标出 12 3 3 15 [ , ] 18 3 1 15 1 1 1 6 A b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 37 - ? ???? ? ① ② 18 3 1 15 12 3 3 15 1 1 1 6 ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 ?????? ? 第 列消元 18 3 1 15 0 1 2.333 5.000 0 1.167 0.944 5.167 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?????? ? 交换②、③列18 1 3 15 0 2.333 1 5.000 0 0.944 1.167 5.167 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 38 - ?????? ? 第2列消元 18 1 3 15 0 2.333 1 5.000 0 0 1.572 3.144 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 消元完成.由于交换过②、 ③两列, 所以改变了方程中 2 x 和3x的次序,最后得到与原方程组等价的上三角形方程组 1 3 2 3 2 2 18 3 15 2.333 5.000 1.572 3.144 x x x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? 回代即可逐步解出: 2 3 1 2.000, 3.000, 1.000. x x x ? ? ? 实际计算和理论分析表明, 列主元高斯消去法具有较好 的数值稳定性,是目前求解中小型稠密线性方程组最常用的 河南农业大学《数值计算方法》教案 汪松玉 - 39 - 方法之一. 在实际计算时,有些具有特殊系数矩阵的方程组,按顺 序消去法,若能保证主元的绝对值最大,或者设舍入误差能 得到控制,则不必选主元.这类矩阵有对称的严格对角占优 阵,对称正定阵等. 定义 5.1 所谓严格对角占优阵,是指 A 的元素满足如下 条件的矩阵: 1 1,2 , . n ii ij j j i a a i n ? ? ? ? ? ? 可以证明,严格对角占优矩阵必定是非奇异阵,而且经 过一步高斯消去法后仍然保持对角占优的优势. 河南农业大学《数值计算方法》教案 汪松玉 - 40 - 5.2 高斯消去法的矩阵形式 和矩阵的 LU 分解 顺序高斯消去法是对方程组(5.0.1)的增广矩阵 11 12 1 21 22 2 1 2 (1) (1) (1) (1) 1 (1) (1) (1) (1) (1) (1) 2 (1) (1) (1) (1) [ , ] n n n n nn n a a a b a a a b b a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A b 进行一系列初等变换来实现的,根据线性代数知识,初等变 换与初等矩阵有密切关系的.所以,顺序消元法有如下形式: 第一步 若11 0 a ? ,令1111 ,( 2,3,4 , ) i i a l i n a ? ? ? ,组成初 河南农业大学《数值计算方法》教案 汪松玉 - 41 - 等下三角阵 21 1 1 1 0 1 0 1 n l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L (5.2.1) 则111 12 1 2 22 2 (1) (1) (2) (1) (2) (2) (2) 2 (1) (1) (2) (2) 1 (2) (2) (2) 0 0 n n nn n n b a a a b a a a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L A b A b 即用 1 i l ? 乘(5.1.9)的第一行分别加到第 ( 2,3, , ) i i n ? ? 行上去, 相当于用 1 L 左乘增广矩阵. 河南农业大学《数值计算方法》教案 汪松玉 - 42 - 第二步 若(2) 22 0 a ? ,令(2) 2 2 (2) 22 , ( 3,4 , ) i i a l i n a ? ? ? ,组成 2 32 2 1 0 0 1 0 1 0 1 n l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L (5.2.2) 用1il?乘(5.1.10)第二行并分别加到第 i( 3,4 , ) i n ? ? 行上去, 相当于用 2 L 左乘 (1) (1) 1[ , ] L A b . 整个消元过程是从左到右、自上而下地将主对角元下 方的 21 31 1 32 42 2 1 n n n n a a a a a a a ? ? ? ? 化为 0,这就相当于 河南农业大学《数值计算方法》教案 汪松玉 - 43 - 用形如(5.2.1),(5.2.2)的初等矩阵 1 2 3 1 , , , , n? ? L L L L 依次左乘 (1) (1) [ , ] A b .经过 1 n ? 步消元后,最后得到 (1) (1) 1 2 2 1[ , ] n n ? ? ? L L L L A b 11 12 1 1 (1) (1) (1) (1) (1) (1) 1, 1 1 (2) (2) (2) (2) (2) 22 2 2, 1 2 2 , 1 ( 1) ( 1) ( 1) 1, 1 1, 1 ( ) ( ) k n k k k n k k k k kk k k kn k k k k k k k n k n n nn n a a a a a b a a a a b a a a b a a b a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ( ) ( ) [ , ] n n ? A b 河南农业大学《数值计算方法》教案 汪松玉 - 44 - 于是有 (1) ( ) 1 2 2 1 (1) ( ) 1 2 2 1 , . n n n n n n ? ? ? ? ? ? ? ? L L L L A A L L L L b b (5.2.3) 注意到 ( 1,2, , 1) i i n ? ? ? L 是可逆的且为主对角元全为 1 的下 三角阵,而有限个单位下三角阵的乘积仍是单位下三角阵, 故可令 1 1 2 1 ( ) n ? ? ? ? L L L L , 11 12 1 (1) (1) (1) (1) 1 (2) (2) (2) 2 22 2 ( ) ( ) , n n n n nn n b a a a b a a a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? U y 则得: 1 (1) (1) ? ? L A b U y ,或1(1) 1 (1) , ? ? ? ? L A U L b y .高斯消 河南农业大学《数值计算方法》教案 汪松玉 - 45 - 去法的回代过程就是求解 ? Ux y . 又由 1 (1) ? ? L A U 可得: (1) ? A LU ,这说明,在()0kkk a ? 的 条件下,通过高斯消去法可以把方程组(5.0.1)的系数矩阵 A 分解为单位下三角阵和上三角阵的乘积, 称为 A 的LU 分解. 一旦实现了矩阵 A 的LU 分解,那么求解方程组 A ? x b 的问题就等价于求解两个三角形方程组: (1) L ? y b ,求y;(2) U ? x y ,求x.定义 5.2 设A为n阶矩阵,若A有分解式 ? A LU ,其中L为下三角矩阵,U 为上三角矩阵,则称之为 A 的LU 分 解或三角分解.特别若 L 是单位下三角矩阵, 则该分解式称为 杜利特尔(Doolittle)分解,如 河南农业大学《数值计算方法》教案 汪松玉 - 46 - 11 12 1 21 22 2 1 2 1 1 1 n n n n nn u u u l u u l l u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A 若U 是单位上三角矩阵,则称之为克劳特(Crout)分解,如A=????????????????????????1112112 2 1 22 21 11 ? ? ? ? ? ? ? n n nn n n u u u l l l l l l 的分解式就是克劳特分解式. 一个矩阵不一定存在三角分解式,当存在时一般不唯 河南农业大学《数值计算方法》教案 汪松玉 - 47 - 一.如2445?????????A,0112???????B,则B不存在三角分解(请自 行用反证法证明), A 的三角分解式不唯一,容易验证 2 0 1 2 4 1 0 3 ? ? ? ? ? ? ? ? ? ? ? ? ? ? A 是一个三角分解式,任取可逆的对角矩 阵2005???????D,可得新的三角分解式 1 2 0 1 2 4 1 0 3 ? ? ?? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ?? ? A D D 2 0 2 0 1/ 2 0 1 2 4 1 0 5 0 1/5 0 3 ? ?? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? 河南农业大学《数值计算方法》教案 汪松玉 - 48 - ? ? ? ? ? ? ? ? ? ? ? ? ? ? 5 / 3 0 1 2 / 1 5 8 0 4 那么,一个矩阵在什么条件下存在唯一的三角分解式呢? 定理 5.1 n 阶矩阵 A 存在唯一的杜利特尔分解和克劳特 分解的充要条件是 A 的顺序主子矩阵 k A (k=1,2,? ,n)非奇异. 推论 若定理 5.1 的条件满足,则A存在唯一的 LDU 分解式 ? A LDU ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 1 1 1 1 1 2 1 12 2 1 2 1 21 ? ? ? ? ? ? ? n n n n n u u u d d d l l l (5.2.4) 其中 L 为单位下三角阵,U 为单位上三角阵, D 为对角阵. 河南农业大学《数值计算方法》教案 汪松玉 - 49 - 反之也真. 设方程组(5.0.1)的系数矩阵 A 非奇异,并有三角分解 ? A LU .求解线性方程组的三角分解法的基本思想是:可以 先把 A 作LU 分解,存储 L 和U ,方程组就化为 x ? LU b 从而使问题转化为求解两个简单的三角形方程组 ? Ly b, ? Ux y (5.2.5) 再根据各个方程组不同的右端向量分别求解(5.2.5)的两个三 角形方程组即可. 下面首先介绍杜利特尔分解法.其步骤为: 第一步 三角分解 对k=1,2,…,n,计算U 和L的元素.为了在计算机上求解, 需找出一般计算公式.设A的杜利特尔分解式 LU A ? 为 河南农业大学《数值计算方法》教案 汪松玉 - 50 - ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? nn n n n n nn n n n n u u u u u u l l l a a a a a a a a a ? ? ? ? ? ? ? ? ? ? ? ? ? 2 22 1 12 11 2 1 21 2 1 2 22 21 1 12 11 1 1 1 (5.2.6) 由矩阵乘法知 A 的每个元素等于 L 的某一行乘U 的某一列, 若能按一定次序每次取出一个这样的等式,使其中只含一个 未求的元素,便可一步步直接算出 L 和U 的元素. 首先容易看出 A 和U 的第一行元素对应相等,即11,(1,2, , ) j j a u j n ? ? ? (5.2.7) 现在 11 u 已求, 故由 A 的第一列元素 11 1 1 u l a i i ? 可求出 L 的第一 列元素 1 1 11 / , ( 2,3, , ) i i l a u i n ? ? ? (5.2.8) 河南农业大学《数值计算方法》教案 汪松玉 - 51 - 是否可以如此继续,依次交替求出U 的第二行、 L 的第 二列, U 的第三行、L 的第三列…?不妨设U 的前 1 ? k 行和 L 的前 1 ? k 列已求, 欲求U 的第 k 行元素 ) ( k j ukj ? , 为此考虑 A 中对应元素 kj a , kj a 为L的第 k 行与U 的第 j 列对应位置元素 乘积之和,即kj a = 1 , 1 [ 1 0 0 ] k k k l l ? ? ? 1 1, j k j kj nj u u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 52 - = 1 1 , 1 1, k j k k k j kj l u l u u ? ? ? ? ? ? 因11,,?kkkll?及jkjuu11,,??分别为 L 中前 1 ? k 列、 U 中前 1 ? k 行元素,故上式唯一的未知数便可求出,即11,(,1, , ) k k j k j k s s j s u a l u j k k n ? ? ? ? ? ? ? ? (5.2.9) 再求 L 的第 k 列元素 ) ( k i lik ? ,考虑 A 中对应元素 ) ( k i aik ? ,与上面类似地可得 1 1 ( ) , ( 1, 2, , ) k ik is sk s ik k k a l u l i k k n u ? ? ? ? ? ? ? ? ? (5.2.10) 编程时对 n k , , 2 , 1 ? ? 循环使用(5.2.9)、(5.2.10)便可依次交替 计算U 的第 k 行、 L 的第 k ( n k , , 2 , 1 ? ? )列,直至全部算出. 河南农业大学《数值计算方法》教案 汪松玉 - 53 - 计算公式为 1 1 1 1 , , 1, , ( ) / , 1, 2, , k k j k j k s s j s k ik ik i s sk k k s u a l u j k k n l a l u u i k k n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.2.11) 第二步 求解过程 (1) 解?Ly b.对k=1,2,…,n 计算 1 1 k k k k j j j y b l y ? ? ? ?? , (5.2.12) (2) 解?Ux y .对k=n,n-1,…,1 计算 河南农业大学《数值计算方法》教案 汪松玉 - 54 - 1 ( )/ n k k k j j k k j k x y u x u ? ? ? ? ? . (5.2.13) 为了便于记忆公式并说明计算U 的行和 L 的列的次序,列 出下面的紧凑格式表: 11 11) ( u a 12 12 ) ( u a 13 13 ) ( u a 14 14 ) ( u a ? ① 1 1) ( y b 21 21) ( l a 22 22 ) ( u a 23 23 ) ( u a 24 24 ) ( u a ? ③ 2 2 ) ( y b 31 31) ( l a 32 32 ) ( l a 33 33 ) ( u a 34 34 ) ( u a ? ⑤ 3 3 ) ( y b 41 41) ( l a 42 42 ) ( l a 43 43 ) ( l a 44 44 ) ( u a ? ⑦ 4 4 ) ( y b ? ? ? ? ? ② ④ ⑥ ⑧ 其中①,②,?表示计算U 的行和 L 的列的次序,计算过程是 河南农业大学《数值计算方法》教案 汪松玉 - 55 - 按照一行一列、二行二列…的次序进行. 由公式(5.2.7)–(5.2.10)可知,上表中U 的第一行等于 A 的第一行,其余元素 j k u 等于 A 的对应元素 j k a 减去左边 L 的 同一行元素与上边U 的同一列元素对应乘积之和;计算 L 的 元素与计算U 的元素方法类似, 只是最后要除以同一列U 的 对角元.例如, 21 21 11 24 24 21 14 43 43 41 13 42 23 33 / , , ( )/ l a u u a l u l a l u l u u ? ? ? ? ? ? 等等.解方程组时,右端常数列b 作为矩阵的最后一列,计算 y 与计算U 的方法类似. 例 用三角分解法解线性方程组 河南农业大学《数值计算方法》教案 汪松玉 - 56 - . 7 17 3 5 3 0 1 0 3 4 2 1 1 0 1 0 0 2 0 1 4 3 2 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? x x x x 解 列表 1 11 ? u 0 12 ? u 2 13 ? u 0 14 ? u 5 1 ? y ① 0 21 ? l 1 22 ? u 0 23 ? u 1 24 ? u 3 2 ? y ③ 1 31 ? l 2 32 ? l 2 33 ? u 1 34 ? u 6 3 ? y ⑤ 0 41 ? l 1 42 ? l 0 43 ? l 2 44 ? u 4 4 ? y ⑦ ② ④ ⑥ 11 11 12 12 13 13 14 14 1 1 1, 0, 2, 0, 5; u a u a u a u a y b ? ? ? ? ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 57 - ; 0 1 0 , 1 1 1 , 0 1 0 11 41 41 11 31 31 11 21 21 ? ? ? ? ? ? ? ? ? a a l a a l a a l , 1 0 0 1 12 21 22 22 ? ? ? ? ? ? ? u l a u , 0 2 0 0 13 21 23 23 ? ? ? ? ? ? ? u l a u , 1 0 0 1 14 21 24 24 ? ? ? ? ? ? ? u l a u ; 3 5 0 3 1 21 2 2 ? ? ? ? ? ? ? b l b y 32 31 12 42 41 12 32 42 22 22 2 1 0 1 0 0 2, 1; 1 1 a l u a l u l l u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? , 2 0 2 2 1 4 23 32 13 31 33 33 ? ? ? ? ? ? ? ? ? ? ? u l u l a u , 1 1 2 0 1 3 24 32 14 31 34 34 ? ? ? ? ? ? ? ? ? ? ? u l u l a u ; 6 3 2 5 1 17 2 32 1 31 3 3 ? ? ? ? ? ? ? ? ? ? ? y l y l b y 河南农业大学《数值计算方法》教案 汪松玉 - 58 - 43 41 13 42 23 43 33 0 0 2 1 0 ; 2 a l u l u l u ? ? ? ? ? ? ? ? ? ? 44 44 41 14 42 24 43 34 3 0 0 1 1 0 1 2, u a l u l u l u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 4 4 41 1 42 2 43 3 7 0 5 1 3 0 6 4. y b l y l y l y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 由表可知 1 0 0 0 1 0 2 0 5 0 1 0 0 0 1 0 1 3 , , 1 2 1 0 0 0 2 1 6 0 1 0 1 0 0 0 2 4 L U ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? y 由U ? x y ,即 河南农业大学《数值计算方法》教案 汪松玉 - 59 - 1 2 3 4 1 0 2 0 5 0 1 0 1 3 0 0 2 1 6 0 0 0 2 4 x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 可得 4 2 4 x ? ,于是 4 3 2 1 2, 6 2 2 2, 3 1 2 1, 5 2 2 1 x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? . 例5.5 用杜利特尔分解法求解方程组 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 2 3 4 2 4 9 16 10 8 27 64 44 16 81 256 190 x x x x x x x x x x x x x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? . 河南农业大学《数值计算方法》教案 汪松玉 - 60 - 解设?ALU ,按紧凑格式表计算 L ,U 以及方程组 b = y L 的解 y ,得(1)1 (2)2 (3)3 (4)4 (2)2 (1)1 (4)2 (9)6 (16)12 (10)8 (1)1 (8)3 (27)6 (64)24 (44)18 (1)1 (16)7 (81)6 (256)24 (190)24 即1123421102612 8 , , 1 3 1 6 24 18 1 7 6 1 0 24 24 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L U y . 再求解 ? Ux y 求得原方程组的解为 河南农业大学《数值计算方法》教案 汪松玉 - 61 - 4 3 2 1 1 , 1, 1, 1. x x x x ? ? ? ? ? ? 接下来我们介绍矩阵的 Crout 分解,并导出计算公式. 设A可分解为下三角阵和单位上三角阵的乘积,即11 12 1 21 22 2 1 2 n n n n nn a a a a a a a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A 11 12 1 21 22 2 1 2 1 1 1 n n n n nn l u u l l u l l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? , 河南农业大学《数值计算方法》教案 汪松玉 - 62 - 则对 1,2, , i n ? ? 有1111,(,1, , ) 1, 2, , ) i ji ji jk ki k i i j ik ik kk ii k l a l u j i i n u a l u u j i i n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.2.14) 仿照杜利特尔分解的记忆方法,我们给出下面的紧凑格式 表: 11 11 ( ) a l 12 12 ) ( u a 13 13 ) ( u a 14 14 ) ( u a ? ② 1 1) ( y b 21 21) ( l a 22 22 ( ) a l 23 23 ) ( u a 24 24 ) ( u a ? ④ 2 2 ) ( y b 31 31) ( l a 32 32 ) ( l a 33 33 ( ) a l 34 34 ) ( u a ? ⑥ 3 3 ) ( y b 41 41) ( l a 42 42 ) ( l a 43 43 ) ( l a 44 44 ( ) a l ? ⑧ 4 4 ) ( y b 河南农业大学《数值计算方法》教案 汪松玉 - 63 - ? ? ? ? ? ① ③ ⑤ ⑦ 其中①,②, ?表示计算 L 的列和U 的行的次序, 计算过程是 按照一列一行、二列二行?的次序进行. 由公式(5.2.14)可知,上表中 L 的第一列等于 A 的第一 列,其余元素 j k u 等于 A 的对应元素 j k a 减去左边 L 的同一行 元素与上边U 的同一列元素对应乘积之和; 计算U 的元素与 计算 L 的元素方法类似, 只是最后要除以同一行 L 的对角元. 解方程组时,右端常数列b 作为矩阵的最后一列,计算 y 与 计算U 的方法类似. 例5.6 用Crout 分解法求解例 5.5 的线性方程组. 解 按紧凑格式,并对 b 作统一处理,有(1)1 (2)2 (3)3 (4)4 (2)2 河南农业大学《数值计算方法》教案 汪松玉 - 64 - (1)1 (4)2 (9)3 (16)6 (10)4 (1)1 (8)6 (27)6 (64)4 (44)3 (1)1 (16)14 (81)36 (256)24 (190)1 即1123421201364,,166143114 36 24 0 1 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L U y 再求解 ? Ux y 求得原方程组的解为 4 3 2 1 1 , 1, 1, 1. x x x x ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 65 - 5.3 解三对角线性方程组的追赶法 利用三角分解,易于导出一些特殊方程组的解法. 但在 一些实际问题中,例如常微分方程数值解、求热传导方程及 三次样条插值函数等问题中都会遇到系数矩阵是三对角矩 阵的方程组,并且系数矩阵大都有对角占优性质,因此不必 选主元,就能保证算法的稳定性. 设有方程组 ? Ax d ,其中 A 为三对角矩阵,即1112222333111,nnnnnnbcdabcdabcabcdab?????????????????????????????????????????Ad(5.3.1) 河南农业大学《数值计算方法》教案 汪松玉 - 66 - 关于三对角矩阵 A ,有如下性质: 定理 5.2 当三对角矩阵(5.3.1)满足条件 1 1 0 0 , 0,( 2,3, , 1) 0 k k k k k n n b c b a c a c k n b a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.3.2) 时,其行列式det ( ) 0 ? A ,即A是非奇异的. 由于方程组中零系数占多数,因此在进行顺序消元过程 时有必要修改并简化消元计算的公式,免去大量不必要的计 算.作A的Crout 分解,并把d 组作为 A 的最后一列统一处 理,可得 河南农业大学《数值计算方法》教案 汪松玉 - 67 - 1 2 2 3 , n n l a l a a l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L 1 2 1 1 1 , 1 n u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? U 1 2 n y y y ? ? ? ? ? ? ? ? ? ? ? ? ? ? y 其中 1 1 0 0 , 1,2, , ) / , ( 1,2, , 1) 0. i i i i i i i i i i i i l b a u y d y a l i n u c l i n u y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.3.3) 解方程 ? Ux y ,得: 河南农业大学《数值计算方法》教案 汪松玉 - 68 - 1, 1, 2, 2,1. n n i i i i x y x y u x i n n ? ? ? ? ? ? ? ? ? ? ? (5.3.4) 上述计算过程也称为求解三对角线性方程组的追赶法, 它由两组递推公式组成:其消元过程是按 i 从小到大的次序 来计算 , , i i i l y u ,(5.3.3)称为追的过程;回代过程求 i x 的次序 正相反,(5.3.4)称为赶的过程. 容易算出追赶法的运算量仅 为8n . 解三对角方程组(5.3.1)的追赶法,它分为"追"和"赶" 两个过程: (1) 追的过程(消元过程):按(5.3.3)的顺序计算系数: 1 2 n l l l ? ? ? ? 、 1 2 n y y y ? ? ? ? 及121nuuu?????河南农业大学《数值计算方法》教案 汪松玉 - 69 - (2) 赶的过程(回代过程):按5.3.4)的逆序求出解: 1 1 n n x x x ? ? ? ? ? . 注 追赶法具有计算量少,方法简单,算法稳定等优点, 在计算机上实现时,只需用三个一维数组分别存放方程组 (5.3.1)的系数矩阵 A 的三条对角线上的元素{ i i i a b c , 两组工作单元保存{ }, { } i i u y . 例5.7 用追赶法求解三对角方程组 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 10 4 2 4 1 1 4 1 1 4 3 2 1 x x x . 解23123,1, n a a c c ? ? ? ? ? ? 10 , 4 , 2 , 4 3 2 1 3 2 1 ? ? ? ? ? ? d d d b b b . 河南农业大学《数值计算方法》教案 汪松玉 - 70 - (追) 4 15 1 , 4 56 0 1 15 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L 1 1 - 0 4 4 1 - , 15 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? U 1 2 6 5 3 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? y . (赶) 3 3 2 2 2 3 3, 6 5 4 15 3 2 x y x y u x ? ? ? ? ? ? ? ? , 1 1 1 2 1 2 1 4 2 1 x y u x ? ? ? ? ? ? 河南农业大学《数值计算方法》教案 汪松玉 - 71 - 5.4 对称正定矩阵的平方根法和 LDLT 分解 在很多实际问题中,归结为求解方程组 b Ax ? ,其中 n n R A ? ? 恰好局有特殊性质,为对称正定阵,即A.计算机 上求解这种类型方程组的有效方法之一就是平方根法,也就 是利用对称定矩阵三角分解的一种分解法. 定义 设A为n阶实方阵,如果 A 满足条件 (1) A AT ? ; (2) 对任意非零向量 n R x ? ,则有 0 T x Ax ? ; 则称 A 为对称正定矩阵. 关于对称正定矩阵有以下性质. 性质 1 如果 , n n R A ? ? 为对称正定阵,则(1) A 非奇异矩阵,且1?A亦是对称正定阵; 河南农业大学《数值计算方法》教案 汪松玉 - 72 - (2) 记kA为A的顺序主子阵,kA亦是对称正定阵),,2,1(nk??;(3) A 的特征值 ) , , 2 , 1 ( 0 ) ( n i A i ? ? ? ? ; (4) A 的顺序主子式都大于零,即). , , 2 , 1 ( 0 ) det( n k Ak ? ? ? 性质 2 (1)设nnRA??为对称矩阵,且A特征值 ( ) 0 i A ? ? ( 1, , i ? ? ) n ,则A为对称正定矩阵. (2) 设nnRA??为对称矩阵且 det( k A ) ) , , 2 , 1 ( 0 n k ? ? ? 则A为对称正定矩阵. 当A为n阶对称正定矩阵时,则A可分解成一个下三角 阵及其转置矩阵的乘积,即T?ALL (5.4.1) 河南农业大学《数值计算方法》教案 汪松玉 - 73 - 其中 L 是下三角阵. 事实上,将A作LU 分解, 11 12 1 21 22 2 1 2 1 1 1 n n n n nn u u u l u u l l u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A LU , 并写出其分块形式,有(1,2, , ) k k k k n ? ? ? A L U ,其中 11 12 1 21 22 2 1 2 , k k k k k kk a a a a a a a a a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A 21 1 2 1 1 , 1 k k k l l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? L 河南农业大学《数值计算方法》教案 汪松玉 - 74 - 11 12 1 22 2 k k k kk u u u u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? U 由于 A 的各阶主子行列式都大于零,所以有 1 det det det 0, ( 1,2, , ). k k k k ii i u k n ? ? ? ? ? ? ? ? A L U (5.4.2) 故有 0, ( 1,2, , ) ii u i n ? ? ? . 进一步将U 分解为 DR ,其中 河南农业大学《数值计算方法》教案 汪松玉 - 75 - 11 22 0 , 0 nn u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? D 1 12 11 11 2 22 1, 1, 1 1 1 1 n n n n n n u u u u u u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? R . 可以证明上述分解是唯一的.再由 T ? A A,可得 T T R DL = LDR , 故有 , T T ? ? R L L R,即T?ALDL (5.4.3) 其中 L 为单位下三角矩阵,D 为对角矩阵.称(5.4.3)为对称正 河南农业大学《数值计算方法》教案 汪松玉 - 76 - 定矩阵 A 的TLDL 分解. 若A是对称正定矩阵, 0, ( 1,2, , ) ii u i n ? ? ? ,将D进一 步分解为: ? D DD ,其中 11 22 nn u u u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? D 则(5.4.3)又可以表示成 ( )( ) T T T ? ? ? A LDDL LD LD LL 记?LLD ,得分解式 T ? A LL (5.4.4) 河南农业大学《数值计算方法》教案 汪松玉 - 77 - 这里 L 是下三角矩阵,并且对角元全是正数. 通常称(5.4.4)为对称正定矩阵 A 的乔累斯基(Cholesky) 分解,或LLT 分解. 当方程组(5.0.1)的系数矩阵 A 为对称正定矩阵时,可以 用乔累斯基分解法求解,即利用 A 的TLL 分解或 T LDL 分解 式来求解,分别又称为平方根法和 T LDL 分解法.用平方根法 解方程组的计算过程分两个步骤: (1) 计算分解式 T ? A LL .设11 11 21 1 21 22 22 2 1 2 n n T n n nn nn l l l l l l l l l l l l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A LL 河南农业大学《数值计算方法》教案 汪松玉 - 78 - 根据矩阵乘法有 1 2 2 2 1 1 ,( 1,2, , ) k k k k ks ks kk s s a l l l k n ? ? ? ? ? ? ? ? ? ? 以及 1 1 1 , k k jk js ks j s k s jk k k s s a l l l l l l ? ? ? ? ? ? ? ? ( 1, 2, , ) j i i n ? ? ? ? 于是有 1 2 1 1 1 , ( ) , ( 1, 2, , ) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? k kk kk ks s k jk js ks s jk kk l a l a l l l j k k n l (5.4.5) 河南农业大学《数值计算方法》教案 汪松玉 - 79 - 为了便于记忆,我们仿照杜利特尔分解的记忆方法,把元素 设想成如下排列 11 11 ( ) a l 12 12 ( ) a l 13 13 ( ) a l 14 14 ( ) a l ? 21 21 ( ) a l 22 22 ( ) a l 23 23 ( ) a l 24 24 ( ) a l ? 31 31 ( ) a l 32 32 ( ) a l 33 33 ( ) a l 34 34 ( ) a l ? 41 41 ( ) a l 42 42 ( ) a l 43 43 ( ) a l 44 44 ( ) a l ? ? ? ? ? ① ② ③ ④ 其中 , ( , 1,2, , ) ij ji ij ji a a l l i j n ? ? ? ? ,由(5.4.5)式可知, L 的 对角元 ii l 等于 A 的对应元素 ii a 减去一个内积,然后开平方. 内积的每一项仍然是同行元素 ik l 与同列元素 ki l 的乘积,但需 注意的是 ik ki l l ? . L 的非对角元素 ( ) ji l j i ? 的计算完全与 LU 河南农业大学《数值计算方法》教案 汪松玉 - 80 - 分解中 L 的元的计算相同,只是这里的对角元的名称是 ii l . 运算次序是按列进行. (2) 依次求解两个三角形方程组 ? Ly b及T?Lxy,即 ① b Ly ? 求y;②yxLT ? 求x;③11()/ ,( 1,2, , ) i i i ik k ii k y b l y l i n ? ? ? ? ? ? ? ; ④ 1 1, ,2,1) n i i ki k ii k i x y l x l i n n ? ? ? ? ? ? ? ? Cholesky 分解算法 设有 n n R A ? ? 为对称正定阵, 下述算 法计算下三角阵 L 使TLL A ? 且计算 ) ( j i lij ? 覆盖 ij a . 对于 n i , , 2 , 1 ? ? (行格式) 河南农业大学《数值计算方法》教案 汪松玉 - 81 - (1) 对于 1 , , 1 ? ? i j ? , jj j k jk ik ij ij ij a a a a l a / ) ( 1 1 ? ? ? ? ? ? (2) 2 / 1 1 1 2 ) ( ? ? ? ? ? ? i k ik ii ii ii a a l a 本算法大约需要 6 / 2 n 次乘除法. 由对称正定阵 A 分解计算公式(2)可得 ? ? ? ? i k ik ii n i l a 1 2 ) , , 2 , 1 ( , ? 于是, ii n i ii ik a a l ? ? ? ? 1 2 max ) , , 2 , 1 ; , , 2 , 1 ( i k n i ? ? ? ? 上式说明在不选主元的平方根法中,矩阵 L 元素有界, 或者说在分解过程中产生 L 元素 ik l 的数量级不会增长,且 河南农业大学《数值计算方法》教案 汪松玉 - 82 - ) , , 1 ( 0 n i lii ? ? ? ,由此可以说明,平方根算法是一个数值稳 定的方法. 对于正定矩阵的三角分解 T LL A ? 计算最约 6 / 3 n 次乘除 法,大约为 LU 分解计算的一半. 此外,计算过程式稳定的. 事实上,由关系式 2 1 ? ? ? k k k ks s a l 得||? ij ii l a 说明 Cholesky 分解中的量 ij l 能够得到控制, 因此其计算过程 是稳定的. 例5.8 用平方根法求解方程组 河南农业大学《数值计算方法》教案 汪松玉 - 83 - ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 7 8 5 11 0 2 0 2 1 2 1 1 3 2 1 x x x 解设(1)1 (1)1 (2)2 (1)1 (2)1 (0)-2 (2)2 (0)-2 (11) 3 于是,有100110223????????????L.由 ? Ly b 解得 1 2 3 5, 3, 3 y y y ? ? ? , 河南农业大学《数值计算方法》教案 汪松玉 - 84 - 由T?Lxy解得 1 , 5 , 2 3 2 1 ? ? ? ? x x x . 为了避免开方运算,有时直接使用对称矩阵 A 的TLDL 分解来计算,令T?ALDL 1 21 1 21 2 2 1 2 1 1 1 1 1 1 n n n n n d l l l d l l l d ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 根据矩阵乘法可以求出 L 和D的元素,然后将方程组(5.0.1) 即T?LDL x b 转化为两个三角形方程组 ? Ly b, 1 T ? ? L x D y , 由前一方程解出 y ,代入后一方程便可解出 x .具体计算步骤 可参看有关书籍.