线性方程组的数值解法.ppt

n阶线性代数方程组的一般形式为 ： 第三章 线性方程组的数值解法 问题的提出： 写成矩阵-向量形式 若矩阵 非奇异，即 的行列式 ，根据 克莱姆（Gramer）法则，方程组有唯一 解： 其中 为系数矩阵， 为解向量， 为右端常向量。 其中 表示 ， 表示 中第 列换成 后 所得的行列式。 当阶数较高时用这种方法求解是不现实的。 阶行 列式有 项，每项又是 个数的乘积。对较大的 ，其计算量之大，是一般计算机难以完成的。而且 ，这时的舍入误差对计算结果的影响也较大。 例如，求解一个20阶线性方程组，用加减消元法需 3000次乘法运算，而用克莱姆法则要进行 次 运算，如用每秒1亿次乘法运算的计算机要30万年。 线性代数方程组的计算机解法常用方法： 直接法 迭代法 消去法 矩阵三角分解法 直接法：经过有限步算术运算，可求得方程组的 精确解的方法（若在计算过程中没有舍入误差） 迭代法：用某种极限过程去逐步逼近线性方程组 精确解的方法 迭代法具有占存储单元少，程序设计简单，原始 系数矩阵在迭代过程中不变等优点，但存在收敛性 及收敛速度等问题 3.1 消去法 消去法的基本思想：是通过将一个方程乘或除以 某个常数，以及将两个方程相加减，逐步减少方程中 的变元数，最终使每个方程只含一个变元，从而得出 所求的解。 消去法在线性代数中已有详细的讨论，在此只给 出一些说明以及算法的具体描述。 消去法常用方法： 高斯消去法 选主元消去法 高斯约旦消去法 高斯消去法属于直接法，一般由高斯消去法属于直接法，一般由“消元过程消元过程” 和和“回代过程回代过程”两部分组成。先举几个简单实例，两部分组成。先举几个简单实例， 再对一般再对一般n n阶方程组说明高斯消去法的基本思想。阶方程组说明高斯消去法的基本思想。 消去法 3.1 高斯消去法 按自然顺序进行的消元法 例 1 用高斯消元法求解方程组 解 用第一个方程削去后两个方程中的 得 再用第2个方程消去第3个方程中的 得 最后，经过会代求得原方程组的解为 例 2 解方程组 解：消元 回代得 消去法 下面讨论一般下面讨论一般 n n 阶线性方程组的高斯消去法。阶线性方程组的高斯消去法。 记记 为为 ， 和和 的元素的元素 分别记为分别记为 和和 ， ，系数上标，系数上标 代表第代表第1 1次消元之前的状态。次消元之前的状态。 第第1 1次消元时，设次消元时，设 对每行计算乘数对每行计算乘数 用用 乘以第乘以第1 1个方程，加到第个方程，加到第 个方程，消去个方程，消去 第第 2 2个方程到第个方程到第 个方程的未知数个方程的未知数 ，得得 即即: : 其中其中 ： 第第 次消元次消元 时，设第时，设第 次消元已次消元已 完成，即有完成，即有 其中：其中： 设设 ，计算乘数，计算乘数 只要只要 ，消元过程就可以进行下去，直到经，消元过程就可以进行下去，直到经 过过 消元之后，消元过程结束，得消元之后，消元过程结束，得 也即也即 这是一个与原方程组等价的这是一个与原方程组等价的上三角形上三角形方程组。把方程组。把 经过经过 n-1n-1次消元将线性方程组化为上三角形方程组的次消元将线性方程组化为上三角形方程组的 计算过程称为计算过程称为消元过程消元过程。 当当 时，对上三角形方程组自下而上逐步回时，对上三角形方程组自下而上逐步回 代解方程组，计算代解方程组，计算 , ,即即 , ,称为各次消元的主元素，称为各次消元的主元素， 主元素所在的行称为主元素所在的行称为主行主行。 高斯消去法的计算步骤为高斯消去法的计算步骤为： 1 1消元过程消元过程 设设 ，对，对 ，计算，计算 2 2回代过程回代过程 综上所述，高斯消去法的框图如图综上所述，高斯消去法的框图如图3-13-1所示。从所示。从 中可看出高斯消去法的计算机运算和存储方式的特点中可看出高斯消去法的计算机运算和存储方式的特点 ： 1 1按消元规则进行运算后，对角线以下元素为按消元规则进行运算后，对角线以下元素为0 0 。故对于对角线以下元素不用作计算，减小了计算量。故对于对角线以下元素不用作计算，减小了计算量 。 2 2对角线以下元素对回代求解无影响，故可将乘对角线以下元素对回代求解无影响，故可将乘 数放在该处，即数放在该处，即 以节省存储单元。以节省存储单元。 3 3对角线以上元素和常数变换后的元素仍放在原来对角线以上元素和常数变换后的元素仍放在原来 的位置以节省存储单元。的位置以节省存储单元。 4 4回代后的数值仍放在常数项存储单元回代后的数值仍放在常数项存储单元 这时这时 单元中存放的就是输出值单元中存放的就是输出值 定理2 Ax=b 可用高 斯消元法求解的充分必要条件是： 系数矩阵 A 的各阶顺序主子式均不为零。 高斯消元法的条件 定理1 如果在消元过程中A的主元素 (k=1,2,n) ,则可通过高斯消元法求出Ax=b 的解。 引理 A的主元素 (k=1,2,n) 的充要条件 是矩阵A的各阶顺序主子式不为零，即 定理定理：高斯消去法求解：高斯消去法求解 阶线性方程组共需乘除法阶线性方程组共需乘除法 次数近似为次数近似为 。 证明：见书证明：见书P64P64 高斯消去法的计算量 讨讨 论论： 在求解线性方程组时其系数矩阵绝大部分都是非奇在求解线性方程组时其系数矩阵绝大部分都是非奇 异的，但可能出现主元素异的，但可能出现主元素 消去法消去法无法进行 ；或| akk(k) |1时，带来舍入误差的扩散。如何处理 ？ 例 1 解方程组 解法一 用高斯消元法求解（取5位有效数字），用第 一个方程消去第二个方程中的 3.1.2 高斯主元素消元法 因而再回代，得 而精确值为 显然该解与精确值相差 太远，为了控制误差，采用另一种消元过程。 解法二 为了避免绝对值很小的元素作为主元，先交 换两个方程，得到 消去第二个方程中的 得 再回代，解得 结果与准确解非常接近。这个例子告诉我们，在采用高 斯消元法解方程组时，用做除法的小主元素可能使舍入 误差增加，主元素的绝对值越小，则舍入误差影响越大 。固应避免采用绝对值小的主元素，同时选主元素尽量 的大，可使该法具有较好的数值稳定性。 为避免上述错误，可在每一次消元之前增加一为避免上述错误，可在每一次消元之前增加一 个选主元的过程，将绝对值大的元素交换到主对角个选主元的过程，将绝对值大的元素交换到主对角 线的位置。根据交换的方法可分成线的位置。根据交换的方法可分成全选主元全选主元和和列选列选 主元主元两种方法。两种方法。 列主元素消元法 列选主元是当变换到第列选主元是当变换到第k k 步时，从步时，从k k列的列的 及以及以 下的各元素中选取绝对值最大的元素，然后通过行下的各元素中选取绝对值最大的元素，然后通过行交交 换换将其交换到将其交换到 的位置上。交换系数矩阵中的两行（的位置上。交换系数矩阵中的两行（ 包括常数项），相当于两个方程的位置交换了。包括常数项），相当于两个方程的位置交换了。 例：求解线性方程组 解法一：用列主元素消元法，方程组增广矩阵为： 交换交换1 1、3 3行（行（列选主列选主） 消元消元 消元消元 回代计算解为 选主元选主元 全选主元素消元法 全选主元是当变换到第全选主元是当变换到第k k 步时，从右下角步时，从右下角 n-n- k+1k+1阶子阵中选取绝对值大的元素，然后通过阶子阵中选取绝对值大的元素，然后通过行交行交 换换与与列交换列交换将其交换到将其交换到a a kk kk 的位置上，并且保留交 的位置上，并且保留交 换的信息，以供后面调整解向量中分量的次序时换的信息，以供后面调整解向量中分量的次序时 使用。使用。 交换交换1 1、3 3行行 交换交换1 1、3 3列列 回代计算得 从而原方程的解为 消元消元 消元消元 比较而言,Gauss顺序消去法条件苛刻,且 数值不稳定; Gauss全主元消去法工作量偏大 ，算法复杂；对于Gauss-Jordan消去法形式 上比其他消元法简单，且无回代求解，但计 算量大。因此从算法优化的角度考虑， Gauss列主元消去法比较好。 消去法小节 3.2 矩阵三角分解法 3.2.1 3.2.1 矩阵三角分解原理矩阵三角分解原理 应用高斯消去法解应用高斯消去法解 阶线性方程阶线性方程 经过经过 步步 消去后得出一个等价的上三角形方程组消去后得出一个等价的上三角形方程组 ， 对上三角形方程组用逐步回代就可以求出解来。这对上三角形方程组用逐步回代就可以求出解来。这 个过程也可通过矩阵分解来实现。个过程也可通过矩阵分解来实现。 将非奇异阵分解成一个下三角阵将非奇异阵分解成一个下三角阵 和上三角阵和上三角阵 的乘积的乘积 ： 称为称为对矩阵对矩阵 的三角分解的三角分解，又称，又称 分解分解。 高斯消去法解线性代数方程组的一种变形解法。高斯消去法解线性代数方程组的一种变形解法。 杜里特尔分解：杜里特尔分解：把把 分解成一个分解成一个单位下三角阵单位下三角阵 和一个上三角阵和一个上三角阵 的乘积。的乘积。 克劳特分解克劳特分解: : 把把 分解成一个下三角阵分解成一个下三角阵 和一个和一个 单位上三角阵单位上三角阵 的乘积。的乘积。 若矩阵若矩阵 各阶主子式不为零，则可各阶主子式不为零，则可唯一地分唯一地分解成一解成一 个单位下三角阵个单位下三角阵 和一个非奇异的上三角阵和一个非奇异的上三角阵 的的 乘积。乘积。 定定 理理： 证明：略（证明：略（P71P71） 杜里特尔分解 A=LU 杜里特尔分解 杜里特尔分解(一般算法) 由矩阵乘法规则 即即 ： 由此可得 的第一行元素和 的第一列元素 当已得出 的前 行元素和 的前 列元素，则对于 ，由 又可得计算 的第 行元素和 的第 列元 素的公式: 例 求矩阵的LU分解。 u11=2 u12=1 u13=4 所以 a11 a12 a1k a1n u11 u12 u1k u1n 第1框 a21 a22 a2k a2n l21 u22 u2k u2n 第2框 ak1 ak2 akk akn lk1 lk2 ukk ukn 第k框 an1 an2 ank ann ln1 ln2 lnk unn 第n框 按下图所示顺序逐框进行,先求 u 后求 l。 矩阵三角分解算法总结： 有了矩阵 A 的LU分解计算公式，当求解线性方 程组 时，等价于求解 。这时可归结为 利用递推计算相继求解两个三角形方程组（系数矩 阵为三角矩阵）。用顺代，由 求出 ，再 利用回代，由 求出 。 3.2.2 解线性方程组的三角分解法 用杜里特尔三角分解法解线性方程组 的 计算方法： 对于 ，计算 和 。 求解 ，即计算 求解 ，即计算。 计算 和 。 克劳特克劳特分解求方程组步骤分解求方程组步骤：略略 其它矩阵分解法求解特殊的线性方程组： 平方根法：主要用于求解正定矩阵的线性方程组 追赶法：主要用于求解特殊矩阵的三对角方程组 见书 P78P87 用LU 直接分解的方法求解线性方程组的计算量为 ，和高斯消去法的计算量的数量级基本相同。 消去法和矩阵三角分解法比较消去法和矩阵三角分解法比较 当方程组系数矩阵不变，只有右侧向量b变化时， 用LU分解法比消去法速度快。因为右侧向量b的改变 不影响LU的变化。 高斯消去法和LU分解法是等价的，其关键是把一般 方程组变为三角方程组，只是实现途径不同。 设给定 中的向量序列 ， 其中 若对任何 ，都有， 则向量 称为向量序列 的极限， 或者说向量序列收敛于向量 ，记为： 向量收敛的定义： 3.6 3.6 迭代法迭代法 生成向量序列 x(k) ，若 称为迭代格式（1）的迭代矩阵。 则有x* =Bx*+f ， 即x*为原方程组Ax=b 的解，B 迭代法基本思想： 将方程组 Ax=b ( |A|0 ) 转化为与其等价的方程组 x = Bx+f x(k+1) = Bx(k) + f （k=0,1,2,) (1) 取初始向量 x(0)按下列迭代格式 雅可比迭代法 对线性方程组可以建立不同的迭代格式。下面 介绍构造迭代格式的几种常用方法，雅克比迭代 法和高斯塞德尔迭代法。 设方程组 其中 aii(i)0 ( i=1 , 2 , , n) 分别从上式n个方程中分离出n个变量，如下： 等价方程组等价方程组 建立迭代格式 称为雅可比(Jacobi)迭代法又称简单迭代法。 高斯塞德尔迭代法 在 Jacobi 迭代中,用已有的迭代新值代替旧值，即 ： 建 立 迭 代 格 式 或缩写为 称为高斯塞德尔（Gauss Seidel）迭代法。 例1 用雅可比迭代法解方程组 取 计算如下 雅克比迭代格式雅克比迭代格式 Gauss-Seidel 迭代格式为 例2 用GaussSeidel 迭代法解上题 。 取 x(0)=(0,0,0)T 计算如下： 以上介绍的两种迭代法，一般地说，高斯塞德尔以上介绍的两种迭代法，一般地说，高斯塞德尔 比雅克比要好，但也不完全是这样。比雅克比要好，但也不完全是这样。 关于迭代法的一些问题关于迭代法的一些问题： 为了使迭代法计算加速，可采用为了使迭代法计算加速，可采用松弛法松弛法。（略）。（略） 迭代法存在收敛性问题。（略）迭代法存在收敛性问题。（略） 3.3 向量与矩阵的范数 问题的提出 向量范数概念是三维欧氏空间中向量长度概念的 推广,在数值分析中起着重要作用. 为了研究线性方程组近似解的误差估计和迭代法 的收敛性，我们需要对 （n维向量空间）中向量 及 （ 维矩阵空间）中矩阵的“大小”及“距 离”引进某种度量向量与矩阵范数的概念. 平面向量 x 大小：用 距离 反映。 3.3 向量与矩阵的范数 引入范数的目的： 实数大小：用绝对值反映 复数大小：用模反映 高维空间向量 x “大小” 用 反映？ 将度量长度数值的计算方法引入高维空间， 用来反映空间向量的大小，就是范数的概念。 为了研究线性方程组近似解的误差估计和迭为了研究线性方程组近似解的误差估计和迭 代法的收敛性等。代法的收敛性等。 非负性 |x|0 ，等号仅当 x=0 时成立； 齐次性 对任何实数 , | x|=| | |x|； 三角不等式 |x+y| |x| +|y| ； 则称 |x| 为向量 x 的范数。 定义 对任意 n 维向量 x Rn，若对应非负实数 |x| ， 满足 3.3.1 向量的范数 由定义可知，向量 的范数 是按一定规律 与 对应的实数，该实数的值没有确定，但只要满 足这三个条件，这个实数就是向量 的一种范数。 因此定义中的三个条件称为范数公理。 当 时， 向量范数 有如下性质 证：利用条件，有 证： 它们分别叫做向量的-范数、1-范数、2-范数、p - 范数（0p+）。 p -范数是以上范数的统一表达 形式。 常用的向量范数: 满足定义的范数不是唯一的. 设 x = ( x1 , x2 , , xn)T ，常用的向量范数有 对于范数，有 当 时，有 ， 只有当 时，才有 对任意实数 ，因为 ，所以 。 对任意向量 ，有 因此因此 是是n n维空间的一种范数维空间的一种范数 例：范数的证明 不难证明这几种范数满足下述关系 事实上，对 Rn 上任意两种向量范数 |x| , |x| ，都存在与 x 无关的正常数 c1 , c2 使 这种性质称为向量范数的等价性。 向量序列收敛性定理： 向量序列收敛 (即 )的充要条件是 ，其中 为向量的任一种范数。 按不同方式规定的范数,其值一般不同。但在各种 范数下考虑向量序列收敛性的结论是一致的。即向量 序列如对某一种范数收敛则对其它范数也收敛，且有 相同的极限。这样，在研究某一具体问题时，可以选 择一较易简单的范数。 矩阵范数是用于定义矩阵“大小”的量。 由于在大多数与估计误差有关的问题中,矩阵和 向量的乘积经常出现,所以希望引进一种矩阵的 范数,它与向量范数有某种关系。 3.3.2 矩阵的范数 定义：设 ， ，定义矩阵 的范数 对于每一种向量范数 ，相应的矩阵范数 为 其中 是指 的最大可能值。即遍取所有不为0 的 ，比值 中最大的，定义为 的矩阵范数。 矩阵范数的性质： ,且仅当 时， ； 对任意实数 ， ； 对同维方阵 ，有: 相容性条件： 由矩阵范数的定义可直接得到 ，即有 ，称为相容性条件。即所给的 矩阵范数与向量范数是相容的。 证明 p92 矩阵范数与向量范数的关系： 矩阵范数与向量范数密切相关，向量范数有相应 的矩阵范数，即： （如 ） 证明证明 ： 矩阵范数与向量范数的关系： 矩阵范数与向量范数密切相关，向量范数有相应 的矩阵范数，即： （如 ） 常用的矩阵范数有（矩阵范数的求取） 它们分别叫做矩阵的-范数、1-范数。 例4 设则 求相应各范数。 解 验证相容性验证相容性 3.4 方程组的性态 前几节讨论了求解线性代数方程组的直接法.给出 系数矩阵A和自由项b,求未知向量x.实践中,A和b往 往是实验观测数据或是计算所得结果.因此我们处理 的线性方程组 实际上变成了 的关系怎样,是人们十分关心的问题 . 3.4.1 方程组的性态和矩阵的条件数 例 1 解方程组 其中 现用绝对精确的计算（即不带任何舍入误差的计算） 求解，可以看出 此时，我们发现对于两组不同的自由项 它的差只有 而所得解 与 之差却是 两组不同的右端其分量之差不过是 可是解的差高 达 之1860倍像这样的方程组或矩阵A 就叫做病态的 定义1 若方程组 Ax=b 的系数矩阵 A 或右端向量 b 的微小变化（小扰动）引起解产生巨大变化,则称 此方程组为病态方程组; A称为病态矩阵, 否则称为 良态方程组, A 称为良态矩阵。 定理：设 非奇异， ，且 则 为了定量地刻划方程组的“病态”程度，下 面就一般方程组 进行讨论。首先考察右端 项b 的扰动对解的影响。 证： 设x为Ax=b的准确解，当方程组右端有小扰动 b而A准确时,受扰解为 x +x , 即 A(x +x)=b+ b 因为 Ax=b 所以 x=A-1 b 又由 得 因此 所以所以 此不等式表明,解的相对误差不超过b的相对误差的 |A| |A-1| 倍. 若系数矩阵A也有小扰动A ,则还可进一步导出更一 般的误差估计式 定义2 设A 为非奇异矩阵，称数 cond(A)= |A| |A-1| 为A 矩阵的条件数 矩阵的条件数与范数有关，通常使用的条件数有 所以，Cond(A)越大,A的病态程度越严重。 对任何非奇异矩阵A，都有 。 不难证明，条件数具有如下性质 例 求矩阵A 的条件数，其中 解： 于是 从而 所以A是病态的 由于计算条件数涉及到计算逆矩阵，很麻烦。 因此实际计算中一般通过可能产生病态的现象来判断 。 若在消元过程中出现小主元，则 A可能是病态 矩阵。但病态矩阵未必一定有这种小主元。 若解方程组时出现很大的解，则A可能是病