计算方法线性方程组解法.ppt

计算方法 力学系本科生 第三章线性方程组解法 SolutionforLinearAlgebraicEquations 3 1问题的提出 第三章线性方程组解法 n阶线性方程组 3 1问题的提出 3 1问题的提出 线性方程组Ax b 其中A是n维方阵 x是n维未知数向量 b是n维常数向量 3 1问题的提出 如果A是非奇异阵时 方程组有唯一解 且可以用克莱姆 Grammer 法则表示 其中xi是解向量x 的第i个分量 D detA Di是用b代替A的第i列后得到矩阵的行列式 3 1问题的提出 克莱姆方法求解计算量太大 需要计算 n 1 个n阶行列式 共需要 n 1 次乘法运算 3 1问题的提出 求解线性方程组的数值方法有两大类 直接法 directmethods 经过有限次算术运算可求方程组精确解的方法 实际上 由于舍入误差不可避免 一般得不到精确解 适合于求解低阶稠密阵方程组 3 1问题的提出 2 迭代法 iterativemethods 采用极限过程去逐步逼近线性方程组精确解的方法 迭代法需要计算机存储单元较少 对计算机要求不高 程序设计简单 但有收敛性和收敛速度方面的问题 迭代法是求解大型稀疏矩阵方程的重要方法 3 1问题的提出 我们在本章将要学习迭代法有 雅可比 Jacobi 迭代法 高斯 塞德尔 Gauss Seidel 迭代法 超松弛迭代法 Successiveoverrelaxationmethod SOR 3 1问题的提出 追赶法 Forwardeliminationandbackwardsubstitution 我们在本章将要学习直接解法有 高斯消去法 GaussElimination 高斯主元素消去法 GaussEleminationwithpivoting 三角分解法 LUdecomposition 3 1问题的提出 历史注记 线性代数方程组数值解法有着悠久的历史 我国古代数学著作 九章算术 公元1世纪 的 方程 章中就有了较好的线性方程组数值解法 相当于现代对方程组的增广矩阵进行初等变换 消去未知数的方法 中世纪的印度数学家也可以求解线性方程组 例如12世纪的婆什迦罗的著作中 也有求解线性方程组的内容 3 1问题的提出 在欧洲 16世纪的比特奥在其 算术 1559 中采用了与 九章算术 类似的消元法 日本数学家关孝和在其 解伏题之法 一书 1683 中首先采用了类似于现代行列式法求解了三元线性方程组 稍后 莱布尼茨提出关于行列式解线性方程组的思想 1693 1721年马可劳林用行列式展开式的方法给出了二元 三元 四元线性方程组的解法 3 1问题的提出 但他的符号记法不完善 1750年 克莱姆给出了现在比较通用的线性方程组行列式解法 即克莱姆法则 1764年 贝祖用行列式建立了线性方程组的一般理论 但由于当时计算的效率很低 这一理论几乎只有理论的意义 实际上只能求出未知数很少的线性代数方程组的解 只是在20世纪中叶电子计算机问世并投入应用之后 大型线性代数方程组的数值求解才成为可能 3 1问题的提出 如何利用计算机更精确 更有效地求解大型线性方程组 是计算数学中最重要的课题之一 现代计算实践中 常用的线性代数方程组数值解法有直接法和迭代法两大类 直接法是在没有舍入误差的假设下 经过有限次运算就可得出方程组的精确解的方法 如各种消元法 迭代法则采用逐次逼近的方法 即从一个初始值出发 按照一定的计算格式 迭代公式 构造一个向量的无穷序列 其极限才 3 1问题的提出 是方程组的精确解 用有限次运算得不到精确解 迭代法是牛顿最先提出来的 1940年经司威尔提出的松弛法也是一种迭代法 共轭梯度法则是另一种迭代法 是弗莱彻等人于20世纪60年代提出来的 3 1问题的提出 例3 1 3 1问题的提出 精确解为 将方程写为 取 3 1问题的提出 重复以上过程得 3 1问题的提出 如果把原方程写为 构造 3 1问题的提出 例3 2 构造 3 1问题的提出 得 3 1问题的提出 构造 由原方程 3 1问题的提出 得 3 1问题的提出 如果A非奇异 则线性方程组Ax b有唯一解x 将上式化为x Bx f 给出初始向量x0 则有 xk 1 Bxk f k 0 1 2 可以构成一向量序列 xk 若向量序列 xk 收敛于x 则x Bx f 即x 是方程组的解 这种方法称为迭代法 B称为迭代矩阵 3 1问题的提出 构造迭代法的中心问题是建立一个由本次近似值计算下一次近似值的规则 用迭代法求解线性方程组时要解决的问题有 构造一种迭代格式 由xk计算xk 1 证明向量序列 xk 的收敛性 给出初始向量x0 如果序列收敛 证明是原方程组的解 给出估计误差和迭代停止判据 3 1问题的提出 定义 在n维空间中给定一个向量序列 如果对每一个分量 当时都有极限xi 即 则称向量序列有极限 或称收敛于x 3 2雅可比迭代 Jacobiiteration 第五章线性方程组解法 3 2雅可比迭代 最简单的迭代方法是从第i个方程解出未知数xi i 1 2 n 于是雅可比迭代式为 把系数矩阵分解为A U D L 其中U为由A上三角部分构成的上三角阵 L为由A下三角元素构成的下三角阵 D为由A对角线元素构成的对角阵 3 2雅可比迭代 显然 所有aii i 1 2 n不为零上式才有意义 从线性代数知 对于任何系数方阵非奇异的方程组 通过适当交换方程的顺序总可以使所有方程的 3 2雅可比迭代 于是原方程组为 U D L x b 3 2雅可比迭代 上式两边左乘D 1得 x D 1b D 1 U L x Bx f 其中B D 1 U L f D 1b 于是有迭代格式xk 1 Bxk f 例3 3用Jacobi迭代格式解下面方程组 3 2雅可比迭代 解 Jacobi迭代格式为 3 2雅可比迭代 取初始向量x 0 0 0 0 T 各次迭代结果 3 3高斯 塞德尔迭代 Gauss Seideliteration 第五章线性方程组解法 在雅可比迭代中 计算第k 1次迭代近似值时用的是上一次即第k次的近似值 从式 3 3高斯 塞德尔迭代 可以看出 在计算第i个xik 1分量时 前面i 1个分量x1k 1 x2k 1 xi 1k 1已经从上式中计算出来了 于是很自然会想到如果把它们代入用来计算xik 1可能会改进迭代 于是就得到Gauss Seidel迭代格式 写成矩阵形式为 3 3高斯 塞德尔迭代 或 如果 D L 1存在 则 其中 3 3高斯 塞德尔迭代 注记 通常高斯 塞得尔方法比雅可比方法有更快的收敛速度 但不是总这样 对于某些方程组 雅可比迭代收敛 而高斯 塞得尔方法发散 即 并不是任何时候高斯 塞得尔方法都比雅可比方法好 例3 4用Gauss Seidel迭代格式解下面方程组 精确到3位有效数 3 3高斯 塞德尔迭代 解 Gauss Seidel迭代格式如下 3 3高斯 塞德尔迭代 取初始近似值x0 0 0 0 T 各次迭代结果 3 4逐次超松弛迭代法 SOR SuccessiveOverrelaxationMethod 第五章线性方程组解法 逐次超松弛迭代简称SOR方法 是高斯 塞得尔法的一种加速方法 3 4逐次超松弛迭代法 高斯 塞得尔法迭代格式得到 把xik 1改进为xik与的加权平均 即 3 4逐次超松弛迭代法 上式中时 就是高斯 塞得尔方法 为保证迭代过程收敛 要求 当时叫低阶松弛法 当时叫超松弛法 SOR方法收敛时 希望选择一个最佳的使收敛速度最快 目前还没有最佳松弛因子的一般算法理论 实际大都由计算经验通过试算确定的近似值 超松弛迭代式的矩阵形式 3 4逐次超松弛迭代法 直接由分量形式公式写 由 有 所以 证明 3 4逐次超松弛迭代法 由高斯 塞德尔公式推导 3 4逐次超松弛迭代法 高斯 塞德尔迭代公式的矩阵形式是 加权平均 证明 3 4逐次超松弛迭代法 3 5迭代法的收敛性 convergence 第五章线性方程组解法 矩阵的特征值的绝对值最大值称为矩阵A的谱半径 即 3 5迭代法的收敛性 定理 迭代法基本定理 设有方程组x Bx f 对于任意初始向量x0及任意f 迭代公式xk 1 Bxk f收敛的充要条件是 3 5迭代法的收敛性 定理 迭代收敛的充分条件 设有迭代式xk 1 Bxk f 如果 则对于任意初始向量x0 这个迭代过程收敛于方程组x Bx f的唯一解x 并且有事后估计 以及事前估计 3 5迭代法的收敛性 定义 如果对于方阵 有 则称方阵对角占优 定理 如果方程组Ax b的系数阵对角占优 则方程组有唯一解且对任意初始向量x0雅可比迭代和高斯 塞德尔迭代都收敛于真解 3 5迭代法的收敛性 思考题3 1 如何对方程组进行调整 使用Gauss Seidel迭代格式求解时收敛 定理 如果方程组Ax b的系数阵对称正定 则方程组有唯一解且对任意初始向量x0高斯 塞德尔迭代收敛于真解 3 5迭代法的收敛性 Jacobi迭代格式的收敛性 3 5迭代法的收敛性 Jacobi迭代矩阵 特征方程 3 5迭代法的收敛性 由于 所以 注意到 所以 将A的对角线元素乘以后取行列式 令其等于零 就是Jacobi迭代矩阵的特征方程 例3 5讨论用Jacobi迭代格式解方程组的收敛性 3 5迭代法的收敛性 解 Jacobi迭代矩阵的特征方程为 展开得 3 5迭代法的收敛性 Jacobi迭代格式收敛 解得 Gauss Seidel迭代格式的收敛性 3 5迭代法的收敛性 Gauss Seidel迭代矩阵 特征方程 3 5迭代法的收敛性 由于 所以 注意到 所以 将A的对角线以下元素乘以后取行列式 令其等于零 就是Gauss Seidel迭代矩阵的特征方程 3 5迭代法的收敛性 思考题3 2 Gauss Seidel迭代矩阵 至少有1个特征值为零 为什么 例3 6讨论用Gauss Seidel迭代格式解方程组的收敛性 3 5迭代法的收敛性 解 Gauss Seidel迭代矩阵特征方程为 展开得 3 5迭代法的收敛性 Gauss Seidel迭代格式收敛 解得 3 5迭代法的收敛性 作业 1 写出下面方程组的Jacobi迭代格式和Gauss Seidel迭代格式的算法 2 讨论它们的收敛性 3 6高斯消去法 GaussElimination 第三章线性方程组解法 高斯消去法是最古老的数值方法之一 现在仍然是一个很有用的方法 它在计算机上容易实现 3 6高斯消去法 其基本思想是在各个方程之间进行乘法和加减运算 逐步消去方程中的未知数 它分为消去和回代两个过程 给定线性方程组 3 6高斯消去法 第一步消元 令 用乘第一个方程再加到第i个方程上作为第i个新方程 消去x1的项 变为 逐次进行同样过程 最后 经过n 1次消元 得到 3 6高斯消去法 以上这些步骤叫消元过程 3 6高斯消去法 然后 从第n个方程开始 依次解出xn xn 1 x1 3 6高斯消去法 高斯消去法的计算量 3 6高斯消去法 消去过程的计算量 第一步计算乘数mi1 i 2 3 n 需要i 1次除法运算 计算aij 2 i j 2 3 n 需要 i 1 2次乘法运算以及 i 1 2次加减法运算 利用求和公式 得到 3 6高斯消去法 于是消去过程乘除法次数为 消去过程加减法次数为 3 6高斯消去法 计算b n 1 的计算量 乘除法次数为 加减法次数为 3 6高斯消去法 回代计算量 乘除法次数为 加减法次数为 3 6高斯消去法 总的计算量 乘除法次数为 加减法次数为 高斯消去法还有没有办法进行 为什么 3 6高斯消去法 思考题3 3 如果遇到 作业 写出高斯消去法的Fortran程序 3 7高斯主元消去法 GaussEliminationwithpivoting 第三章线性方程组解法 3 7高斯主元消去法 从高斯消去法我们已经看出 为使高斯消去法能顺利进行 必须在每一步消去步都满足条件 但若 相应的系数 计算可能引起很大的舍入误差 为此 需要改进高斯消去法 有两种改进方法 列选主元 完全选主元 3 7高斯主元消去法 列选主元通过交换方程而使得aki i 1 k i i 1 n 中绝对值最大的一个换到 i i 位置而成为第i步的消去主元 完全选主元法就是在系数矩阵的子块 中找出绝对值最大的元素 作为第k 1次消去过程的主元 3 7高斯主元消去法 假设 则k 1行与p行交换 第k 1列与第q列交换 右端也同时交换 在做列交换时 要注意未知量也作交换 即把xk 1与xp交换 3 7高斯主元消去法 列选主元算法 fork 1 2 n 1 找出满足元素的行位置p if error 无唯一解 stoperror if 换行 forj k k 1 n ak j与ap j交换 3 7高斯主元消去法 bk与bp交换 Endfor Endif 消去计算 Endfor 3 7高斯主元消去法 回代计算 fori n 1 n 2 2 1 Endfor 3 7高斯主元消去法 3 8三角分解法 LUdecomposition 第三章线性方程组解法 高斯消去法 3 7三角分解法 第一次消元过程为 用 mi1乘上面第一个方程再加到第i个方程 即可消去第i个方程中的未知量x1 这个过程实际上是给系数矩阵A左乘这样一个下三角阵L1 作第二次消元相当于给系数矩阵A 1 左乘了一个下三角阵L2 3 7三角分解法 作第k次消元相当于系数矩阵A k 1 左乘一个下三角阵Lk 3 7三角分解法 3 7三角分解法 于是第n 1次消元过程可表示为 于是 下三角阵乘积也是一个下三角阵 记 3 7三角分解法 L是单位下三角阵 即对角元素全为1的下三角阵 3 7三角分解法 叫矩阵A的三角分解 或LU分解 如果L为单位下三角阵 则叫Doolittle分解 若U为单位上三角阵 则叫Crout分解 只要A的各顺序主子式不为零 则A可唯一分解成一个单位下三角阵L与一个上三角阵U的乘积 3 7三角分解法 设Ax b A LU 则Ax LUx b 于是令Ux y 则Ly b 3 7三角分解法 这样原来方程能化为两个简单方程组 由求y 然后由求x 3 7三角分解法 现在显式地给出lij和uij的表达式 从第一行可以看出 u1j a1j j 1 2 n 3 7三角分解法 从第一列可以看出ai1 li1u11 i 2 3 n 即li1 ai1 u11 i 2 3 n 从第二行可以看出 a2j l21u1j u2j j 2 3 n 则u2j a2j l21u1j j 2 3 n 从第二列可以看出 ai2 li1u12 li2u22 i 2 3 n 则li2 ai2 li1u12 u22 i 2 3 n 3 7三角分解法 一般地 如果U的前k 1行 L的前k 1行已经求出 则第k k 2 3 n 步 行 3 7三角分解法 列 直到第n步 A全部分解成LU 3 7三角分解法 3 9追赶法 Forwardeliminationandbackwardsubstitution 第三章线性方程组解法 在很多实际问题中 方程组Ax b的系数矩阵A是一个稀疏矩阵 3 6追赶法 假设非零元素集中分布在对角线及其相邻两个次对角线上 且系数阵为对角占优阵 即有 3 6追赶法 把系数矩阵三角分解有 3 6追赶法 利用LU分解公式 写出 3 6追赶法 得到 于是 得到 3 6追赶法 3 6追赶法 3 10其它应用 第三章线性方程组解法 计算行列式值 3 6其他应用 求矩阵逆 分别由 解出xij ij 1 2 n 于是 3 6其他应用 3 11误差分析 Erroranalysis 第三章线性方程组解法 3 6误差分析 定义 设为方程组Ax b的精确解 近似解与精确解之间的误差向量定义为 一种衡量近似解精确度的方法是把近似解代入方程组 看残差向量的大小 于是 3 6误差分析 由范数性质有 而 上面两个式子