解线性方程组的迭代法

1 解线性方程组的迭代法解线性方程组的迭代法 Haha 送给需要的学弟学妹 摘要摘要 因为理论的分析表明 求解病态的线性方程组是困难的 但是实际情况是否如此 需要我们来具体 检验 系数矩阵 H 为 Hilbert 矩阵 是著名的病态问题 因而决定求解此线性方程组来验证上述问Hxb 题 详细过程是通过用 Gauss 消去法 J 迭代法 GS 迭代法和 SOR 迭代法四种方法求解线性方程组 Hxb 关键词关键词 病态方程组 Gauss 消去法 J 迭代法 GS 迭代法 SOR 迭代法 目录 目录 一 问题背景介绍一 问题背景介绍 二 建立正确额数学模型二 建立正确额数学模型 三 求解模型的数学原理三 求解模型的数学原理 1 Gauss 消去法求解原理 2 Jacobi 迭代法求解原理 3 G S 迭代法求解原理 4 SOR 迭代法求解原理 5 Jacobi 和 G S 两种迭代法收敛的充要条件 四 计算过程四 计算过程 一 一 Hilbert 矩阵维数矩阵维数 n 6 时时 1 Gauss 消去法求解 2 Jacobi 迭代法求解 3 G S 迭代法求解 4 SOR 迭代法求解 二 二 Hilbert 矩阵维数矩阵维数 n 20 50 和和 100 时时 1 G S 迭代法求解图形 2 SOR 迭代法求解图形 五 编写计算程序五 编写计算程序 六 解释计算结果六 解释计算结果 1 Gauss 消去法误差分析 2 G S 迭代法误差分析 3 SOR 迭代法误差分析 G S 迭代法与 SOR 迭代法的误差比较 七 心得体会七 心得体会 正文正文 一 问题背景介绍 一 问题背景介绍 理论的分析表明 求解病态的线性方程组是困难的 实际情况是否如此 会出现怎样的现象呢 二 建立正确的数学模型 二 建立正确的数学模型 考虑方程组的求解 其中系数矩阵 H 为 Hilbert 矩阵 Hxb 1 1 2 1 i jn ni j Hhhi jn ij 这是一个著名的病态问题 通过首先给定解 为方便计算 笔者取 x 的各个分量等于 1 再计算出右端 这样的解就明确了 再用 Gauss 消去法 J 迭代法 GS 迭代法和 SOR 迭代法四种方法分 bHx Hxb 别求解将求解结果与给定解比较 而后求出上述四种方法的误差 得出哪种方法比较好 Hxb 三 求解模型的数学原理 三 求解模型的数学原理 1 Gauss 消去法求解原理消去法求解原理 对于 A 非奇异 求解时 可以先将 A 分解成一个下三角矩阵 L 和一个上三角矩阵 U 的乘积 Axb 2 即 就可以通过ALU MERGEFORMAT 1 1 LybUxy 求解出的值 x 接下来就具体讲讲如何将 A 分解成 L 和 U 也就是 Gauss 消去法 欲把一个给定的矩阵 A 分解为一个下三角阵 L 与一个上三角阵 U 的乘积 最自然的做法便是通过一 系列的初等变换 逐步将 A 约化为一个上三角阵 而又能保证这些变换的乘积是一个下三角阵 这可归 结为 对于一个任意给定的向量找一个尽可能简单的下三角阵 使经这一矩阵作用之后的第 n xR x 至第个分量均为零 能够完成这一任务的最简单的下三角阵便是如下形式的初等下三角阵 1k n T kkk LIl e 其中 1 0 0 T kkknk lll 即 1 1 1 1 1 k kk n k L l l 这种类型的初等下三角阵称作 Gauss 变换 而称向量为 Gauss 向量 k l 对于一个给定的向量我们有 1 T n n xxxR 111 T kkkk kknk nk L xxxxx lxx l 由此立即可知 只要取 1 i ik k x likn x 便有 1 0 0 T kk L xxx 当然 这里我们要求0 k x 而后经过多次变换可以得到 MERGEFORMAT 1 2 11 nn L LL AU 从而求出上三角阵 U 而后通过 11 nn LAULL LL 求得下三角阵 MERGEFORMAT 1 3 1 LUA 将 1 2 和 1 3 带入到 1 1 式中求出的值即可 x 2 J 迭代法求解原理迭代法求解原理 考虑非奇异线性代数方程组 Axb 令 MERGEFORMAT 1 4 ADLU 其中 1122 ijnn AaDdiag aaa 3 21 3132 12 1 0 0 0 0 nnn n a aaL aaa 12131 232 1 0 0 0 0 n n nn aaa aa U a 那么 1 4 式和合并后可以写成Axb xBxg 其中 MERGEFORMAT 1 5 11 BDLUgD b 若给定初始向量 000 012 T n xxxx 并代入 1 4 式右边 又可得到一个向量 一次类推 有 2 x MERGEFORMAT 1 6 1 1 2 kk xBxgk 这就是所谓的 Jacobi 迭代法 其中叫做 Jacobi 迭代法的迭代矩阵 叫做 Jacobi 迭代法的常数项 Bg 3 GS 迭代法求解原理迭代法求解原理 注意到 Jacobi 迭代法中各分量的计算顺序是没有关系的 先算那个分量都一样 现在 假设不按 Jacobi 迭代格式 而是在计算的第一个分量用的各个分量计算 但当计算的第二个分量时 k x 1k x k x 2 k x 因已经算出 用它代替 其他分量仍用 类似地 计算时 因都已算出 1 k x 1 1 k x 1k i x k l x 11 kk l xx 用它们代替其他分量仍用的分量 于是有 11 11 kk l xx 1k x MERGEFORMAT 1 7 11 1 1 2 kkk xD LxD Uxgk 我们称这种迭代格式为 Gauss Seidel 迭代法 简称为 G S 迭代法 它的一个明显的好处是在编写程序是存 储量减少了 如果存在 G S 迭代法可以改写成 1 DL MERGEFORMAT 1 8 11 1 kk xDLUxDLb 我们把叫做 G S 迭代法的迭代矩阵 而把叫做 G S 迭代法的常数项 1 1 LDLU 1 DLb 4 SOR 迭代法求解原理迭代法求解原理 我们知道 G S 迭代法的迭代格式为 111 11 kkk xD LxD UxD b 现在令则有 1 kk xxx MERGEFORMAT 1 9 1 kk xxx 这就是说 对 G S 迭代法来说 可以看作在向量上加上修正项而得到的 若修正项的前面加上 1k x k xx 一个参数 便得到松弛迭代法的迭代格式 MERGEFORMAT 1 10 1 111 1 1 kk kkk xxx xD LxD UxD b 4 用分量形式表示即为 MERGEFORMAT 1 11 1 11 11 1 in kkkk iiijjijji jj i xxb xb xg 其中叫做松弛因子 当时 相应的迭代法叫做超松弛迭代法 当时 叫做低松弛迭代法 1 1 当时 就是 G S 迭代法 我们把超松弛迭代法简称为 SOR 迭代法 1 因为存在 所以 1 10 式可以改写为 1 1 ID L 1 1 kk xL xDLb 其中 1 1LDLDU 叫做松弛迭代法矩阵 而 SOR 迭代法收敛的充要条件是 MERGEFORMAT 1 12 1 L 由 1 12 式知 SOR 迭代法的谱半径依赖于 当然会问 能否适当选取使收敛速度最快 这就是 选择最佳松弛因子的问题 经过相关计算可知 随着从 0 增加 减少 直至 L MERGEFORMAT 1 13 2 2 11 opt B B 为Jacobi 迭代矩阵 时 达到极小 L MERGEFORMAT 1 14 2 2 11 11 opt B L B 再增加时 开始增加 因此 称为最佳松弛因子 L opt 5 Jacobi 和和 G S 两种迭代法收敛的充要条件两种迭代法收敛的充要条件 Jacobi 迭代法和 G S 迭代法两种迭代法有一个共同的特点 那就是新的近似解是已知近似解的 k x 1k x 线性函数 并且只与有关 即它们都可以表示成如下形式 1k x MERGEFORMAT 1 15 1 kk xMxg 事实上 对 Jacobi 迭代法 有 11 MDLUgD b 对 G S 迭代法 有 11 MDLUgDLb 故要求出上述两个迭代法中有确定的解 且与相对误差较小 就必须说明用上述两种迭代 k x 1k x k x 法求解时收敛 下面就给出关于上述两种迭代法收敛的证明原理 解方程组的单步线性定常迭代法 1 15 收敛的充分必要条件是其迭代矩阵的谱半径小于 1 即M MERGEFORMAT 1 16 1 M 从上述的 1 16 可知 迭代序列收敛取决于迭代矩阵的谱半径 而与初始向量的选取和常数项无关 四 计算过程 四 计算过程 方程组的求解 其中系数矩阵 H 为 Hilbert 矩阵 Hxb 1 1 2 1 i jn ni j Hhhi jn ij 一 一 求解 我们暂时选择系数矩阵 H 的维数 所以Hxb 6n 5 11111 1 23456 111111 234567 111111 345678 111111 456789 111111 5678910 111111 67891011 H 令 x 的各分量都为 1 即 1 1 1 1 1 1 T x 根据得Hxb 2 450 1 593 1 218 0 996 0 846 0 737 bHx 而后我们接下来就用 Gauss 消去法 J 迭代法 GS 迭代法和 SOR 迭代法四种迭代法求解 MERGEFORMAT 1 17 11111 1 23456 111111 2 450 234567 1 593 111111 1 218 345678 1111110 996 4567890 846 111111 0 737 5678910 111111 67891011 x 的解 x 1 Gauss 消去法求解消去法求解 因为 6 11111 1 23456 111111 234567 111111 345678 111111 456789 111111 5678910 111111 67891011 H 所以由 1 2 式知 H 分解的上三角矩阵 3333 444 55 6 10 50 3330 250 20 167 0 0830 0830 0750 0670 06 5 556 108 333 109 524 109 921 10 3 571 107 143 109 921 10 2 268 105 669 10 1 432 10 U 由 1 3 式求出 H 的下三角矩阵 1 1 0 51 0 33311 0 250 91 51 0 20 81 71421 0 1670 7141 7862 7782 51 LA U 再通过 1 1 求出 x 的值 1 3 5 6 2 45 0 368 0 033 2 063 10 7 937 10 1 432 10 yLb 1 1 000 1 000 1 000 1 000 1 000 1 000 xUy 所以 1 0001 0001 0001 0001 0001 000 T x 7 2 Jacobi 迭代法求解迭代法求解 将系数矩阵 H 用 1 4 方法分解成 1 1 1 1 1 1 3 5 7 9 11 Ddiag 0 1 0 2 11 0 34 111 0 456 1111 0 5678 11111 0 678910 L 11111 0 23456 1111 0 4567 111 0 678 11 0 89 1 0 10 0 U 所以由 1 5 式知 MERGEFORMAT 1 18 1 11111 0 23456 33313 0 24527 55555 0 34678 77777 0 45689 99999 0 567810 1111111111 0 678910 BDLU 在用 Jacobi 迭代法求解前 我们先计算迭代矩阵 B 的谱半径 B 的特征值为 8 4 309 0 38 0 932 0 997 1 1 所以 MERGEFORMAT 1 19 1 4 3091 16 i Bmaxi 由 1 16 知迭代矩阵 B 发散 所以无法用 Jacobi 迭代法解 x 的值 3 G S 迭代法求解迭代法求解 由 1 8 式知 1 1 00 50 3330 250 20 167 00 750 250 2250 20 179 00 1040 8680 1350 1310 124 00 0530 0790 910 0920 091 00 0310 0520 0630 9320 07 00 0190 0360 0460 0520 945 LDLU g 1 2 45 1 104 0 626 0 406 0 283 0 207 DLb 通过计算得迭代矩阵的特征值为 1 L 0 0 495113 0 916133 0 99482 0 99985 0 999998 所以 11 0 9999981 16 i Lmaxi 由 1 16 式知迭代矩阵收敛 1 L 令初始值 0 000000 T x 将其代入 1 8 式中 直到 5 1 10 0 1 2 kk max xxk 得 0 9991 0130 9541 0371 0300 966 T k x 9 此时即为的解 k xHxb 4 SOR 迭代法求解迭代法求解 根据 1 18 和 1 19 知 MERGEFORMAT 1 20 4 309 B 再代入 1 13 式 得 2 2 0 1080 452 11 i B 因为为一个虚数 从而说明其最佳松弛因子不存在 故取 1 5 所以 1 1 0 5 0 750 50 3750 30 25 1 1251 18800 0560 0750 08 0 859 0 3520 750 2070 1810 162 0 4540 0900 9650 0510 058 0 290 1210 0960 090 9140 083 0 1420 0230 0130 030 0390 LDLDU 956 其特征值为 0 053183 0 330636 0 896308 0 991614 0 999777 0 999995 故的谱半径为L 0 9999951L 由 1 12 知 SOR 迭代法收敛 令初始值 0 000000 T x 将其代入 1 10 式中 直到 5 1 10 0 1 2 kk max xxk 得 0 9991 0080 973 1 023 1 0120 984 T k x 此时即为的解 k xHxb 二 二 现在逐步增大问题的维数 因为由 一 可知 四种方法只能有 Gauss 消去法 G S 迭代法和 SOR 迭代法求解 故下面只列出这三种求解方法 Hxb 1 当 当 n 20 时 时 Gauss 消去法消去法 1 1 通过计算知道 此时上三角矩阵 U 的第十八行全部变成 0 了 从而导致结果无法计算 故此方法无法求 解 Hxb 因为 n 20 时 Gauss 消去法求解无法算出结果 故以下计算不用此方法了 10 G S 迭代法迭代法 2 2 01020 0 95 1 1 05 x3i i SOR 迭代法迭代法 3 3 01020 0 9 1 1 1 1 2 xi i 2 当 当 n 50 时时 G S 迭代法 迭代法 1 1 02040 0 9 1 1 1 x3i i SOR 迭代法迭代法 2 2 02040 0 95 1 1 05 xi i 3 当 当 n 100 时时 G S 迭代法 迭代法 1 1 11 050100 0 9 1 1 1 x3i i SOR 迭代法 迭代法 2 2 050100 0 95 1 1 05 xi i 五 编写计算程序 五 编写计算程序 一 一 n6 A0 X0 H A Ai j 1 ij 1 j0 n1 for i0 n1 for A X X Xi 0 1 i0 n1 for X AH A X1X X bA X1 H A 1 1 2 1 3 1 4 1 5 1 6 1 2 1 3 1 4 1 5 1 6 1 7 1 3 1 4 1 5 1 6 1 7 1 8 1 4 1 5 1 6 1 7 1 8 1 9 1 5 1 6 1 7 1 8 1 9 1 10 1 6 1 7 1 8 1 9 1 10 1 11 b 2 450 1 593 1 218 0 996 0 846 0 737 1 Gauss 消去法求解消去法求解 12 Gauss A L1identity n L0 Li i 1 i0 n1 for Lj k Aj k Ak k jk1 n1 for AL A L1L L1 k0 n2 for M0L1 1 M1A M MGauss A LM0 UM1 L 1 0 5 0 333 0 25 0 2 0 167 0 1 1 0 9 0 8 0 714 0 0 1 1 5 1 714 1 786 0 0 0 1 2 2 778 0 0 0 0 1 2 5 0 0 0 0 0 1 U 1 0 0 0 0 0 0 5 0 083 0 0 0 0 0 333 0 083 5 55610 3 0 0 0 0 25 0 075 8 33310 3 3 57110 4 0 0 0 2 0 067 9 52410 3 7 14310 4 2 26810 5 0 0 167 0 06 9 92110 3 9 92110 4 5 66910 5 1 43210 6 yL 1 b y 2 45 0 368 0 033 2 06310 3 7 93710 5 1 43210 6 x1U 1 y x1 1 000 1 000 1 000 1 000 1 000 1 000 2 Jacobi 迭代法求解迭代法求解 13 J G A D0 L0 U0 Di j Ai j Li j 0 Ui j 0 ijif Li j Ai j ij if Ui j Ai j ij if j0 n1 for i0 n1 for M0D M1L M2U M MJ G A D1M0 L1M1 U1M2 D1 1 0 0 0 0 0 0 0 333 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 143 0 0 0 0 0 0 0 111 0 0 0 0 0 0 0 091 L1 0 0 5 0 333 0 25 0 2 0 167 0 0 0 25 0 2 0 167 0 143 0 0 0 0 167 0 143 0 125 0 0 0 0 0 125 0 111 0 0 0 0 0 0 1 0 0 0 0 0 0 U1 0 0 0 0 0 0 0 5 0 0 0 0 0 0 333 0 25 0 0 0 0 0 25 0 2 0 167 0 0 0 0 2 0 167 0 143 0 125 0 0 0 167 0 143 0 125 0 111 0 1 0 B1D1 1 L1U1 B1 0 1 5 1 667 1 75 1 8 1 833 0 5 0 1 25 1 4 1 5 1 571 0 333 0 75 0 1 167 1 286 1 375 0 25 0 6 0 833 0 1 125 1 222 0 2 0 5 0 714 0 875 0 1 1 0 167 0 429 0 625 0 778 0 9 0 1eigenvals B1 1 4 309 0 38 0 932 0 997 1 1 max 1 4 309 3 G S 迭代法求解迭代法求解 部分程序在上面 2 Jacobi 迭代法求解 的里面 14 B2D1L1 1 U1 B2 0 0 0 0 0 0 0 5 0 75 0 104 0 053 0 031 0 019 0 333 0 25 0 868 0 079 0 052 0 036 0 25 0 225 0 135 0 91 0 063 0 046 0 2 0 2 0 131 0 092 0 932 0 052 0 167 0 179 0 124 0 091 0 07 0 945 2eigenvals B2 2 0 0 495113 0 916133 0 99482 0 99985 0 999998 max 2 0 999998 G S D1 L1 U1 BD1L1 1 U1 gD1L1 1 b x 0 b 0 1 A0 k0 kk1 x k B x k 1 g Ax k x k 1 max A 10 5 while M0 x k M1 M2k M MG S D1 L1 U1 x3M0 M1 kM2 x3 0 999 1 013 0 954 1 037 1 03 0 966 9 99510 6 k997 4 SOR 迭代法求解迭代法求解 部分程序在上面 2 Jacobi 迭代法求解 的里面 BD1 1 L1U1 B 0 1 5 1 667 1 75 1 8 1 833 0 5 0 1 25 1 4 1 5 1 571 0 333 0 75 0 1 167 1 286 1 375 0 25 0 6 0 833 0 1 125 1 222 0 2 0 5 0 714 0 875 0 1 1 0 167 0 429 0 625 0 778 0 9 0 15 eigenvals B 4 309 0 38 0 932 0 997 1 1 max 4 309 1 2 11 2 10 1080 452i 1 5 B1D1 L1 1 1 D1 U1 B1 0 5 1 125 0 859 0 454 0 29 0 142 0 75 1 188 0 352 0 09 0 121 0 023 0 5 0 0 75 0 0 096 0 013 0 375 0 056 0 207 0 965 0 09 0 03 0 3 0 075 0 181 0 051 0 914 0 039 0 25 0 08 0 162 0 058 0 083 0 956 1eigenvals B1 1 0 053183 0 330636 0 896308 0 991614 0 999777 0 999995 1max 1 10 999995 G S D1 L1 U1 BD1 L1 1 1 D1 U1 g D1 L1 1 b x 0 b 0 1 A0 k0 kk1 x k B x k 1 g Ax k x k 1 max A 10 5 while M0 x k M1 M2k M MG S D1 L1 U1 xM0 M1 kM2 16 x 0 999 1 008 0 973 1 023 1 012 0 984 1010 6 k5 635103 二 二 Gauss 消去法消去法 1 1 前面部分算法和 一 中的 1 Gauss 消去法 类似 此处就不列出了 U 141516171819 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 2 46 10 4 2 515 10 4 2 555 10 4 2 582 10 4 2 598 10 4 2 605 10 4 3 416 10 5 3 66 10 5 3 872 10 5 4 054 10 5 4 21 10 5 4 342 10 5 3 993 10 6 4 537 10 6 5 05 10 6 5 528 10 6 5 971 10 6 6 377 10 6 3 909 10 7 4 78 10 7 5 665 10 7 6 549 10 7 7 419 10 7 8 267 10 7 3 173 10 8 4 249 10 8 5 439 10 8 6 717 10 8 8 06 10 8 9 447 10 8 2 1 10 9 3 149 10 9 4 43 10 9 5 927 10 9 7 62 10 9 9 486 10 9 1 107 10 10 1 915 10 10 3 026 10 10 4 461 10 10 6 23 10 10 8 329 10 10 4 587 10 12 9 532 10 12 1 74 10 11 2 885 10 11 4 445 10 11 6 463 10 11 1 726 10 13 4 579 10 13 9 998 10 13 1 906 10 12 3 288 10 12 5 254 10 12 4 476 10 14 1 341 10 13 3 116 10 13 6 139 10 13 1 077 10 12 1 735 10 12 6 543 10 15 3 662 10 14 1 183 10 13 2 929 10 13 6 11 10 13 1 13 10 12 0003 198 10 15 1 274 10 14 3 468 10 14 009 038 10 15 4 453 10 14 1 396 10 13 3 411 10 13 000000 0000 2 638 10 15 3 942 10 14 000004 161 10 15 G S 迭代法迭代法 2 2 此处程序基本上 一 中的 3 G S 迭代法求解 基本类似 此处就不累赘了 SOR 迭代法迭代法 3 3 此处程序基本上 一 中的 4 SOR 迭代法求解 基本类似 此处就不累赘了 六 解释计算结果 六 解释计算结果 一 一 Gauss 消去法误差分析 消去法误差分析 虽然保留三位小数的情况下