线性方程组直接解法.ppt

第三章 线性方程组直接解法 第三章目录 1 Gauus消元法 2 主元素法2 1引入主元素法的必要性2 2列主元素法2 3全主元素法2 4解三对角方程组的追赶法 3 矩阵分解法3 1Gauss消去法的矩阵形式3 2矩阵的三角分解3 3直接三角分解法 4 平方根法与改进的平方根法 5 矩阵求逆 6 方程组的性态和条件数 设n阶线性方程组 其矩阵形式为 Ax b 2 2 其中 在科学研究和工程技术中所提出的计算问题中 线性方程组的求解问题是基本的 常见的 很多问题如插值函数 最小二乘数据拟合 构造求解微分方程的差分格式等 都包含了解线性方程组问题 因此 线性方程组的解法在数值计算中占有较重要的地位 求解Ax b 曾经学过高斯 Gauss 消元法 克莱姆 Cramer 法则 矩阵变换法等 但已远远满足不了实际运算的需要 主要体现两个方面 一是运算的快速和准确 其次是方程组的个数增大时的计算问题 如何建立能在计算机上可以实现的有效而实用的解法 具有极其重要的意义 我们也曾指出过 Cramer法则在理论上是绝对正确的 但当n较大时 在实际计算中却不能用 如果线性方程组Ax b的系数行列式不为零 即det A 0 则该方程组有唯一解 线性方程组的数值解法 解线性方程组的数值方法大致分为两类 请注意 由于在计算中某些数据实际上只能用有限位小数 即不可避免地存在着舍入误差的影响 因而即使是准确解法 也只能求到近似解 直接法在求解中小型线性方程组 100个 特别是系数矩阵为稠密型时 是常用的 非常好的方法 直接法 指假设计算过程中不产生含入误差 经过有限步四则运算可求得方程组准确解的方法 2 迭代法 从给定的方程组的一个近似值出发 构造某种算法逐步将其准确化 一般不能在有限步内得到准确解 这一章介绍计算机上常用的直接法 它们都是以Gauss消元法为基本方法 即先将线性方程组化为等价的三角形方程组 然后求解 1Gauss消元法 Gauss消元法是最基本的一种方法 下例说明其基本思想 例1 解线性方程组 解 消去x1 进行第一次消元 首先找乘数 以 12乘第一个方程加到第二个方程 以18乘第一个方程加到第三个方程上可得同解方程组 例1 续 上述Gauss消元法的基本思想是 先逐次消去变量 将方程组化成同解的上三角形方程组 此过程称为消元过程 然后按方程相反顺序求解上三角形方程组 得到原方程组的解 此过程称为回代过程 再消一次元得 二次消元后将方程化为倒三角形式 然后进行回代容易解出 x3 3 x2 2 x1 1 我们的目的 是要总结归纳出一般情况下的n阶线性方程组的消元公式和回代求解公式 从而得到求解n阶线性方程组的能顺利在计算机上实现的行之有效的算法 为能更清楚地得到算法 下面以4阶线性方程组为例总结求解步骤 并且很容易地可推广至一般的n阶线性方程组 可以检查 分别以 li1乘第一个方程加到第i个方程上可以完成第一次消元 得同解方程组 变化以后的方程组系数及右边的常数项可总结出如下的计算公式 Gauss消元法的基本步骤3 4阶 以方程组中第i个方程减去第二个方程乘li2 i 3 4 完成第二次消元 上标为3的系数和右端项可由下面公式计算 第三步 消元 4阶方程组需进行3次消元 将上述A 3 X b 3 中最后一个方程中的x3消为零 然后可回代求解 由于A 4 为上三角形 所以可按变量的逆序逐步回代求原方程组的解 上述消元 回代求解过程很容易推广到一般的n阶线性方程组 经过上述消元步骤 得到同解的上三角形方程组 A 4 x b 4 Gauss消元法的消元过程1 2 n阶 一般地 设n阶方程组 消元过程为 第k步消元后同解方程组中上标为k 1的元素的计算公式见下屏 照此消元下去 完成n 1次消元后 可将原方程组化成同解的上三角形方程组如下 Gauss消元法的回代过程 n阶 回代过程 逐步回代求得原方程组的解 Gauss消元法的计算量 由于在计算机中作乘除运算量所需时间远大于作加减运算所需时间 故只考虑作乘除运算量 由消元法步骤知 第k次消元需作n k次除法 作 n k n k 1 次乘法 故消元过程中乘除法运算量为 所以Gauss消去法的乘除法总运算量为 Gauss法与Cramer法则的计算量比较 Gauss消元法的乘除法总运算量为 与我们曾经介绍的Cramer法则的乘除法总运算量 n2 1 n n相比 由下表可知 当阶数越高时 Gauss消元法所需乘除法次数比Cramer法则要少得多 Gauss消元法的优缺点 但其计算过程中 要求akk k 称为主元素 均不为零 因而适用范围小 只适用于从1到n 1阶顺序主子式均不为零的矩阵A 计算实践还表明 Gauss消元法的数值稳定性差 当出现小主元素时 会严重影响计算结果的精度 甚至导出错误的结果 Gauss消元法简单易行 2主元素法 2 1引入主元素的必要性对线性方程组AX b 若其系数行列式det A 0 则该方程组有唯一解 但是这一条件不能保证所有主元素都不等于零 只要某一主元素等于零 就不能用Gauss消元法求解该方程组 即使所有主元素不等于零 但某一主元素的绝对值很小时 Gauss消元法也是不适用的 如下例 例2 例2 续1 解 为减小误差 计算过程中保留3位有效数字 按Gauss消元法步骤 第一次消元后得同解方程组 第二次消元后得同解方程组 回代得解 x3 2 02 x2 2 40 x1 5 80 容易验证 方程组 3 8 的准确解为 x1 2 60 x2 1 00 x3 2 0 显然两种结果相差很大 若在解方程组前 先交换方程的次序 如将 2 8 交换一行与三行改写成如下所示 再用Gauss消元法 顺序消元后得同解方程组 回代得解 x3 2 00 x2 1 00 x1 2 60 与准确解相同 例2 续2 例2两种解法的误差分析 在例2中 对 2 8 的方程进行顺序消元时 主元a 1 11 0 50 a 2 22 0 100都比较小 以它们作除数就增长了舍入误差 而导致计算结果不准确 产生上述现象的原因在于舍入误差 当 a k kk 很小时 进行第k次消元 要用 a k kk 作除数 这样就可能增大舍入误差造成溢出停机 或者导致错误的结果 为了在计算过程中 抑制舍入误差的增长 应尽量避免小主元的出现 如例2的第二种解法 通过交换方程次序 选取绝对值大的元素作主元基于这种思想而导出主元素法 2 2列主元素法 为简便起见 对方程组 2 1 用其增广矩阵 表示 并直接在增广矩阵上进行运算 列主元素法的具体步骤如下 列主元素法 如此经过n 1步 增广矩阵 2 9 被化成上三角形 然后由回代过程求解 在上述过程中 主元是按列选取的 故称为列主元法 例2中的第二种解法就是按列主元法进行的 2 3全主元素法 经过n k次消元后 得到与方程组 2 1 同解的上三角形方程组 再由回代过程求解 例6 计算过程保留三位小数 2 4解三对角方程组的追赶法 在很多问题中 需要解如下形式的三对角方程组 三对角方程组的系数矩阵为三对角阵 对于这种特殊而又简单的方程组 用前面介绍的方法求解由于有大量的零元素既占内存又浪费计算时间 显然很不经济 充分注意到三对角方程组的特点 根据顺序消元的思想导出一个简便的算法 追赶法 首先进行顺序消元 且每步将主元系数化为1 将方程组化为 其中系数按下式计算 然后回代求解 得 上述追赶法能进行到底 追赶法举例 用追赶法解下列三对角方程组 例4 解 首先将方程组化为 先追 然后回代 赶 求解 x5 0 x4 30 7 x3 6 7 x2 12 7 x1 0 可以看出 追赶法本质上还是顺序消元法 但由于计算过程中只涉及系数矩阵的非零元 因此大大节约了计算机内存与计算量 按乘除法次数进行比较 Gauss消元法约为n3 3 全主元法为n3 2 而追赶法仅为5n 3次 可见追赶法是求解三对角方程组的非常好的方法 3矩阵分解法 如果用矩阵形式表示 Gauss消元法的消元过程是对方程组 2 1 的增广矩阵 A b 进行一系列的初等行变换 将系数矩阵A化成上三角矩阵的过程 也等价于用一串初等变换阵去左乘增广矩阵 因此 消元过程可以通过矩阵运算来实现 紧接下屏 3 1Gauss消元法的矩阵形式 事实上 Gauss消元法的第一次消元相当于用初等矩阵 第二次消元相当于用初等矩阵 第k次消元相当于用初等矩阵 Gauss消元法的矩阵形式 经过n 1步消元后得到 因为Lk k 1 2 n 1 均为非奇异阵 故它们的逆矩阵存在 容易求出 这说明 在的条件下 消元过程实际上是把系数矩阵A分解成单位下三角阵与上三角矩阵的乘积的过程 Gauss消元法的矩阵形式 续 杜利特尔 Doolittle 分解 LU分解 事实上 只要A满足一定条件 由上述结论 它一定可以分解成两个三角形矩阵的乘积 即 A LU 上述分解称为杜利特尔 Doolittle 分解 也称为LU分解 当系数矩阵完成三角分解后 对于求解方程组 Ax b 消元过程相当于分解A LU及求解三角形方程组Ly b 回代过程则是求解另一个三角形方程组Ux y 因此 解线性方程组问题可转化为矩阵的三角分解问题 其中 L为单位下三角矩阵 U为上三角矩阵 3 2矩阵的三角分解 正如Gauss消元法要在一定条件下才能进行到底一样 矩阵A也必须满足一定条件才能进行三角分解 设A为n阶方阵 若A的顺序主子式Ai i 1 2 n 均不为零 则矩阵A存在唯一的Doolittle分解 定理2 1 下面讨论如何对A进行LU分解 1行1列 由于两个矩阵相等就是它们的对应元素都相等 因此通过比较A与LU的对应元素 即可得到直接计算L U的元素的公式 A LU即 紧接下屏 对A进行LU分解 由矩阵乘法规则及比较 2 11 两端的元素 得 即可由 2 12 求出U的第一行 L的第一列 下面讨论一般情况 即 U的第i行 L的第j列 一般情况下 可由 对A进行LU分解 r行r列 计算过程应按U第1行 L第1列 第1框 U第2行 L第2列 第2框 的顺序 具体分解步骤见下屏 得到计算uij和lij的公式 对A进行LU分解的具体步骤 1 计算U的第1行 L的第1列 亦称为计算第1框 2 计算U的第r行 L的第r列 r 2 n 即第r框 矩阵A的LU分解举例 解 按分解公式 2 13 一框一框分解 每框计算时先行后列 所以 例5 3 3直接三角分解法 若线性方程组Ax b的系数矩阵A完成三角分解 A LU 那么解方程组Ax b等价于求解两个三角形方程组Ly b Ux y 即由 再由 可解得 直接三角分解法 续 容易看出 式 2 14 与式 2 13 的运算规律相同 或者由 解线性方程组举例 例6 解 按表2 3计算 三角分解法的几点说明 1 用三角分解法求线性方程的乘除法运算量也是n3 3数量级 由于在求出uij lij和yi后 aij和bi无需保留 故上机计算时 可把L U和y存在A b所占的单元 回代时x取代y 整个计算过程中不需要增加新的存贮单元 3 完成A LU分解后可以较容易地求出行列式 A 的值 2 从三角分解法的推导及例中可以看出 系数矩阵的三角分解与右端项无关 因而在计算多个系数矩阵为A而右端不同的线性方程组系时 用三角分解法更为简便 如可用于求逆矩阵 三角分解法的几点说明 续 6 分解法的优点除上述2 3外 还有 a 可求解A2z b 因为算A2计算量大 可用 b 可根据A的形状设计算法 当A为大型稀疏 且非零元素有规律如带状 三对角等 作分解时能充分利用A的特点 L U能保持A的原状 即L与A的下三角相同 U与A的上三角形状相同 4 三角分解法一般也可采用选主元的技术 以使算法更具数值稳定性 5 也可以把矩阵A分解成一个下三角矩阵与一个单位上三角矩阵的乘积 矩阵的这种分解称为克劳特 Crout 分解 解线性方程组举例 例 下述矩阵能否分解为LU 其中L为单位下三角阵 U为上三角阵 若能分解 那么分解是否惟一 4平方根法与改进的平方根法 对称正定矩阵概念 工程实际问题的计算中 线性方程组的系数矩阵常常具有对称正定性 即其各阶顺序主子式及全部特征值均大于零 矩阵的这一特性使它的三角分解也有更简单的形式 从而导出一些特殊的解法 5 A的所有顺序主子式均为正数 即det Ak 0 k 1 2 n 4 A的所有特征值 0 3 A的对角线元素aii 0 i 1 2 n 2 A是非奇异阵 且A 1也是对称正定阵 1 A的所有顺序主子矩阵Ak k 1 2 n 也是对称正定阵 若A为对称正定矩阵 则 4 1平方根法 定理2 2 证明 首先A可直接作LU分解 记为A LU1 其中 定理2 2 续 其中U0是单位上三角阵 则由A的对称性可得 再由唯一性得U0 LT 所以A LDLT 并且这种分解是唯一的 又由于U1的对角线元素Ukk就是Gauss消元法的主元素 由于LD1 2是下三角阵 因此有推论 Choleskg分解1 推论 若A是对称正定矩阵 则存在唯一的主对角线元素都为正的下三角阵L 使得 A LLT 矩阵的这种分解称为Choleskg分解 用比较法可以导出L的计算公式 设 比较A与LLT的相应元素 即由A LLT可得 计算顺序按列进行 因此 Choleskg分解2 当完成矩阵A的Cholesky分解后 求解方程组AX b就转化求 求解对称正定线性方程组的方法称为平方根法 也称为Cholesky分解法 这种方法无需选主元 计算过程也是稳定