计算方法——共轭梯度法求解线性方程组.pdf

计算方法上机报告 1 计算方法上机报告 1 共轭梯度法求解线性方程组 1.1 算法原理及程序框图 当线性方程组 Ax=b 的系数矩阵 A 是对称正定矩阵是，可以采用共轭梯度法对该方程组进行求解，可以证明，式(1)所示的 n 元二次函数 TT1()2f xx Axb x(1)取得极小值点 x*是方程 Ax=b 的解。共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小化的问题。从任意给定的初始点出发，沿一组关于矩阵 A 的共轭方向进行线性搜索，在无舍入误差的假定下，最多迭代 n 次（其中 n 为矩阵 A 的阶数），就可求得二次函数的极小点，也就求得线性方程组 Ax=b 的解。其迭代格式为公式(2)。(1)()()kkkkxxd(2)共轭梯度法中关键的两点是迭代格式(2)中最佳步长k和搜索方向 d(k)的确定。其中k可以通过一元函数f(x(k)+d(k)的极小化来求得，其表达式为公式(3)；取 d(0)=r(0)=b-Ax(0)，则 d(k+1)=r(k+1)+kd(k)，要求 d(k+1)满足(d(k+1),Ad(k)=0，可得k的表达式(4)。TTkkkkkrddAd(3)(1)T()()T()kkkkkrddd AA(4)经过一系列的证明和简化，最终可得共轭梯度法的计算过程如下，计算程序框图如图 1。(1)给定初始计算向量 x(0)即精度0；(2)计算 r(0)b-Ax(0)，取 d(0)=r(0)；(3)for k=0 to n-1 do TTkkkkkrrdAd；1kkkkxxd；1 共轭梯度法求解线性方程组 2 11kkrbAx；若1|kr或1kn，则输出近似解 x(k+1)，停止；否则，转；12222|kkkrr；11kkkkdrd；end do 图 1 共轭梯度法求解线性方程组程序框图 1.2 程序使用说明 共轭梯度法求解线性方程组的 matlab 程序见附录 1，该程序可以求解系数矩阵为对称正定矩阵的线性方程组。在使用该程序时，可将程序复制到 matlab 命令窗口中直接运行或者复制到编辑窗口中保存运行，运行时刻根据提示输入，直至得到结果。开始运行程序时，会出现提示“请选择系数矩阵、右端项以及系数矩阵阶数的输计算方法上机报告 3 入方式：从文件中输入数据输入 1，从命令窗口输入数据请输入 2。”此时需要选择数据输入的方式（方式 1 和方式 2），即文件输入（选择 1）或者手动输入（选择 2）。当线性方程组 Ax=b 的系数矩阵的阶数较大时，将该系数矩阵 A、右端项 b 以及系数矩阵的阶数n保存为txt格式放在D盘的根目录下并分别命名为data_A、data_b和data_n，并输入 1 选择文件输入。若系数矩阵 A 的阶数较小，使用手动输入工作量小时，可在命令框中输入 2 选择手动输入。选择手动输入时需要输入系数矩阵 A、右端项 b 以及系数矩阵的阶数 n 这三个量，其输入格式与 matlab 中矩阵、列向量和数的输入格式要求相同。A=aijn n的输入格式为a11,a12,a1n;a21,a22,a2n;an1,an2,ann，b=(bi)T的输入格式为b1;b2;bn，n 直接输入对应的数值即可。定义完需要求解的线性方程组之后，接下来会出现提示“输入计算要求的精度 epsilon：”，按照提示输入精度值，如要求的精度为 10-6时输入 1e-6 即可。以上数据的输入全部完成，每次按提示输入完成时按 Enter 键继续下一步。在程序运行过程中，屏幕上会显示残差|r|2随着迭代步数的变化散点图，程序运行完成之后，matlab 命令窗口会显示线性方程的解 x，同时该解会被保存到 D 盘根目录下名为 data_x_result 的 Excel 文件中，方便之后的数据处理（注意注意在求解一个新方程组之前查看 D 盘根目录，将上次得到的文件删除，以免影响本次计算的结果）。1.3 算例计算结果(1)数值分析课本第 29 页的例题 2.2.2，下面采用共轭梯度法来求解线性方程组 Ax=b，其中 918927184504590126927459135A，12168b 由于该方程组的系数矩阵以及右端项都比较简单，因此采用从 matlab 命令窗口手动输入的方式来输入数据，取计算精度为 10-6，运行过程及结果如图 2 和图 3（由于迭代的初始值为随机产生，因此每次得到的残量图会有所不同，但最终都趋于 0）：1 共轭梯度法求解线性方程组 4 图 2 命令窗口显示的运行结果 图 3 残量随迭代步数的变化(2)数值分析课本第 113 页的计算实习题 3.2，用共轭梯度法求解线性方程组Ax=b，其中 2112112112A，1001b 11.522.530510152025迭代步数残量计算方法上机报告 5 矩阵 A 的阶数取 200 进行求解。由于该线性方程组的系数矩阵阶数比较大，而且具有一定的规律，因此首先用matlab 编程将系数矩阵、右端项以及阶数保存在 D 盘根目录的三个文件中（生成系数矩阵，右端项以及阶数的程序见附录 2），然后运行共轭梯度法程序进行方程组的求解。最终的运行结果为图 4 和图 5。程序运行之后 D 盘根目录下生成的文件如图 6 所示。图 4 命令窗口显示的运行结果 图 5 残量随迭代步数的变化 图 6 程序运行后 D 盘根目录下生成的文件02040608010000.20.40.60.811.21.41.6迭代步数残量