4
1 1
1 0 ，b 2 1 0 1 1 2
计算方法上机报告
矩阵 A 的阶数取 200 进行求解。 由于该线性方程组的系数矩阵阶数比较大，而且具有一定的规律，因此首先用 matlab 编程将系数矩阵、右端项以及阶数保存在 D 盘根目录的三个文件中（生成系数 矩阵， 右端项以及阶数的程序见附录 2） ， 然后运行共轭梯度法程序进行方程组的求解。 最终的运行结果为图 4 和图 5。程序运行之后 D 盘根目录下生成的文件如图 6 所示。
18 9 27 9 1 18 45 0 45 2 ，b A 9 16 0 126 9 27 45 9 135 8
由于该方程组的系数矩阵以及右端项都比较简单，因此采用从 matlab 命令窗口手 动输入的方式来输入数据，取计算精度为 10-6，运行过程及结果如图 2 和图 3（由于迭 代的初始值为随机产生，因此每次得到的残量图会有所不同，但最终都趋于 0） ：
k
k
r d k T k d Ad
k T k
(3)
r ( k 1)T Ad ( k ) d ( k )T Ad ( k )
(4)
经过一系列的证明和简化，最终可得共轭梯度法的计算过程如下，计算程序框图 如图 1。 (1) 给定初始计算向量 x(0)即精度>0； (2) 计算 r(0) = b -Ax(0)，取 d(0) = r(0)； (3) for k =0 to ① k ②
|| r ||2 k k 22 ； || r ||2
d
k 1
r
k 1
k d ；
k
end do
图 1 共轭梯度法求解线性方程组程序框图
1.2 程序使用说明 共轭梯度法求解线性方程组的 matlab 程序见附录 1，该程序可以求解系数矩阵为 对称正定矩阵的线性方程组。在使用该程序时，可将程序复制到 matlab 命令窗口中直 接运行或者复制到编辑窗口中保存运行，运行时刻根据提示输入，直至得到结果。 开始运行程序时，会出现提示“请选择系数矩阵、右端项以及系数矩阵阶数的输
计算方法上机报告
计算方法上机报告
1 共轭梯度法求解线性方程组
1.1 算法原理及程序框图 当线性方程组 Ax = b 的系数矩阵 A 是对称正定矩阵是，可以采用共轭梯度法对该 方程组进行求解，可以证明，式(1)所示的 n 元二次函数 1 f ( x ) x T Ax bT x (1) 2 取得极小值点 x*是方程 Ax = b 的解。共轭梯度法是把求解线性方程组的问题转化为求 解一个与之等价的二次函数极小化的问题。从任意给定的初始点出发，沿一组关于矩 阵 A 的共轭方向进行线性搜索，在无舍入误差的假定下，最多迭代 n 次（其中 n 为矩 阵 A 的阶数） ，就可求得二次函数的极小点，也就求得线性方程组 Ax = b 的解。其迭 代格式为公式(2)。
x( k 1) x( k ) k d ( k )
(2)
(k)
共轭梯度法中关键的两点是迭代格式(2)中最佳步长k 和搜索方向 d
(k) (k)
的确定。其
中k 可以通过一元函数f(x +d )的极小化来求得，其表达式为公式(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)。
3
1 共轭梯度法求解线性方程组
图 2 命令窗口显示的运行结果
25
20
15
残量
10 5 0 1
1.5
2
2.5
3
迭代步数
图 3 残量随迭代步数的变化
(2) 《数值分析》课本第 113 页的计算实习题 3.2，用共轭梯度法求解线性方程组 Ax = b，其中
2 1 1 2 A
2
计算方法上机报告
入方式：从文件中输入数据输入 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=[aij]n×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，其中
图 4 命令窗口显示的运行结果
1.6 1.4 1.2 1
残量
0.8 0.6 0.4 0.2 0 0 20 40 60 80 100
迭代步数
图 5 残量随迭代步数的变化
图 6 程序运行后 D 盘根目录下生成的文件
5
k T
n-1 do
k
r r ； k T k d Ad
xkLeabharlann 1 x k d ；
k k
1
1 共轭梯度法求解线性方程组
③
r
k 1
b Ax
k 1
k 1
；
④ 若 || r k 1 || 或 k 1 n ，则输出近似解 x(k+1)，停止；否则，转⑤； ⑤ ⑥