第五章 求解线性方程组的迭代法

本章研究的对象是n阶线性代数方程组 在自然科学和工程技术中 许多问题的解决常常归结为线性代数方程组的求解 例如 电学中的网络问题 船体放样中三次样条求解问题 求解非线性方程组问题 常微分方程及偏微分方程数值求解问题等 都可归结为求解线性代数方程组 有关线性代数方程组解的存在唯一性及解的结构等理论问题 线性代数已作了详细讨论 本书只介绍线性代数方程组的两类求解方法 直接法和迭代法 直接法是指经过有限步计算后求得方程组精确解的方法 第五章求解线性方程组的直接法 1 1 若用矩阵和向量的记号来表示 若A是非奇异矩阵 则方程组有唯一解 记D detA 应用Cramer法则可得 Di是用b代替中D第i列而得到的相应行列式 在实际中 用Cramer法则求解 1 1 并不可行 P98 因此 用Cramer法则求解 1 1 的乘除次数为 n阶行列式共有n 项 每项n个因子 所以计算一个n阶行列式共需 n 1 n 次乘法 应用Cramer法则求解方程组需要计算 n 1 个行列式 另外还需n次除法 当n 10时 N 359251210 当n 20时 N 9 7073 1020 P98 2 1高斯消去法 2 1 1顺序高斯消去法 高斯消去法是一个古老的方法 但实践证明它仍是目前计算机上一种常用有效的方法 其基本思想是通过初等变换将方程组 1 1 转化为一个等价的三角形方程组 1 3 这个过程称为消元 得到 1 3 后 就可逐个求出这个过程称为回代 从 1 3 中的最后一个方程直接可求出 代入倒数第二个方程便有 一般有 具有过程如下 2 1 1 1顺序高斯消去法的计算过程 P100 为了统一符号 把方程组 1 1 改写成下面形式 用矩阵表示即为 其中 则 1 4 化为 将 1 4 第1行乘后加到第i行 i 2 3 n 其中 用矩阵符号表示为 P101 其中 将 1 6 中第2个方程乘后加到第i行 i 3 4 n 类似地 若 则 1 6 化为 其中 用矩阵符号表示为 其中 P102 若上述过程可以继续进行下去 重复n 1次后 可得到等价的三角形方程组 用矩阵符号表示为 其中 P102 消元过程中元素的计算公式可以归纳为 从上面消元过程我们可以看到 要顺序高斯消去法能进行下去的充分必要条件是A矩阵的所有顺序主子式不为零 P102 回代过程算法 P102 P102 2 1 2顺序高斯消去法求解算法 上述算法中 为了不增加存储量 消元过程直接在系数矩阵A中进行 的引用主要是为了控制溢出 2 1 3高斯消去法的计算量 1 消元过程的计算量 第1次消元 计算mi1 i 2 3 n 需要n 1次除法 计算 i j 2 3 n 需要次乘法 计算 i 2 3 n 需要次乘法 因此第一次消元的乘除次数为 加减法次数为 P103 第2次消元等价于一个阶方程组的一次消元 第2次消元后乘除次数为 加减法次数为 因此消元过程的计算量为 乘除法次数 加减法次数 P103 2 回代过程的计算量 从而顺序高斯消去法的计算量为 乘除法次数 加减法次数 P103 在上述消元过程中 用到了矩阵初等变换 即将某一行减去另一行的常数倍 在后面的列主元和全主元高斯消去法中我们还会用到矩阵的行交换与列交换 今对相应初等变换的矩阵表示说明如下 2 1 4顺序高斯消去法的矩阵解释 矩阵表示 利用单位坐标向量ei i 1 2 n 交换单位矩阵I的i j两行所得矩阵记为Pij 单位矩阵I可以表示成 1 交换矩阵的两行 两列的变换 P104 用Pij左乘矩阵A 其作用就是交换A的第i行和第j行 用Pij右乘矩阵A 其作用就是交换A的第i列和第j列 2 矩阵的第i i 2 3 n 行减去mi1乘第1行的变换 这种变换相当于用矩阵 左乘矩阵A 高斯消去法的消元过程 从矩阵运算的角度来看 从 1 5 化为 1 7 相当于系数矩阵和左端项左乘一个初等变换矩阵M1 即 类似地 记 P105 则 P105 从而 注意到 均为下三角矩阵且对角线元素为1 因此 为对角线元素为1的下三角矩阵 记 则 即系数矩阵A的所有顺序主子式不为零 其中L为对角线元素为1的单位下三角矩阵 只要 由顺序高斯消去法我们看到 U为上三角矩阵 P111 定理1 1 如果n阶矩阵A的所有顺序主子式不为零 则A有的LU分解 第一部分 则A一定可以分解成 1 1 2列主元高斯消去法 主元高斯消去法根据主元选取范围的不同 又分为列主元消去法和全主元消去法 P105 在顺序高斯消去法的消元过程中 若出现 则消元无法进行 即使 就会导致其它元素量级的巨大增长和舍入误差的扩散 引起解的失真 用它作为除数 先看一个例子 主元高斯消去法是为了控制舍入误差而提出来的一种算法 例1 方程组 该方程组精确解为 解利用顺序高斯消去法取4位有效数进行运算得 将第1个方程的 1764倍加到第2个方程 我们得 P105 通过回代我们得 x2 1 001 x1 10 00 与准确解相差很大 这是由于用0 003000作除数 将第1个方程的 0 0005760倍加到第2个方程 我们得 我们交换一下方程的位置 再利用顺序高斯消去法求解 回代 我们得 x1 10 00 x2 1 000 因此 为使这种不稳定现象发生的可能性减至最小 必须在每次消元前选择一个相对较大的系数作为主元 1 1 2 1列主元高斯消去法计算过程 第一步 在系数矩阵的第1列元素中选出绝对值最大的元素称之为第1列的主元 即 P106 小主元是不稳定的根源 如果i1 1 则交换第1行和第i1行元素及相应的右端项 交换完后再进行消元得到 1 12 第二步 在系数矩阵第2列以下n 1个元素中选主元 即 如果i2 2 则对调 1 12 中的第2行与第i2行元素及相应的右端项 交换完后再进行消元得 一般地 第k步 求Ik使 如果ik k 则交换第k行与第Ik行及相应的右端项 经消元后得 列主元消去法除了每步需要按列选主元并可能进行矩阵行交换外 其消元过程与顺序高斯消去法的过程一样 归纳以上过程有以下算法 P106 2 2 2列主元高斯消去算法 1 输入A b及 P106 3 回代 4 打印 i 1 2 n 2 2 3列主元高斯消去法的计算量 列主元高斯消去法乘除及加减次数与顺序高斯消去法相同 增加了选主元时n n 1 次比较以及交换方程次序所需时间 1 1 2 4列主元高斯消去法的矩阵解释 记 则说明第1个主元在第i1行 交换第1行与第i1行 P107 相当于用左乘矩阵 然后再进行消元 用矩阵表示即为 同样道理 我们有 由此得到等价的三角形方程组为 这时 例3解方程组 它的准确解为 解第1步 首先选取第1列绝对值最大的主元 对于本方程组 绝对值最大的主元为 18 交换第1行和第2行得 P107 将第1个方程的2 3倍加到第2个方程 第1个方程的1 18倍加到第3个方程得 P108 第2步 选取第2列 第2个元素以下 绝对值最大的主元 对于本方程组 绝对值最大的主元为7 6 交换第2行和第3行消元后得 把上面方程进行回代 就可逐步解出 2 3全主元高斯消去法 但找主元和交换行列次序要花费大量的机器时间 因此有时我们不用全主元高斯消去法而用列主元高斯消去法 P108 全主元高斯消去法与列主元高斯消去法类似 所不同的是选取主元的范围不同 全主元高斯消去法在整个系数矩阵中找绝对值最大的元素作为主元 这样可控制舍入误差的增长 将舍入误差控制在一个最小的范围 2 3 1全主元高斯消去法计算过程 若j1 1 则交换第1列和第j1列 左乘 第一步 求i1 j1使得 并交换x1与 交换完后再进行消元得到 1 13 如果i1 1 则交换第1行和第i1行及相应的右端项 左乘 P108 记 上式变为 一般地 第k步 k 1 2 n 1 求ik jk使得 如果ik k 则交换第k行和第ik行及相应的右端项 若jk k 则交换第xk列和第jk列 交换完后再进行消元得到 记 上式变为 归纳上述过程有以下算法 2 3 2全主元高斯消去法算法 P109 2 3 3全主元高斯消去法的计算量 全主元高斯消去法乘除及加减次数与顺序高斯消去法相同 增加了选主元时 比较以及交换方程次序所需的时间 P110 2 3 4全主元高斯消去法的矩阵解释 记 由上述过程我们知道 主元在第i1行 第j1列 交换第1行与第i1行 第1列与第j1列 然后再消元即可得到A 2 用矩阵表示即为 P110 同理我们得 两边同乘 注意到 由于 由 1 14 知 1 15 通过 1 17 可求得Z 由 1 16 最终求得 2 2LU分解法 而从消元过程我们可以看到从A 1 转化到A n 的过程实际是进行若干次初等变换的过程 由顺序高斯消去法的矩阵解释我们得到 从而 注意到 P110 前面介绍过的顺序高斯消去法可以把方程组 转化为等价的一个上三角形方程组 均为下三角矩阵且对角线元素为1 因此 为对角线元素为1的下三角矩阵 由顺序高斯消去法我们看到 只要方程 的系数矩阵A的所有顺序主子式不为零 其中L为对角线元素为1的单位下三角矩阵 则A一定可以分解成 记 U为上三角矩阵 如果n阶矩阵A的所有顺序主子式不为零 设 P111 定理1 1 证明 则A有唯一的LU分解 第二部分 其中L L1为单位下三角矩阵 U U1为上三角矩阵 因为A非奇异 于是有 故存在 上式左端为单位下三角矩阵 右端为上三角矩阵 从而上式两边均为单位矩阵 从而证得 如果能实现这种分解 则求解方程组 就相当简单 注意到A LU 因此 则 用到上 下 三角阵的乘积 逆阵仍为上 下 三角矩阵 U U1 L L1 A的三角形矩阵分解的唯一性获证 记 i 2 n 解 我们得 我们得 解 2 2 1直接LU分解法 下面将详细讨论矩阵A的LU分解 设A 即 要两个矩阵相等 其对应元素应相等 利用矩阵乘法 由矩阵的第 行对应相等得 P111 由矩阵的第 列对应相等得 这就求出了U矩阵的第1行和L矩阵的第1列元素 一般地 设U矩阵的前k 1行和L矩阵的前k 1列已经求出 则由 我们得 又由 我们得 综上所述 A的LU分解公式如下 k 1 2 n 上述计算公式有如下特点 U的元素按逐行求 L的元素按逐列求 先求U的第k行元素 然后求L的第k列元素 U和L的元素一行一列交叉进行 例4P114 练习P126EX2 2 1 1直接LU分解法算法 P113 每次将计算结果和仍存放在矩阵A的相应元素和所占的单元内 不必占用新的单元 只要将上述算法中出现的和相应变换为和即可 2 1 2直接LU分解法的计算量 乘除次数 LU分解计算U的计算量 计算L的计算量 总的计算量为 求解的计算量 求解的计算量 P114 例4用LU分解求解方程组 解先把系数矩阵进行LU分解得 求解 得y1 0 y2 3 y3 0 5 求解 解得x1 1 x2 1 x3 1 故AX b的解为 P114 用直接LU分解法求解方程组所需要的计算量仍为 和高斯消去法乘除法次数相同 但LU分解把对系数矩阵的计算和对右端项的计算分开了 这就使得求解系数矩阵相同而右端项不同的一系列方程组变得特别方便 2 2 2列主元LU分解法 假若第k 1步的分解已完成 在进行第k步分解时 为避免出现小的作除数 计算 P115 选取行号 使 若 则对调A的第k行与第行 再求 于是完成了U的第k行和L的第k列的计算 算法如下 P115 列主元LU分解算法 2 3对称正定矩阵的平方根法和LDLT分解法 当A是对称正定矩阵时 存在一个实的非奇异下三角矩阵L1使 且当限定L1的对角线为正时 这种分解是唯一的 这种分解称为矩阵的Cholesky分解 定理1 2设A是对称正定矩阵 则A有如下分解 其中 L是单位下三角阵 D为对角阵 且这种分解是唯一的 P116 证明因为A对称正定 从而A有唯一的LU分解 因为 故上式可化为 将右端三角矩阵分别记为L D和R 因A对称有 A LDR RTDLT 因为A的LU分解唯一 故L RT LT R A LDLT 此分解显然是唯一的 定理1 3n阶对称正定矩阵A一定有Cholesky分解 当限定L1的对角线为正时 矩阵的Cholesky分解唯一 P116 证明 使 又记 它是一个非奇异下三角阵 且对角线元素为正 故A有Cholesky分解 分解的唯一性显而易见 下边推导Cholesky分解的计算公式 设A LLT 即 其中 i j 1 2 n 第1步 由矩阵乘法有 且 故求得 i 2 3 n 一般地 设L矩阵的前k 1列元素已求出 则 i 1 2 n 第k步 由矩阵乘法得 P117 由于分解公式 1 18 中的每一步都有开方运算 故也称Cholesky方法为平方根法 注意到 所以 这说明在分解过程中元素的平方不会超过A的最大对角元 因而舍入误差的放大受到限制 所以平方根法求解对称正定方程组时可以不考虑选主元的问题 可以证明 若用顺序高斯消去法求解对称正定方程组AX b 则有 k 1 2 n 其中是第k步顺序高斯消去法过程所得到的元素 这说明高斯消去法求解对称正定方程组也可以不用选主元 从运算的角度看 平方根法是有利的 用平方根法求解AX b所需的乘除次数为 另外还有n次开平方运算 乘除次数只是高斯消去法的一半左右 例5用平方根法解方程组 解 由平方根分解法可得 P117 求解 得 求解 得 平方根求解算法 P118 为避免平方根法的开方运算 利用定理1 2 也可以对A作LDLT分解 其分解方法可借助于前面所说的LU分解法 从定理1 1我们可以看到 当A的所有顺序主子式不为零且对称时 A的LU分解法存在唯一 并且 U DLT 其中 D diag u11 u22 unn 利用LU分解法先计算出U的第k行 由U的第k行我们得到D的第k个对角元素及L的第k列 因此LDLT分解的乘除法计算次数为 大概只是高斯消去法的一半 与平方根法大致相同 但避免了平开运算 例6用分解求方程组 解由分解法得 P119 解 得 解 得 最后解 得 这样就得到 LDLT分解算法 P120 P126EX3 P126EX3 2 5向量与矩阵范数 用直接法求解线性方程组AX b时 由于有舍入误差 只能得到近似解 为了进行解的误差分析和以后要介绍的迭代法的收敛性分析 我们首先介绍向量范数和矩阵范数的概念 P120 定义1 1 设X Rn 表示定义在Rn上的一个实值函数 称之为X的范数 它具有下列性质 3 三角不等式 即对任意两个向量X Y Rn 恒有 1 非负性 即对一切X Rn X 0 0 2 齐次性 即对任何实数a R X Rn 向量范数是n维欧几里德空间中长度概念的推广 P121 2 4 1向量的范数 设X x1 x2 xn T 则有 4 更一般的p范数为 1 向量的1范数 2 向量的2范数 3 向量的 范数 容易验证 上面几种范数的定义满足范数的三个条件 常用的范数 P121 定义1 2 在Rn上定义的任一向量范数都与范数等价 即存在正数M与m M m 对一切X Rn 不等式 成立 定理1 4Rn上定义的任何两个范数都是等价的 证 对常用范数 容易验证下列不等式 设给定Rn中的向量序列 即 其中 若对任何i i 1 2 n 都有 则向量 称为向量序列 的极限 或者说向量序列 依坐标收敛于向量 记为 定义2 3 定理1 5向量序列 Xk 依坐标收敛于X 的充要条件是 向量序列依范数收敛与依坐标收敛是等价的 P121 2 4 2矩阵的范数 1 当A 0时 0 当A 0时 0 2 对任意实数和任意A 有 3 对任意两个n阶矩阵A B有 对任意两个n阶矩阵A B 有 定义 则称 A 为矩阵A的范数 P122 设 若 的非负实数满足下列四个条件 设是中矩阵序列 A是中矩阵 则 其中和分别表示和A的分量 定理2 6 定义1 4对于给定的向量范数 与矩阵范数 如果 A Rn n和 X Rn 满足 则称此矩阵范数与向量范数是相容的 P122 并称其为由向量 诱导出的矩阵范数 现在考察由向量范数 1 2 产生的矩阵范数 P122 定理2 7 设A为n阶方阵 Rn中已定义了向量范数 则称为矩阵A的范数或模 记为 定理2 8设n阶方阵A aij n n 则 与相容的矩阵范数是 与相容的矩阵范数是 其中 1为矩阵ATA的最大特征值 与相容的矩阵范数是 证明略 这3种范数中 1范数与范数的计算比较简单 2范数要计算矩阵的特征值 比较复杂 但2范数有一些好的性质 常用于理论分析 P122 P123 此外 在上的一个常用且易于计算的矩阵范数为 通常称作Frobenius范数 它是上的2范数的自然推广 A的范数与A的特征值之间的关系 矩阵A的任一特征值绝对值不超过A的范数 同向量范数一样 矩阵范数也有下面定理 定理1 6 设是中矩阵序列 A是中矩阵 则 其中和分别表示和A的分量 2 5 3谱半径 P123 定义1 5 矩阵A的绝对值最大的特征值称为A的谱半径 记为 定理1 9 定理1 10设A是任意阶矩阵 由A的各次幂所组成的矩阵序列 收敛于零 即的充分必要条件是 P124 证明从略 读者可参阅冯康等编的 数值计算方法 求解时 A和b的误差对解x有何影响 设A精确b有误差 得到的解为 即 绝对误差放大因子 相对误差放大因子 线性方程组的性态和解的误差分析 1 4 4方程右端误差对解的影响 P124 P124 在这种情形 数就是误差的放大率 即解的相对误差不超过