数值分析实验报告

实验五实验五 解线性方程组的直接方法解线性方程组的直接方法 实验实验 5 1 主元的选取与算法的稳定性 问题提出问题提出 Gauss 消去法是我们在线性代数中已经熟悉的 但由于计算机 的数值运算是在一个有限的浮点数集合上进行的 如何才能确保 Gauss 消去法 作为数值算法的稳定性呢 Gauss 消去法从理论算法到数值算法 其关键是主 元的选择 主元的选择从数学理论上看起来平凡 它却是数值分析中十分典型 的问题 实验内容实验内容 考虑线性方程组 nnn RbRAbAx 编制一个能自动选取主元 又能手动选取主元的求解线性方程组的 Gauss 消去过程 实验要求实验要求 1 取矩阵 则方程有解 14 15 15 7 68 168 168 16 bA T x 1 1 1 取 n 10 计算矩阵的条件数 让程序自动选取主元 结果如何 2 现选择程序中手动选取主元的功能 每步消去过程总选取按模最小或 按模尽可能小的元素作为主元 观察并记录计算结果 若每步消去过程总选取 按模最大的元素作为主元 结果又如何 分析实验的结果 3 取矩阵阶数 n 20 或者更大 重复上述实验过程 观察记录并分析不 同的问题及消去过程中选择不同的主元时计算结果的差异 说明主元素的选取 在消去过程中的作用 4 选取其他你感兴趣的问题或者随机生成矩阵 计算其条件数 重复 上述实验 观察记录并分析实验结果 思考题一 思考题一 Vadermonde 矩阵 设 n i i n n i i n i i n i i n nnn n n n x x x x b xxx xxx xxx xxx A 0 0 2 0 1 0 0 2 2 2 22 1 2 11 0 2 00 1 1 1 1 其中 nkkxk 1 0 1 01 1 对 n 2 5 8 计算 A 的条件数 随 n 增大 矩阵性态如何变化 2 对 n 5 解方程组 Ax b 设 A 的最后一个元素有扰动 10 4 再求解 Ax b 3 计算 2 扰动相对误差与解的相对偏差 分析它们与条件数的关系 4 你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗 日或牛顿插值法的原因吗 相关 MATLAB 函数提示 zeros m n 生成 m 行 n 列的零矩阵 ones m n 生成 m 行 n 列的元素全为 1 的矩阵 eye n 生成 n 阶单位矩阵 rand m n 生成 m 行 n 列 0 1 上均匀分布的随机矩阵 diag x 返回由向量 x 的元素构成的对角矩阵 tril A 提取矩阵 A 的下三角部分生成下三角矩阵 triu A 提取矩阵 A 的上三角部分生成上三角矩阵 rank A 返回矩阵 A 的秩 det A 返回方阵 A 的行列式 inv A 返回可逆方阵 A 的逆矩阵 V D eig A 返回方阵 A 的特征值和特征向量 norm A p 矩阵或向量的 p 范数 cond A p 矩阵的条件数 L U P lu A 选列主元 LU 分解 R chol X 平方根分解 Hi hilb n 生成 n 阶 Hilbert 矩阵 5 15 1 实验过程 实验过程 5 1 1 程序 程序 function x gauss n r n input 请输入矩阵 A 的阶数 n A diag 6 ones 1 n diag ones 1 n 1 1 diag 8 ones 1 n 1 1 b A ones n 1 for i 1 4 p input 条件数对应的范数是 p 范数 p pp cond A p end pause m n size A nb n 1 Ab A b r input 请输入是否为手动 手动输入 1 自动输入 0 r for i 1 n 1 if r 0 pivot p max abs Ab i n i ip p i 1 if ip i Ab i ip Ab ip i disp Ab pause end end if r 1 i i ip input 输入 i 列所选元素所处的行数 ip Ab i ip Ab ip i disp Ab pause end pivot Ab i i for k i 1 n Ab k i nb Ab k i nb Ab k i pivot Ab i i nb end disp Ab pause end x zeros n 1 x n Ab n nb Ab n n for i n 1 1 1 x i Ab i nb Ab i i 1 n x i 1 n Ab i i end 5 1 2 实验结果如下 实验结果如下 1 按照实验要求一 取矩阵 A 的阶数 n 10 且自动选取主元 程序结果运行 如下 2 现选择程序中手动选取主元的功能 观察并记录计算结果 选取绝对值最大的元素为主元 程序运行开始如第一问的截图也是求范数故这里不在给出 The answer is 1 1 1 1 1 1 1 1 1 1 选取绝对值最小的元素为主元 The answer is 1 0e 003 INF 0 007 0 0057 0 0410 0 0303 0 3430 0 2577 2 7290 2 0463 2 7308 取矩阵 A 的阶数 n 20 手动选取主元 选取绝对值最大的元素为主元 The answer is 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 选取绝对值最小的元素为主元 The answer is 1 0e 007 Inf 0 0000 0 0000 0 0000 0 0000 0 0000 0 0000 0 0003 0 0002 0 0022 0 0016 0 0175 0 0131 0 1398 0 1049 1 1185 0 8389 8 9478 6 7109 8 9478 修改程序如下 function x gaussong n r n input 请输入矩阵A的阶数 n A hilb n b A ones n 1 for i 1 4 p input 条件数对应的范数是p 范数 p pp cond A p end pause m n size A nb n 1 Ab A b r input 请输入是否为手动 手动输入1 自动输入0 r for i 1 n 1 if r 0 pivot p max abs Ab i n i ip p i 1 if ip i Ab i ip Ab ip i disp Ab pause end end if r 1 i i ip input 输入i列所选元素所处的行数 ip Ab i ip Ab ip i disp Ab pause end pivot Ab i i for k i 1 n Ab k i nb Ab k i nb Ab k i pivot Ab i i nb end disp Ab pause end x zeros n 1 x n Ab n nb Ab n n for i n 1 1 1 x i Ab i nb Ab i i 1 n x i 1 n Ab i i end 所求范数为 自动输入结果为 ans 1 0000 1 00001 0000 1 0000 1 0002 0 9996 1 0007 0 9993 1 0004 0 9999 选取绝对值最大的元素为主元结果为 The answer is NaN NaN NaN NaN NaN Inf Inf Inf 281 3945 283 3708 选取绝对值最小的元素为主元结果为 The answer is NaN NaN NaN Inf 5 8976 1 9243 2 0291 4 9972 23 4548 11 1012 5 1 35 1 3 对实验结果进行分析 对实验结果进行分析 5 1 3 1 对实验要求一的结果进行分析 对于 Gauss 消去法就是用行的初 等变换将原线性方程组系数矩阵转化为简单形式 从而进行求解 缺点是迭代 次数可能较多 效率不高 且在消去过程中不可以将主元素很小的做除数 否 则将导致其他元素数量级的严重增长和舍入误差的扩散 使得计算解不可靠 5 1 3 2 对实验要求二的结果进行分析 通过每次选取最大或最小的主元 可以发现取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结 果比较准确 即选取绝对值小的主元时结果产生了较大的误差 条件数越大产 生的误差就越大 所以应尽量避免很小的数作为除数 5 1 3 3 对实验要求三的结果进行分析 此要求是对要求一和要求二的一 个延续 通过实验结果可以看出若采用很小的数作为主元迭代次数越多导致的 结果越不可靠 甚至出现错误 5 1 3 4 对实验要求四的结果进行分析 对新矩阵进行实验发现依然符合 上述规律 可以知道 在进行迭代时主元的选择与算法的稳定性有密切的联系 选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差 较小 条件数越大对用 gauss 消去法解线性方程组时 对结果产生的误差就越 大 5 1 45 1 4 实验总结 实验总结 1 在用 gauss 消去法解线性方程组时 主元的选取与算法的稳定性有密切的联 系 选取适当的主元有利于得出稳定的算法 2 在算法的过程中 选取绝对值较大的主元比选取绝对值较小的主元更有利于 算法的稳定 选取绝对值最大的元素作为主元时 得出的结果相对较准确较稳 定 3 条件数越小 对用这种方法得出的结果更准确 4 在算除法的过程中要尽量避免使用较小的数做为除数 以免发生结果数量级 加大 使大数吃掉小数 产生舍入误差 实验实验 5 2 线性代数方程组的性态与条件数的估计 问题提出 理论上 线性代数方程组的摄动满足 bAx b b A A AA Acond x x 1 1 矩阵的条件数确实是对矩阵病态性的刻画 但在实际应用中直接计算它显然不 现实 因为计算通常要比求解方程还困难 1 AbAx 实验内容 Matlab 中提供有函数 condest 可以用来估计矩阵的条件数 它给 出的是按 1 范数的条件数 首先构造非奇异矩阵 A 和右端 使得方程是可以精 确求解的 再人为地引进系数矩阵和右端的摄动 使得充分bA 和bA 和 小 实验要求 1 假设方程 Ax b 的解为 x 求解方程 以 1 范数 给出bbxAA 的计算结果 x xx x x 2 选择一系列维数递增的矩阵 可以是随机生成的 比较函数 condest 所需机器时间的差别 考虑若干逆是已知的矩阵 借助函数 eig 很容易给出 cond2 A 的数值 将它与函数 cond A 2 所得到的结果进行比较 3 利用 condest 给出矩阵 A 条件数的估计 针对 1 中的结果给出 的理论估计 并将它与 1 给出的计算结果进行比较 分析所得结果 注 x x 意 如果给出了 cond A 和的估计 马上就可以给出的估计 A 1 A 4 估计著名的 Hilbert 矩阵的条件数 nji ji hhH jinnji 2 1 1 1 5 25 2 实验过程如下 实验过程如下 5 2 1 15 2 1 1 实验要求一的程序如下 实验要求一的程序如下 function n jisuan n a fix 100 rand n 1 x ones n 1 b a x data rand n 0 00001 datb rand n 1 0 00001 A a data B b datb x0 get A B x1 norm x0 x 1 norm x 1 function x get A B m n size A nb n 1 AB A B for i 1 n 1 pivot AB i i for k i 1 n AB k i nb AB k i nb AB k i pivot AB i i nb end end x zeros n 1 x n AB n nb AB n n for i n 1 1 1 x i AB i nb AB i i 1 n x i 1 n AB i i End 5 2 1 25 2 1 2 实验要求一程序运行结果如下 实验要求一程序运行结果如下 系数矩阵 a 为 701814295253 63982473465 80289074421 962620992338 538830595879 897467437769 加扰动后的系数矩阵 A 为 70 18 14 29 52 53 63 98 2 47 34 65 80 28 90 7 44 21 96 26 20 99 23 38 53 88 30 59 58 79 89 74 67 43 77 69 b 值为 236 309 270 302 367 419 加扰动后的 b 值为 236 309 270 302 367 419 data 的值为 4 E 064 E 069 E 063 E 064 E 061 E 06 5 E 063 E 067 E 069 E 069 E 062 E 06 7 E 068 E 064 E 067 E 066 E 066 E 06 5 E 071 E 074 E 064 E 062 E 066 E 06 6 E 067 E 062 E 067 E 068 E 063 E 06 5 E 079 E 066 E 062 E 066 E 065 E 06 datb 的值为 4 E 06 4 E 07 2 E 07 3 E 06 1 E 07 3 E 06 xx 的值为 0 1 0 1 1 0 x0 的值为 1 E 06 x1 的值为 6 8990e 007 5 2 1 35 2 1 3 实验结果为 实验结果为 的计算结果为 6 8990e 007 x xx x x 5 2 2 15 2 2 1 实验要求二的程序如下 实验要求二的程序如下 function cond2 A B A A V1 D1 eig B V2 D2 eig B 1 cond2A sqrt max max D1 sqrt max max D2 end for n 10 10 100 n n A fix 100 randn n condestA condest A cond2 A condA2 cond A 2 pause end 5 2 2 25 2 2 2 实验要求二的程序运行结果如下 实验要求二的程序运行结果如下 NcondestAcond2AcondA2 101 46E 0242 42 204 59E 021 23E 021 23E 02 304 05E 0279 79 402 21E 034 26E 024 26E 02 503 02E 034 08E 024 08E 02 604 75E 037 78E 027 78E 02 704 69E 035 14E 025 14E 02 805 47E 034 89E 024 89E 02 905 66E 035 50E 025 50E 02 1004 47E 035 06E 025 06E 02 5 2 3 15 2 3 1 实验要求三的程序如下 实验要求三的程序如下 function bijiao n a fix 100 rand n 1 x ones n 1 b a x data rand n 0 00001 datb rand n 1 0 00001 A a data B b datb xx geshow A B x1 norm xx x 1 norm x 1 x2 cond A 1 norm inv A norm xx x norm xx x norm A norm datb norm B datx abs x1 x2 5 2 3 25 2 3 2 实验要求三的程序运行结果如下 实验要求三的程序运行结果如下 给出对的估计是 7 8125e 007 x xx x x 的理论结果是 3 7297e 007 x xx x x 结果相差 3 0828e 007 5 2 4 15 2 4 1 实验要求四的程序如下 实验要求四的程序如下 for n 4 11 n n Hi hilb n cond1Hi cond Hi 1 cond2Hi cond Hi 2 condinfHi cond Hi inf pause end 5 2 4 25 2 4 2 实验要求四的程序运行结果如下 实验要求四的程序运行结果如下 ncond1Hicond2HicondinfHi 42 84E 041 55E 042 84E 04 59 44E 054 77E 059 44E 05 62 91E 071 50E 072 91E 07 79 85E 084 75E 089 85E 08 83 39E 101 53E 103 39E 10 91 10E 124 93E 111 10E 12 103 54E 131 60E 133 54E 13 111 23E 155 22E 141 23E 15 讨论 讨论 线性代数方程组的性态与条件数有着很重要的关系 既矩阵的条件数是刻 画矩阵性质的一个重要的依据 条件数越大 矩阵 病态 性越严重 在解线 性代数方程组的过程中较容易产生比较大的误差 则在实际问题的操作过程中 我们必须要减少对条件数来求解 把条件数较大的矩阵化成条件数较小的矩阵 来进行求解 实验总结 在本次实验中 使我们知道了矩阵条件数对线性代数方程组求解的影响 条件数越大 对最后解的影响的越大 hilbert 矩阵是一个很 病态 的矩阵 他 的条件数随着阶数的增加而增大 每增加一阶 条件数就增大一个数量级 在 求解的过程中要尽量避免 hilbert 矩阵 实验六实验六 解线性方程组的迭代法解线性方程组的迭代法 实验实验 6 1 病态的线性方程组的求解 问题提出问题提出 理论的分析表明 求解病态的线性方程组是困难的 实际情况 是否如此 会出现怎样的现象呢 实验内容实验内容 考虑方程组 Hx b 的求解 其中系数矩阵 H 为 Hilbert 矩阵 nji ji hhH jinnji 2 1 1 1 这是一个著名的病态问题 通过首先给定解 例如取为各个分量均为 1 再计算出右端 b 的办法给出确定的问题 实验要求实验要求 1 选择问题的维数为 6 分别用 Gauss 消去法 J 迭代法 GS 迭代法 和 SOR 迭代法求解方程组 其各自的结果如何 将计算结果与问题的解比较 结论如何 2 逐步增大问题的维数 仍然用上述的方法来解它们 计算的结果如 何 计算的结果说明了什么 3 讨论病态问题求解的算法 由于本实验用到迭代法 故先给出迭代法相关资料 由于本实验用到迭代法 故先给出迭代法相关资料 迭代法也称辗转法 是一种不断用变量的旧值递推新值的过程 跟迭代法相对应的是直 接法 或者称为一次解法 即一次性解决问题 迭代法又分为精确迭代和近似迭代 二分 法 和 牛顿迭代法 属于近似迭代法 迭代算法是用计算机解决问题的一种基本方法 它利 用计算机运算速度快 适合做重复性操作的特点 让计算机对一组指令 或一定步骤 进 行重复执行 在每次执行这组指令 或这些步骤 时 都从变量的原值推出它的一个新值 迭代是数值分析 中通过从一个初始估计出发寻找一系列近似解来解决问题 一般是 解方程或者方程组 的过程 为实现这一过程所使用的方法统称为迭迭代代法法 Iterative Method 一般可以做如下定义 对于给定的线性方程组 这里的同为fBxx fBx 矩阵 任意线性方程组都可以变换成此形式 用公式 括号中为上标 fkBxkx 1 代表迭代 k 次得到的 x 初始时 k 0 逐步带入求近似解的方法称为迭代法 或称一阶定常 迭代法 如果 k 趋向无穷大时存在 记为 称此迭代法收敛 显然就是此 limkx x x 方程组的解 否则称为迭代法发散 跟迭代法相对应的是直接法 或者称为一次解法 即 一次性的快速解决问题 例如通过开方解决方程 一般如果可能 直接解法总是43 x 优先考虑的 但当遇到复杂问题时 特别是在未知量很多 方程为非线性时 我们无法找 到直接解法 例如五次以及更高次的代数方程没有解析解 参见阿贝尔定理 这时候或许 可以通过迭代法寻求方程 组 的近似解 最常见的迭代法是牛顿法 其他还包括最速下 降法 共轭迭代法 变尺度迭代法 最小二乘法 线性规划 非线性规划 单纯型法 惩 罚函数法 斜率投影法 遗传算法 模拟退火法等等 利用迭代算法解决问题 需要做好 以下三个方面的工作 1 确定迭代变量 在可以用迭代算法解决的问题中 至少存在一个直接或间接地不断由旧值递推出新值 的变量 这个变量就是迭代变量 2 建立迭代关系式 所谓迭代关系式 指如何从变量的前一个值推出其下一个值的公式 或关系 迭代关 系式的建立是解决迭代问题的关键 通常可以顺推或倒推的方法来完成 3 对迭代过程进行控制 在什么时候结束迭代过程 这是编写迭代程序必须考虑的问题 不能让迭代过程无休 止地重复执行下去 迭代过程的控制通常可分为两种情况 一种是所需的迭代次数是个确 定的值 可以计算出来 另一种是所需的迭代次数无法确定 对于前一种情况 可以构建 一个固定次数的循环来实现对迭代过程的控制 对于后一种情况 需要进一步分析出用来 结束迭代过程的条件 实验过程 实验过程 6 1 1 1 实验要求一的实验要求一的 gauss 实现程序如下 实现程序如下 Gauss 消去法程序 function x gaussong n r A hilb n b A ones n 1 for i 1 4 p input 条件数对应的范数是 p 范数 p pp cond A p end pause m n size A nb n 1 Ab A b r input 请输入是否为手动 手动输入 1 自动输入 0 r for i 1 n 1 if r 0 pivot p max abs Ab i n i ip p i 1 if ip i Ab i ip Ab ip i disp Ab pause end end if r 1 i i ip input 输入 i 列所选元素所处的行数 ip Ab i ip Ab ip i disp Ab pause end pivot Ab i i for k i 1 n Ab k i nb Ab k i nb Ab k i pivot Ab i i nb end disp Ab pause end x zeros n 1 x n Ab n nb Ab n n for i n 1 1 1 x i Ab i nb Ab i i 1 n x i 1 n Ab i i end 6 1 1 2 实验要求一的实验要求一的 gauss 实现程序运行结果如下 实现程序运行结果如下 6 1 2 1 实验要求一的实验要求一的 J 程序如下 程序如下 function n JD n A hilb n for i 1 n a0 i 1 x i 0 end b A a0 for i 1 100 y x for j 1 n x j b j A j j for k 1 j 1 x j x j A j k y k A j j end for k j 1 n x j x j A j k y k A j j end end end x 6 1 2 2 实验要求一的实验要求一的 J 程序运行结果如下 程序运行结果如下 6 1 3 1 实验要求一的实验要求一的 GS 程序如下 程序如下 function n GS n A2 hilb n for i 1 n a02 i 1 x2 i 0 end b2 A2 a02 for i 1 for j 1 n x2 j b2 j A2 j j for k 1 j 1 x2 j x2 j A2 j k x2 k A2 j j end for k j 1 n x2 j x2 j A2 j k x2 k A2 j j end end end x2 6 1 3 2 实验要求一的实验要求一的 GS 程序运行结果如下 程序运行结果如下 6 1 4 1 实验要求一的实验要求一的 SOR 程序如下 程序如下 function n ss SOR n ss A3 hilb n for i 1 n a03 i 1 x3 i 0 end b3 A3 a03 for i 1 for j 1 n rc x3 j x3 j b3 j A3 j j for k 1 j 1 x3 j x3 j A3 j k x3 k A3 j j end for k j 1 n x3 j x3 j A3 j k x3 k A3 j j end x3 j 1 ss rc ss x3 j end end x3 6 1 4 2 实验要求一的实验要求一的 SOR 程序运行结果如下 程序运行结果如下 注意由于运行结果迭代过程很长故不体现在实验报告中 6 2 1 实验要求二的实验要求二的 gauss 程序运行结果如下 程序运行结果如下 x 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9975 6 2 2 实验要求二的实验要求二的 J 程序运行结果如下 程序运行结果如下 6 2 3 实验要求二的实验要求二的 GS 程序运行结果如下 程序运行结果如下 6 2 4 实验要求二的实验要求二的 SOR 程序运行结果如下 程序运行结果如下 实验要求二的结果分析 实验要求二的结果分析 选择问题的维数为 20 时 1 用 Gauss 消去法求得的解与精确解相差很大 说明能否得到优秀的解取决于算法的稳定 性 如果算法不够稳定产生的结果将变的无法理喻 2 取初始向量为 0 用 J 迭代方法迭代发散 无法求解 3 取初始向量为 0 用 GS 迭代方法迭代不发散 能求得解 但收敛非常缓慢 而且迭代 次数越多 与准确解的偏差就越大 说明 GS 迭代适合迭代次数少的 但是通常我们无法 得知需要迭代的次数 4 取初始向量为 0 用 SOR 迭代方法迭代不发散 能求得解 但同样收敛非常缓慢 总之 从上面的结果可以看出当病态问题的阶数升高时作为直接法的 Gauss 消去法又能求 解变成不能求解 而 GS 和 SOR 迭代法在阶数升高时仍能求解 但在阶数较低时直接法能 求得精确解而迭代发却总存在一定的误差 可见直接法与迭代法各有各的优势与不足 关于病态问题的求解 主要的方法是对原方程作某些预处理 以降低系数矩阵的条件数 可以采取将系数矩阵 A 的每一行本别乘上适当常数的方法 即找到可逆的对角阵和 1 D 使方程组化为 理论上最好选择对角阵满足 2 DAxb 1 1212 D AD yDb xD y 12 D D 1212 min cond D ADcond D AD 补充 补充 本实验用到了 Gauss 消去法 J 迭代法 GS 迭代法 SOR 迭代 故对其重 新理解学习了一下 Gauss 消去法 消去法 1 将其增广矩阵化为行阶梯形 2 若最后有形如 0 0 0 a a0 的行则无解 3 若含有自由变量则有无穷组解 4 原方程有唯一解 采用回代求解 至于有无穷组解的方程组的求解 需将其化为行最简形矩阵 其方法称为高斯 若尔当 消元法 J 迭代法 迭代法 一 算法理论 迭代格式的引出是依据迭代法的基本思想 构造一个向量系列 Jacobi n x 使其收敛至某个极限 则就是要求的方程组的准确解 X X 迭代Jacobi 将方程组 nnnnnn nn nn bxaxaxa bxaxaxa bxaxaxa 2211 22222121 11212111 1 在假设 改写成0 ii a nnnnnnn nn nn gxbxbxbx gxbxbxbx gxbxbxbx 112211 223231212 113132121 2 如果引用系数矩阵 及向量 nnn n aa aa A 1 111 0 0 1 1 n n b b B nnn g g g b b b x x X 111 方程组 1 和 2 分别可写为 及 这样就得到了bAX gBXX 迭代格式用迭代解方程组时 就可任意取Jacobi 0 1 gBXX kk JacobibAX 初值带入迭代可知式 然后求 但是 比较大的时 0 XgBXX kk 1k k X limn 候 写方程组 1 和 2 是很麻烦的 如果直接由 能直接得到 AbB 就是矩阵与向量的运算了 那么如何得到 呢gBg 实际上 如果引进非奇异对角矩阵0 ii a nn a a D 0 0 0 11 将分解成 要求的解 实质上就有 ADDAA bAX DXXDAAX 而是非奇异的 所以存在 从而有D 1 DXADAXDX 我们在这里不妨令 就得到迭代bDAXDX 11 ADIB 1 bDg 1 Jacobi 格式 现在考虑迭代法的计算程序gBXX kk 1 Jacobi float a 3 3 10 2 1 2 10 1 1 2 5 float b 3 3 15 10 分别代表的X 系数和等号右边的常数项 即 1052 15102 3210 333231 232211 131211 xxx xxx xxx 先输入方程 运行 main 函数 如果 first 不为 null 则执行 if 括号里的 否则执 行 else 里面的 最后会调用方法 sum 在 sum 中 sum a m n x n sum y b m sum a m m 之后运行 for 循环 最后输出结果 算法结束 二 算法框图 开始 读入 ba sumnxnmasum GS 迭代法 迭代法 main isumiy 输出 结束 调用 方法 sum 1 高斯 塞德尔迭代法 公式的矩阵形式 首先将高斯 塞德尔迭代法的公式表示为矩阵形式 为此设 ULDA 这里 是系数矩阵的对角部分 是严格下三角部分 是严格上三角部分 则高DAL U 斯 塞德尔迭代法的公式可表示为 1 0 1 1 1 mbUXLXDX mmm 用矩阵 乘等式两边得D 1 0 1 mbUXXLD mm 再用矩阵 乘等式两边得 1 LD 1 0 1 1 mbLDUXLDX mm 其中矩阵 称为高斯 塞德尔迭代矩阵 ULDB 1 1 由此可见 高斯 塞德尔迭代法是一般迭代法中迭代矩阵为 的特殊情ULDB 1 1 形 需要指出的是 由于矩阵 难于计算 所以式 2 多用在理论分析中 ULD 1 2 高斯 塞德尔迭代法计算框图 见图 高斯 塞德尔迭代法计算框图 松弛法 其中 ii k i k i n ij k jij i j iiji ii k i a r x xakxab a x 1 1 1 1 1 1 1 ij k jij ij k jiji k i xaxabr 1 1 相当于在的基础上加个余项生成 下面令 希望 k i x 1 k i x ii k i k i k i a r xx 1 1 通过选取来加速收敛 这就是松弛法 实验七实验七 非线性方程求根非线性方程求根 实验 7 1 迭代法 初始值与收敛性 实验目的 初步认识非线性问题的迭代法与线性问题迭代法的差别 探讨迭代 法及初始值与迭代收敛性的关系 问题提出 迭代法是求解非线性方程的基本思想方法 与线性方程的情况一样 其构造方法可以有多种多样 但关键是怎样才能使迭代收敛且有较快的收敛速 度 实验内容 考虑一个简单的代数方程 01 2 xx 针对上述方程 可以构造多种迭代法 如 1 7 1 2 1 nn xx 2 7 1 1 1 n n x x 3 7 1 1 nn xx 在实轴上取初始值 x0 请分别用迭代 7 1 7 3 作实验 记录各算法的迭 代过程 实验要求 1 取定某个初始值 分别计算 7 1 7 3 迭代结果 它们的收敛性如何 重复选取不同的初始值 反复实验 请自选设计一种比较形象的记录方式 如 利用 Matlab 的图形功能 分析三种迭代法的收敛性与初值选取的关系 2 对三个迭代法中的某个 取不同的初始值进行迭代 结果如何 试分析迭 代法对不同的初值是否有差异 3 线性方程组迭代法的收敛性是不依赖初始值选取的 比较线性与非线性问 题迭代的差异 有何结论和问题 实验过程 实验过程 7 1 程序 程序 7 1 1 对于第一个迭代方程的程序 对于第一个迭代方程的程序 保存为 diedai71 clear clc a 1 5 b 2 5 y00 0 x00 input 请输入第一个函数的初值 x00 x linspace a b 80 y0 x 计算直线y x y1 diedai7f1 x 计算迭代函数y f x clear y y y0 y1 plot x y linewidth 1 legend y x y f1 title x n 1 x n 2 1 输出标题 hold on plot a b 0 0 k 0 0 a b k axis a b a b 画坐标轴 z for i 1 15 画蛛网图 迭代过程为n 15次 xt 1 x00 yt 1 y00 决定始点坐标 xt 2 diedai7f1 xt 1 决定终点坐标 yt 2 diedai7f1 xt 1 diedaiplot71 xt yt 0 6 画蛛网图 if i 5 pause 按任意键逐次观察前5次迭代的蛛网图 end x00 xt 2 y00 yt 2 将本次迭代的终点作为下次的始点 z z xt 1 保存迭代点 end 保存为 diedai7f1 function y diedai7f1 x y x x 1 保存为 diedaiplot72 function out diedaiplot72 x y p 画一次迭代的蛛网图 改变p调节箭头的大小 u 1 0 v 1 y 2 y 1 画出始点 x 1 y 1 终点 x 2 y 2 的有向折 线段 u 2 eps v 2 eps X x 1 x 1 Y y 1 y 2 h quiver X Y u v p set h color red hold on u 1 x 2 x 1 v 1 0 u 2 eps v 2 eps h quiver x 1 x 2 y 2 y 2 u v p set h color red plot x 1 x 1 x 2 y 1 y 2 y 2 r 7 1 1 对于第一个迭代方程的程序运行结果 对于第一个迭代方程的程序运行结果 注 由于可以估计出方程的根 故选取根附近的值开始迭代 注 由于可以估计出方程的根 故选取根附近的值开始迭代 当 x 1 5 时程序运行如下图 当 x 1 时程序运行如下图 当 x 0 5 时程序运行如下图 当 x 0 5 时程序运行如下图 当 x 0 时程序运行如下图 7 2 1 对于第二个迭代方程的程序 对于第二个迭代方程的程序 保存为 diedai72 clear clc a 0 1 b 6 5 y00 0 x00 input 请输入第二个函数的初值 x00 x linspace a b 80 y0 x 计算直线y x y1 diedai7f2 x 计算迭代函数y f x clear y y y0 y1 plot x y linewidth 1 legend y x y f1 title x n 1 x n 2 1 输出标题 hold on plot a b 0 0 k 0 0 a b k axis a b a b 画坐标轴 z for i 1 15 画蛛网图 迭代过程为n 15次 xt 1 x00 yt 1 y00 决定始点坐标 xt 2 diedai7f2 xt 1 决定终点坐标 yt 2 diedai7f2 xt 1 diedaiplot72 xt yt 0 6 画蛛网图 if i 5 pause 按任意键逐次观察前5次迭代的蛛网图 end x00 xt 2 y00 yt 2 将本次迭代的终点作为下次的始点 z z xt 1 保存迭代点 end 保存为迭代 7f2 function y diedai7f2 x y 1 1 x 保存为 diedaiplot72 function out diedaiplot72 x y p 画一次迭代的蛛网图 改变p调节箭头的大小 u 1 0 v 1 y 2 y 1 画出始点 x 1 y 1 终点 x 2 y 2 的有向折线段 u 2 eps v 2 eps X x 1 x 1 Y y 1 y 2 h quiver X Y u v p set h color red hold on u 1 x 2 x 1 v 1 0 u 2 eps v 2 eps h quiver x 1 x 2 y 2 y 2 u v p set h color red plot x 1 x 1 x 2 y 1 y 2 y 2 r 7 2 2 对于第二个迭代方程的程序运行结果 对于第二个迭代方程的程序运行结果 注 由于可以估计出方程的根 故选取根附近的值开始迭代 注 由于可以估计出方程的根 故选取根附近的值开始迭代 当 x 1 5 时程序运行如下图 当 x 5 时程序运行如下图 当 x 2 5 时程序运行如下图 当 x 0 5 时程序运行如下图 7 3 1 对于第三个迭代方程的程序 对于第三个迭代方程的程序 保存为diedai73 clear clc a 0 b 2 y00 0 x00 input 请输入第三个函数的初值 x00 x linspace a b 80 y0 x 计算直线y x y1 diedai7f3 x 计算迭代函数y f x clear y y y0 y1 plot x y linewidth 1 legend y x y f1 title x n 1 x n 2 1 输出标题 hold on plot