数值分析 第三章 基于MATLAB的科学计算—线性方程组1

科学计算科学计算 理论 方法理论 方法 及其基于 的程序实现与分析及其基于 的程序实现与分析 三 三 解线性方程组 线性矩阵方程 解线性方程组 线性矩阵方程 解线性方程组是科学计算中最常见的问题 所说的解线性方程组是科学计算中最常见的问题 所说的 最常见最常见 有两方面的有两方面的 含义 含义 问题的本身是求解线性方程组 问题的本身是求解线性方程组 许多问题的求解需要或归结为线性方程组的求解 许多问题的求解需要或归结为线性方程组的求解 关于线性方程组关于线性方程组 BAxBAx 1 其其求解方法有两类 求解方法有两类 直接法 直接法 高斯消去法 高斯消去法 GaussianGaussian EliminationElimination 间接法 各种迭代法 间接法 各种迭代法 IterationIteration 高斯消去法高斯消去法 引例引例 考虑如下 梯形 线性方程组 考虑如下 梯形 线性方程组 5 0 141 3 1 5 322 1 1 2 200 430 121 12 143 22 3 32 321 3 2 1 3 32 321 x xx xxx x x x x xx xxx 高斯消去法的求解思路 把一般的线性方程组 化成 上或下 高斯消去法的求解思路 把一般的线性方程组 化成 上或下 梯形的形式 梯形的形式 高斯消去法 高斯消去法 示例示例 考虑如下线性方程组 考虑如下线性方程组 30601 5 1 291 01 2 001 2 2 111 3060129 501 2 001 2 2 1 3 2 1 321 321 321 x x x xxx xxx xxx 第一个方程的两端乘第一个方程的两端乘加到第二个方程的两端 第一个方程的两端乘加到第二个方程的两端 第一个方程的两端乘 1 2 加到第三个方程的两端 得 加到第三个方程的两端 得 30600 3 1 1100 01 0 001 0 0 111 3 2 1 x x x 第二个方程的两端乘第二个方程的两端乘加到第三个方程的两端 得加到第三个方程的两端 得 001 0 10 60600 3 1 10100 01 0 001 0 0 111 3 2 1 x x x 从上述方程组的第三个方程依此求解 得从上述方程组的第三个方程依此求解 得 600 300001 0 31000 24011 3 32 321 x xx xxx 高斯消去法的不足及其改进 高斯消去法的不足及其改进 高斯 全 列 主元素消去法高斯 全 列 主元素消去法 在上例中 由于建模 计算等原因 系数在上例中 由于建模 计算等原因 系数 2 0012 001 而产生而产生 0 00050 0005 的误差 实的误差 实 际求解的方程组为际求解的方程组为 30601 5 1 291 01 2 0005 2 2 111 3 2 1 x x x 70 450 9 3014 2 2565 3 2 1 x x x 注 注 数值稳数值稳定定的算法的算法 高斯列主元素消去法就是在消元的每一步选取 列 主元素高斯列主元素消去法就是在消元的每一步选取 列 主元素 一列中绝对值最大的元取做主元素 高斯列主元素消去法是数值稳一列中绝对值最大的元取做主元素 高斯列主元素消去法是数值稳 定的方法 定的方法 列主元素消去法的基本思想 在每轮消元之前 选列主元素列主元素消去法的基本思想 在每轮消元之前 选列主元素 绝对值最大的元素 绝对值最大的元素 使乘数使乘数 1 ik l 列主元素消去法的步骤 设已经完成第列主元素消去法的步骤 设已经完成第 1 1 步到第步到第步的按列步的按列1 k 选主元 交换两行 消元计算 得到矩阵选主元 交换两行 消元计算 得到矩阵 2 2 2 2 2 2 2 22 1 1 1 1 1 1 1 12 1 11 k n k nn k nk k k k kn k kk nk nk kk baa baa baaa baaaa bA 第第 步计算如下 对于步计算如下 对于 k1 2 1 nk 1 1 选列主元素 即确定选列主元素 即确定 使使 0 i max 0 k ik nik k ki aa 2 2 如果如果 则方程组解不唯一 或者 则方程组解不唯一 或者接近奇异矩阵 接近奇异矩阵 0 0 k ki aA 停止运算 停止运算 3 3 如果如果 则交换 则交换第第 行与第行与第 行元素 行元素 ki 0 bA 0 ik 4 4 消元计算 消元计算 2 1 1 1 nkkjiblbb aalalaa k kik k i k i k kk k ikik k kjik k ij k ij 5 5 回代计算 回代计算 n ij iijijii nnnn nniaxabx abx 1 1 2 2 1 完全主元素消去法即是每次选主元时 依次按行 列选取绝对完全主元素消去法即是每次选主元时 依次按行 列选取绝对 值最大的元素作为主元素 然后交换两行 两列 再进行消元计算值最大的元素作为主元素 然后交换两行 两列 再进行消元计算 完全主元素消去法的步骤 设已经完成第完全主元素消去法的步骤 设已经完成第 1 1 步到第步到第步的选步的选1 k 主元 交换行和列 消元计算 得到矩阵主元 交换行和列 消元计算 得到矩阵 2 2 2 2 2 2 2 22 1 1 1 1 1 1 1 12 1 11 k n k nn k nk k k k kn k kk nk nk kk baa baa baaa baaaa bA 第第 步计算选主元素的范围为步计算选主元素的范围为 即确定 即确定使使knjik 00 j i max 00 k ij njik k ji aa 第第 步计算如下 对于步计算如下 对于 k1 2 1 nk 1 1 选主元素 即确定选主元素 即确定使使 00 j i max 00 k ij njik k ji aa 2 2 如果如果 则方程组解不唯一 或者 则方程组解不唯一 或者接近奇异矩阵 接近奇异矩阵 0 00 k ji aA 停止运算 停止运算 3 3 如果如果 则交换 则交换第第 行与第行与第 行元素 如果行元素 如果ki 0 kk bA 0 ik 则交换 则交换第第列与第列与第 列元素 列元素 kj 0 kk bA 0 jk 4 4 消元计算 消元计算 2 1 1 1 nkkjiblbb aalalaa k kik k i k i k kk k ikik k kjik k ij k ij 5 5 回代求解回代求解 注注 完全主元消去法是解低阶稠密矩阵方程组的有效方法 完全主元消去法是解低阶稠密矩阵方程组的有效方法 但完全主元消去法解方程组 在选主元素时要化费较多的计算机时但完全主元消去法解方程组 在选主元素时要化费较多的计算机时 间 行主元消去法与列主元消去法运算量大体相同 实际计算时 间 行主元消去法与列主元消去法运算量大体相同 实际计算时 用列主元消去法即可满足一定的精度要求用列主元消去法即可满足一定的精度要求 对同一数值问题 用不同的计算方法 所得结果的精度大不一对同一数值问题 用不同的计算方法 所得结果的精度大不一 样样 对于一个算法来说 如果计算过程中舍入误差能得到控制 对计对于一个算法来说 如果计算过程中舍入误差能得到控制 对计 算结果影响较小 则称此算法是数值稳定的 否则 如果计算过程算结果影响较小 则称此算法是数值稳定的 否则 如果计算过程 中舍入误差增长迅速 计算结果受舍入误差影响较大 则称此算法中舍入误差增长迅速 计算结果受舍入误差影响较大 则称此算法 为数值不稳定的为数值不稳定的 因此 我们解数值问题时 应选择和使用数值稳定因此 我们解数值问题时 应选择和使用数值稳定 的算法 否则如果使用数值不稳定的算法 就可能导致计算失败的算法 否则如果使用数值不稳定的算法 就可能导致计算失败 例例 用高斯列主元消去法解方程组用高斯列主元消去法解方程组 0653 0542 132 321 321 321 xxx xxx xxx 解解 0 0 1 653 542 321 bA 1 0 0 13 10 13 20 23 51 1 0 0 321 542 653 31 rr 2 3 1 100 010 001 1 0 0 2 100 2 310 2 101 所以 方程组的解为所以 方程组的解为 T x 2 3 1 高斯列主元素消去法的 高斯列主元素消去法的 MATLAB 实现 实现 意 意xA b 为为 bAx 1 例例 LinearEquiation02 m open LinearEquiation02 LinearEquiation02 一个典型的例子一个典型的例子 HilbertHilbert 矩阵矩阵 1 1 ji H 注注 非奇异矩阵的条件数非奇异矩阵的条件数 分解 分解 Factorization 高斯消去法 高斯消去法 Doolittle 分分 解 解 高斯消去法的消元过程 从代数运算的角度看就是用一个下三角矩阵高斯消去法的消元过程 从代数运算的角度看就是用一个下三角矩阵左左L 乘方程组的系数矩阵 且乘积的结果为上三角矩阵乘方程组的系数矩阵 且乘积的结果为上三角矩阵 即 即U yUx bLy bLUxAx LUA UAL 1 可通过直接用可通过直接用 A 元素计算矩阵元素计算矩阵 A 的三角分解矩阵的三角分解矩阵 L 和和 U 这种这种 直接计算直接计算 A 的三角分解的方法有实用上的好处的三角分解的方法有实用上的好处 下面利用矩阵乘法规下面利用矩阵乘法规 则来确定三角矩阵则来确定三角矩阵 L 和和 U nn n n nn u uu uuu ll ll l A 222 11211 21 3231 21 1 1 1 1 第一步 利用第一步 利用 A 的第一行 第一列元素确定的第一行 第一列元素确定 U 的第一行 的第一行 L 的的 第一列元素第一列元素 由矩阵乘法 由矩阵乘法 2 1 0 0 0 0 0 1 1211 njuuuua j T ijjjj 3 2 0 0 0 1 111111211 niulullla i T iiiii 得到得到 3 73 7 2 1 11 njau jj 3 2 1111 niual ii 设已经计算出设已经计算出 U 的第的第 1 1 至至 r 1 1 行元素 行元素 L 的第的第 1 1 至至 r 1 1 列元素 列元素 现在要计算现在要计算 U 的第的第 r 行元素及行元素及 L 的第的第 r 列元素列元素 第第 r 步 利用步 利用 A 的第的第 r 行 第行 第 r 列剩下的元素确定列剩下的元素确定 U 的第的第 r 行 行 L 的第的第 r 列元素列元素 由矩阵乘法 有由矩阵乘法 有 3 2 nr T jjjjrrrrrj uuullla 0 0 0 0 1 21121 1 1 1 nrrjuul rjkj r k rk 得得 U 的第的第 r 行元素为行元素为 3 83 8 1 1 3 2 1 r k kjrkrjrj nrnrrjulau 由由 T rrrriiiiir uuullla 0 0 0 0 1 21121 2 1 1 1 nrriulul rrirkr r k ik 得得 3 93 9 1 1 3 2 1 1 nrinruulal rrkr r k ikirir 例例 5 5 用用 LU 分解法求解方程组分解法求解方程组 7 1 3 542 774 322 3 2 2 x x x 解解 对系数矩阵对系数矩阵 A 进行进行 LU 分解 分解 1 2 3 2 2 3121131211 lluuu 由由 有 有 jjj ulau 12122 1 3 2322 uu 2 2212313232 uulal6 233213313333 ululau 因此因此 6 13 322 121 12 1 LUA 解方程组解方程组 得 得 bLy 627 521 3 213121 yyyyyy 解方程组解方程组 得 得 yUx 22 323 23 5 1 321323 xxxxxx 6 LU 分解的分解的 MATLAB 实现 实现 或或 AluUL AluPUL 例例 A rand 5 A rand 5 L U L U P lu A P lu A A rand 5 L U P lu A L P L 当当是是主对角占优的三对角矩阵主对角占优的三对角矩阵时 基于时 基于 Doolittle 分解可得到分解可得到A 解这类方程组的解这类方程组的追赶法追赶法 Cholesky 分解分解 Cholesky Factorization 对称正定矩阵对称正定矩阵的的 CholeskyCholesky 分解和以分解和以为系数矩阵地的线性方为系数矩阵地的线性方AA 程组的程组的改进的平方根法改进的平方根法 设设 阶方程组阶方程组 是对称正定矩阵是对称正定矩阵 Positive Positive DefiniteDefinite nbAx A Matrix Matrix 则 则有三角分解有三角分解A nn n n nn u uu uuu ll ll l LUA 222 11211 21 3231 21 1 1 1 1 再将再将分解为分解为U nn n n u uu uuu U 222 11211 0 22 2 11 1 11 12 22 11 1 1 1 UD u u u u u u u u u n n nn 则则 0 UDLA 1 1 对称正定矩阵有唯一的分解对称正定矩阵有唯一的分解 T LDLA 这是由于这是由于 且对称阵 且对称阵 则有 则有 0 UDLA TTT DLUA 0 AAT TT DLULDU 00 再利用三角分解的唯一性 得再利用三角分解的唯一性 得 因此 对称正定矩阵有唯一的因此 对称正定矩阵有唯一的LU T 0 分解分解 T LDLA 2 2 是正定对角阵 即是正定对角阵 即 Dniuii 2 1 0 由于由于对称正定的充要条件是对称正定的充要条件是对称正定 其中对称正定 其中是是 阶可逆阶可逆A T BABBn 方阵方阵 取取 就推知 就推知是正定对角阵是正定对角阵 T LB D 因此因此的对角元素的对角元素 记 记 其中 其中D 2 1 0niuii 2 1 2 1 DDD nn u u u D 22 11 2 1 则则 TT LDLDLDLDA 2 1 2 1 2 1 2 1 3 3 乔莱斯基 乔莱斯基 CholeskyCholesky 分解 分解 将将记为记为 则 则称为称为 CholeskyCholesky 分解分解 利用利用 CholeskyCholesky 直直 2 1 LDL T LLA 接分解公式 推导出的解方程组方法 称为接分解公式 推导出的解方程组方法 称为 CholeskyCholesky 方法或平方根方法或平方根 法法 4 4 解方程组的平方根法 解方程组的平方根法 CholeskyCholesky 方法 方法 由由 CholeskyCholesky 分解 有分解 有 3 103 10 T nn n n nnnn LL l ll lll lll ll l A 222 12111 21 2221 11 利用矩阵乘法 逐步确定利用矩阵乘法 逐步确定 的第的第 行元素行元素 Li 2 1 ijlij 由由 当 当时 时 有分解公式 对于 有分解公式 对于 1 1 j k jjijjkikij llllakj 0 jk l ni 2 1 3 113 11 1 1 2 1 2 1 1 1 2 1 i k ikiiii j k jjjkikijij lalijlllal 将对称正定矩阵将对称正定矩阵作作 CholeskyCholesky 分解后 则解方程组分解后 则解方程组就转化为解两个三角方程就转化为解两个三角方程AbAx 组组 yxLbLy T 例例 7 7 用用 CholeskyCholesky 方法解方程组方法解方程组 25 7 6 4 5 375 2 1 75 2 25 4 1 114 x 解解 对系数矩阵作对系数矩阵作 CholeskyCholesky 分解得到分解得到 1 5 12 5 05 02 15 15 0 25 0 2 5 375 2 1 75 225 4 1 114 T LL 解解 得 得 bLy T y 1 5 3 2 解解 得 得 yxL T T x 1 1 1 choleskycholesky 分解的分解的 MATLABMATLAB 的实现 的实现 L chol A L chol A 3 追赶法追赶法 在许多实际问题中 如 常微分方程两点边值问题 三次样条在许多实际问题中 如 常微分方程两点边值问题 三次样条 插值方法等 往往遇到线性方程组插值方法等 往往遇到线性方程组的求解 其中的求解 其中fAx 3 133 13 T n nn nnn fff ba cba cba cb A 1 111 222 11 称具有公式 称具有公式 3 133 13 形式的系数矩阵 形式的系数矩阵为为三对角阵三对角阵 称相应的线性方称相应的线性方A 程组为程组为三对角方程组三对角方程组 TridiagonalTridiagonal LinearLinear Systems Systems 具有这种形具有这种形 式的方程组在实际问题中是经常遇到的 而且往往是对角占优式的方程组在实际问题中是经常遇到的 而且往往是对角占优 DiagonallyDiagonally DominantDominant 的 的 满足条件 满足条件 A 0 11 cb 1 3 2 0 0 nicacab iiiii 0 nn ab 这类方程组这类方程组的解存在唯一 的解存在唯一 非奇异 可以直接利用高斯消非奇异 可以直接利用高斯消fAx A 去法或直接分解法 而其解答可以用极其简单的递推公式表示出来 去法或直接分解法 而其解答可以用极其简单的递推公式表示出来 即下面介绍的追赶法即下面介绍的追赶法 追赶法通常是数值稳定的追赶法通常是数值稳定的 对对作作 LU 分解 分解 DoolitleDoolitle 分解 可以发现分解 可以发现 L U 具有非常具有非常A 简单的形式简单的形式 n n n l u ul ul m m m ULA 1 22 11 3 2 1 1 1 1 由矩阵乘积 得由矩阵乘积 得 nnnnn n nn nnn lumlm u ulumlm ul ba cba cba cb A 11 1 221212 11 111 222 11 比较等式两端 得到比较等式两端 得到 3 143 14 3 2 3 2 1 2 1 1 1 11 niumbl nilam bl nicu iiii iii ii 因为上述分解因为上述分解 则方程组 则方程组的求解转化为解两个简单的求解转化为解两个简单LUA fAx 的三角方程组的三角方程组和和 从而得到求解方程组的算法公式 从而得到求解方程组的算法公式 fLy yUx 先解先解 即 即fLy 3 153 15 3 2 111 niymfyfy iiii 再解再解 即 即yUx 3 163 16 1 2 1 1 nnilxuyxlyx iiiiinnn 这种把三对角方程组的解用递推公式 这种把三对角方程组的解用递推公式 3 143 14 3 153 15 3 163 16 表示出来的方法形象化地叫做追赶法 其中 表示出来的方法形象化地叫做追赶法 其中 3 143 14 3 153 15 是关于下标 是关于下标 由小到大的递推公式称为追的过程 而 由小到大的递推公式称为追的过程 而 1616 i 却是下标却是下标 由大到小的递推公式称为赶的过程 一追一赶构成了求解由大到小的递推公式称为赶的过程 一追一赶构成了求解i 的追赶法的追赶法 fAx 例例 9 9 用追赶法解三对角方程组用追赶法解三对角方程组 2 3 1 41 141 14 x 解解 系数矩阵分解得到系数矩阵分解得到 7333 3 175 3 14 12667 0 125 0 1 UL 解解 得 得 fLy T y 8667 2 25 3 1 解解 得 得 yUx T x 7679 0 0714 1 5179 0 调用函数调用函数 LU Factorization m 解例解例 9 9 输入输入 A 4 1 0 1 4 1 0 1 4 b 1 3 2 x L U index LU Factorization A b 得到方程组的解及相应的得到方程组的解及相应的 LU 分解矩阵 分解矩阵 x 0 5179 L 1 0000 0 0 U 4 0000 1 0000 0 1 0714 0 2500 1 0000 0 0 3 7500 1 0000 0 7679 0 0 2667 1 0000 0 0 3 7333 为了对线性方程组的直接法作出误差分析 为了讨论方程组迭代法为了对线性方程组的直接法作出误差分析 为了讨论方程组迭代法 的收敛性 需要对向量和矩阵的大小进行度量 进而引入了的收敛性 需要对向量和矩阵的大小进行度量 进而引入了 范数范数 用于度量用于度量 量量 的大小的概念的大小的概念 引言引言 实数的绝对值 实数的绝对值 是数轴上的点是数轴上的点到原点到原点的距离 的距离 aa0 复数的模 复数的模 是平面上的点是平面上的点到原点到原点的距离 的距离 22 babia ba 0 0 还有其他刻画复数大小的方法 准则 如还有其他刻画复数大小的方法 准则 如 ba ba max 向量的内积 范数及向量的内积 范数及维空间距离的度量维空间距离的度量n 令令是一数域 是一数域 是是上的向量空间 如果函数上的向量空间 如果函数PPnP 有如下性质 有如下性质 xy PPP nn 共轭对称性 共轭对称性 xyPn yxxy 非负性 非负性 xPn xx 0 xxx 00 线性性 线性性 xyzPn abP axbyzaxzbyz 则称则称是是上的一个向量内积 上的一个向量内积 innerinner productproduct 向量空间 向量空间 xy Pn 上的向量内积通常用符号上的向量内积通常用符号表示 定义了内积的向量空间表示 定义了内积的向量空间Pn xy 称为内积空间 称为内积空间 innerinner productproduct spacespace 记做记做表示 表示 Pn Pn 例 例 容易验证函数 容易验证函数 QQP Tn n Q 0 xyPn xyx Qy T 定义了定义了上的一个内积 上的一个内积 Pn 令令是一数域 是一数域 是是上的向量空间 如果函数上的向量空间 如果函数有有PPnP x PP n 如下性质 如下性质 非负性 非负性 xPn x 0 xx 00 齐次性 齐次性 xPn aP axax 三角不等式 三角不等式 xyPn xyxy 则称则称是是上的一个向量范数 上的一个向量范数 normnorm 向量空间 向量空间上的范数通上的范数通 xPnPn 常用符号常用符号表示 定义了范数的向量空间表示 定义了范数的向量空间称为赋范空间 称为赋范空间 normednormed xPn spacespace 记做 记做表示 表示 Pn 例 例 容易验证函数 容易验证函数 QQP Tn n Q 0 xPn xx Qx T 定义了定义了上的一个范数 这样定义的范数称为由内积 诱上的一个范数 这样定义的范数称为由内积 诱Pn 导的范数 导的范数 例 例 上常用的向量范数 上常用的向量范数 RC nn xxxxn T 12 RC nn 范数 范数 xxk k n 1 1 范数 范数 xx x T 2 范数 范数 xxxx k n n max 1 12 令令是一数域 是一数域 是是上的向量空间 如果实值函数上的向量空间 如果实值函数PPnP 有如下性质 有如下性质 d xy PPR nn 对称性 对称性 xyPn d yxd xy 非负性 非负性 xyPn d xy 0 d xy 0 xy 三角不等式 三角不等式 xyzPn d xz d xy d yz 则称则称是是上的一个距离 函数 上的一个距离 函数 distancedistance functionfunction 或 或 d xy Pn 度量 度量 metricmetric 定义了度量的向量空间 定义了度量的向量空间称为度量空间 称为度量空间 metricmetric Pn spacespace 记做 记做表示 表示 Pd n 例 例 4 4 上常用的 由范数诱导的 度量 上常用的 由范数诱导的 度量 RC nn xy RC nn 范数诱导的度量 范数诱导的度量 d xyxy 1 范数诱导的度量 范数诱导的度量 d xyxy 2 范数诱导的度量 范数诱导的度量 d xyxy 矩阵的范数矩阵的范数 矩阵矩阵是线性映射 当是线性映射 当时为线性变换 时为线性变换 的的APn m nm PP mn 一种表现形式 因此 除了可以把矩阵看做向量而定义其范数外 一种表现形式 因此 除了可以把矩阵看做向量而定义其范数外 更为基本 更为重要的是表征其线性映射的算子范数 更为基本 更为重要的是表征其线性映射的算子范数 operatoroperator normnorm 以 以的情况为例 的情况为例 APn n x Ax A x 0 sup 其中 右端的范数其中 右端的范数是赋范空间是赋范空间中向量的范数 中向量的范数 Pn 由矩阵算子范数的定义 容易证明 对映像大小的估计 由矩阵算子范数的定义 容易证明 对映像大小的估计 不等式 不等式 AxA x xPn 称满足不等式 的矩阵范数是与对应的向量范数相容的 称满足不等式 的矩阵范数是与对应的向量范数相容的 例 例 常用的矩阵范数 常用的矩阵范数 范数 列范数 范数 列范数 Aajn j n ij i n 1 1 1 12 max 范数 谱范数 范数 谱范数 AA A T 2 max 范数 行范数 范数 行范数 Aain i n ij i n max 1 1 12 上述三种范数是如下定义的矩阵上述三种范数是如下定义的矩阵 范数的特例 范数的特例 p 由向量的 由向量的 范数 范数 定义 定义 pxx pk p k n p 1 1 1 p AAx p x p sup 0 范数 范数 FrobeniusFrobenius Aa Fij j n i n 2 11 1 2 例例 设设 求范数 求范数 43 21 A F AAAA 21 解解 6 42 31max 1 A 7 43 21max A 4772 5 3016941 F A 由于由于 所以 所以 2014 1410 AAT 因此 因此 2211522115 22115max AAT 4650 5 22115 2 A 向量和矩阵范数的向量和矩阵范数的 MATLABMATLAB 实现实现 在在 MATLABMATLAB 中 中 norm norm 函数求向量与矩阵的范数 其命令格式函数求向量与矩阵的范数 其命令格式 为为 norm X p norm X p 当当 X X 为向量或矩阵时 为向量或矩阵时 norm X p norm X p 表示向量或矩阵表示向量或矩阵 X X 的的 p p 范数 范数 例如 例如 norm X 1 norm X 1 表示表示 X X 的的 1 1 范数 范数 norm X 2 norm X 2 表示表示 X X 的的 2 2 范数 范数 norm X inf norm X inf 表示表示 X X 的的范数范数 对于矩阵对于矩阵 X X norm X fro norm X fro 表示表示 X X 的的 F F 范数范数 矩阵的条件数与直接法的误差分析矩阵的条件数与直接法的误差分析 解线性方程组的直接法产生误差的主要原因 解线性方程组的直接法产生误差的主要原因 1 不同的算法及 不同的算法及 舍入误差的影响 舍入误差的影响 2 方程组本身固有的问题 病态或良态 本节我 方程组本身固有的问题 病态或良态 本节我 们将分析方程组的状态并估计算法的误差 即原始数据扰动对解的们将分析方程组的状态并估计算法的误差 即原始数据扰动对解的 影响影响 考虑考虑 阶线性方程组阶线性方程组 其中 其中为非奇异矩阵为非奇异矩阵 nbAx A 由于由于 或 或 的数值是测量得到的 或者是计算的结果 在第 的数值是测量得到的 或者是计算的结果 在第Ab 一种情况下一种情况下 或 或 常带有某些观测误差 在后一种情况 常带有某些观测误差 在后一种情况 或 或 AbAb 又包含有舍入误差又包含有舍入误差 因此我们处理的实际矩阵是因此我们处理的实际矩阵是 或 或 AA bb 下面我们来研究数据下面我们来研究数据 或 或 的微小误差对解的影响 的微小误差对解的影响 首先考虑一首先考虑一Ab 个例子个例子 例例 设方程组设方程组 即 即 它的精确解为 它的精确解为 bAx 00001 8 8 00001 6 2 62 2 1 x x T 1 1 现在考虑系数矩阵和右端项的微小变化对方程组解的影响 即现在考虑系数矩阵和右端项的微小变化对方程组解的影响 即 考察方程组考察方程组 00002 8 8 99999 5 2 62 2 1 x x 其解变为其解变为 这种现象的出现是完全由方程组的性态决定的这种现象的出现是完全由方程组的性态决定的 T 2 10 定义 定义 如果矩阵如果矩阵或常数项或常数项 的微小变化 引起方程组的微小变化 引起方程组解的巨解的巨AbbAx 大变化 则称此方程组为大变化 则称此方程组为病态方程组病态方程组 矩阵 矩阵称为称为病态矩阵病态矩阵 相对 相对A 于方程组而言 否则称方程组为于方程组而言 否则称方程组为良态方程组良态方程组 矩阵 矩阵称为称为良态矩良态矩A 阵阵 需要一种能刻划矩阵和方程组需要一种能刻划矩阵和方程组 病态病态 程度的标准程度的标准 线性方程组的误差分析线性方程组的误差分析 设线性方程组为设线性方程组为 bAx 其中其中 且 且非奇异非奇异 为精确解 为精确解 为解的误差 记为解的误差 记 nn RA n Rbx A xx xxx 设设为为的误差 的误差 为为 的误差的误差 下面分别讨论下面分别讨论与与 的关的关A Ab bx A b 系系 b 有误差而有误差而 A 无误差情形无误差情形 将带有误差的右端项和带误差的解向量代入方程组 则将带有误差的右端项和带误差的解向量代入方程组 则 bbxxA 由于由于 而得到 而得到 从而 从而 bAx bAx 1 bAx 1 另一方面 由另一方面 由 3 24 3 24 式取范数 有式取范数 有 即 即 可得可得xAb 0 1 b b A x 定理定理 设设是非奇异矩阵 是非奇异矩阵 且 且 则有误差 则有误差A0 bAxbbxxA 估计式估计式 b b Acond x x 其中其中称为称为方阵方阵 A A 的条件数的条件数 AAAcond 1 注注 1 1 解的相对误差是右端项解的相对误差是右端项 的相对误差的的相对误差的 cond A cond A b 倍 倍 2 2 如果条件数越大 则解的相对误差就可能越大 如果条件数越大 则解的相对误差就可能越大 3 3 条件数成了刻划矩阵的病态程度和方程组解对条件数成了刻划矩阵的病态程度和方程组解对或或 扰动扰动Ab 的敏感程度的敏感程度 定义定义 7 7 称条件数很大的矩阵为称条件数很大的矩阵为 病态病态 矩阵 称病态矩阵矩阵 称病态矩阵 对应的方程组为病态方程组对应的方程组为病态方程组 反之 则称矩阵为良态矩阵 对应的方反之 则称矩阵为良态矩阵 对应的方 程组为良态方程组程组为良态方程组 A 及及 b 都有误差的情形都有误差的情形 定理定理 7 7 设方程组设方程组 为非奇异矩阵 为非奇异矩阵 及及 都有误差 若都有误差 若0 bAxAAb 的误差的误差非常小 使非常小 使 则有误差估计式 则有误差估计式AA 1 1 AA 3 26 3 26 b b A A A A Acond Acond x x 1 例例 已知方程组已知方程组中中bAx 2 10 60 37 60 47 12 13 6 1 5 1 4 1 5 1 4 1 3 1 4 1 3 1 2 1 bA 若若时 估计解的相对误差时 估计解的相对误差 T b 001 0 001 0 001 0 解解 由于由于由式 由式 3 253 25 有误差 有误差 001 0 108 3333 10 12 13 2 bb 估计估计 0186 0 10 12 13 001 0 2015 2 x x 比右端项的相对误差比右端项的相对误差扩大了扩大了 20152015 倍倍 x x 6 102308 9 b b 例例 下列希尔伯特 下列希尔伯特 HilbertHilbert 矩阵是一族著名的病态矩阵 矩阵是一族著名的病态矩阵 12 1 1 1 1 1 13 12 1 12 11 nnn n n Hn 用用 MATLABMATLAB 函数计算条件数函数计算条件数 输入输入 n H forfor n 3 8n 3 8 cond hilb n cond hilb n endend 得到得到 3 3 至至 8 8 阶希尔伯特矩阵的条件数分别为 阶希尔伯特矩阵的条件数分别为 524 0568524 0568 1 5514e 0041 5514e 004 4 7661e 0054 7661e 005 1 4951e 0071 4951e 007 4 7537e 0084 7537e 008 1 5258e 0101 5258e 010 由此可见 随着由此可见 随着 的增加 的增加 的病态可能越严重的病态可能越严重 常出现在常出现在n n H n H 数据拟合和函数逼近中数据拟合和函数逼近中 对于病态方程组 通常的方法无法得到它的准确解 需要采用对于病态方程组 通常的方法无法得到它的准确解 需要采用 一些特殊的处理方法一些特殊的处理方法 病态方程组的处理病态方程组的处理 对于病态方程组 可采用高精度的算术运算 如双精度或扩充对于病态方程组 可采用高精度的算术运算 如双精度或扩充 精度 以改善或减轻方程组的病态程度精度 以改善或减轻方程组的病态程度 如果用无限精度运算即不存如果用无限精度运算即不存 在舍入误差的话 即使条件数很大 也没有病态可言在舍入误差的话 即使条件数很大 也没有病态可言 我们也可对病我们也可对病 态方程组作预处理 使改善方程组系数矩阵的条件数态方程组作预处理 使改善方程组系数矩阵的条件数 例例 设方程组设方程组 即 即bAx 6 666 10 3 2000 1 10 210 310 4 10