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

1、计算方法上机报告计算方法上机报告1 共轭梯度法求解线性方程组1.1 算法原理及程序框图当线性方程组 Ax = b 的系数矩阵 A 是对称正定矩阵是，可以采用共轭梯度法对该方程组进行求解，可以证明，式(1)所示的 n 元二次函数f ( x) = 1 xT Ax - bT x(1)2取得极小值点 x*是方程 Ax = b 的解。共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小化的问题。从任意给定的初始点出发，沿一组关于矩阵 A 的共轭方向进行线性搜索，在无舍入误差的假定下，最多迭代 n 次（其中 n 为矩阵 A 的阶数），就可求得二次函数的极小点，也就求得线性方程组 Ax =

2、 b 的解。其迭代格式为公式(2)。5kx(k +1) = x(k ) +a d (k )(2)共轭梯度法中关键的两点是迭代格式(2)中最佳步长ak 和搜索方向 d (k)的确定。其中ak 可以通过一元函数f(a) = f(x(k)+ad(k)的极小化来求得，其表达式为公式(3)；取 d (0)= r(0) = b-Ax(0)，则 d(k+1) = r(k+1) +bkd(k)，要求 d(k+1)满足 (d(k+1) , Ad(k) = 0，可得bk 的表达式(4)。r(k )T d (k )ak = d (k )T Ad (k )(3)b =- r(k +1)T Ad (k )kd (k )

3、T Ad (k )(4)经过一系列的证明和简化，最终可得共轭梯度法的计算过程如下，计算程序框图如图 1。(1) 给定初始计算向量 x(0)即精度e >0； (2) 计算 r(0) = b -Ax(0)，取 d(0) = r(0)；(3) fork =0ton-1dor(k )T r(k )ak = d (k )T Ad (k ) ；kx(k +1) = x(k ) +a d (k ) ；r(k +1) = b + Ax(k +1) ；若| r(k +1) |£ e 或k +1 = n ，则输出近似解 x(k+1)，停止；否则，转；| r(k +1) |2| r(k) |2bk

4、= 2 ；2kd (k +1) = r(k +1) + b d (k ) ；end do图 1 共轭梯度法求解线性方程组程序框图1.2 程序使用说明共轭梯度法求解线性方程组的 matlab 程序见附录 1，该程序可以求解系数矩阵为对称正定矩阵的线性方程组。在使用该程序时，可将程序复制到 matlab 命令窗口中直接运行或者复制到编辑窗口中保存运行，运行时刻根据提示输入，直至得到结果。开始运行程序时，会出现提示“请选择系数矩阵、右端项以及系数矩阵阶数的输入方式：从文件中输入数据输入 1，从命令窗口输入数据请输入 2。”此时需要选择数据输入的方式（方式 1 和方式 2），即文件输入（选择 1）或者

5、手动输入（选择 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=

6、(bi)T 的输入格式为b1;b2;bn，n 直接输入对应的数值即可。定义完需要求解的线性方程组之后，接下来会出现提示“输入计算要求的精度 epsilon：”，按照提示输入精度值，如要求的精度为 10-6 时输入 1e-6 即可。以上数据的输入全部完成，每次按提示输入完成时按 Enter 键继续下一步。在程序运行过程中，屏幕上会显示残差|r|2 随着迭代步数的变化散点图，程序运行完成之后，matlab 命令窗口会显示线性方程的解 x，同时该解会被保存到 D 盘根目录下名为 data_x_result 的 Excel 文件中，方便之后的数据处理（注意在求解一个新方程组之前查看 D 盘根目录，将上

7、次得到的文件删除，以免影响本次计算的结果）。1.3 算例计算结果(1) 数值分析课本第 29 页的例题 2.2.2，下面采用共轭梯度法来求解线性方程组 Ax = b，其中æ 9189-27 öæ 1 öç 18450-45 ÷ç 2 ÷A = ç÷ ， b = ç÷ç 901269 ÷ç16 ÷ç -27-459135 ÷ç 8 ÷èøèø由于该方程组的系数

8、矩阵以及右端项都比较简单，因此采用从 matlab 命令窗口手动输入的方式来输入数据，取计算精度为 10-6，运行过程及结果如图 2 和图 3（由于迭代的初始值为随机产生，因此每次得到的残量图会有所不同，但最终都趋于 0）：图 2 命令窗口显示的运行结果2520残量15105011.522.53迭代步数图 3 残量随迭代步数的变化(2) 数值分析课本第 113 页的计算实习题 3.2，用共轭梯度法求解线性方程组Ax = b，其中æ -21öæ -1öç 1-21÷ç 0 ÷0 ÷÷÷÷ç÷çA = ç÷ ， b = çç1-21 ÷çç÷çç1-2 ÷ç -1÷èøèø矩阵 A 的阶数取 200 进行求解。由于该线性方程组的系数矩阵阶数比较大，而且具有一定的规律，因此首先用matlab 编程将系数矩阵、右端项以及阶数保存在 D 盘根目录的三个文件中（生成系数矩阵，右端项以及阶数的程序见附录 2），然后运行共轭梯度法程序进行方程