数值代数实验报告
一、实验名称：用共轭梯度法解线性方程组。
二、实验目的：进一步熟悉理解掌握共轭梯度法解法思路，提高matlab 编程能力。 三、实验要求：已知线性方程矩阵，应用共轭梯度法在相关软件编程求解线性方程组的解。 四、实验原理：
1．共轭梯度法：
考虑线性方程组
Ax b =
的求解问题，其中A 是给定的n 阶对称正定矩阵，b 是给定的n 维向量，x 是待求解的n 维向量.为此,定义二次泛函
()2T T x x Ax b x ϕ=-.
定理1 设A 对称正定,求方程组Ax b =的解，等价于求二次泛函()x ϕ的极小值点. 定理1表明，求解线性方程组问题就转化为求二次泛函()x ϕ的极小值点问题.求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量0x ，确定一个下山方向0p ，沿着经过点0x 而方向为0p 的直线00x x p α=+找一个点
1000x x p α=+，
使得对所有实数α有
()()00000x p x p ϕαϕα+≤+，
即在这条直线上1x 使()x ϕ达到极小.然后从1x 出发，再确定一个下山的方向1p ，沿着直
线11x x p α=+再跨出一步，即找到1α使得()x ϕ在2111x x p α=+达到极小：
()()11111x p x p ϕαϕα+≤+.
重复此步骤，得到一串
012,,,αααL 和 012,,,p p p L ，
称k p 为搜索方向，k α为步长.一般情况下，先在k x 点找下山方向k p ，再在直线
k k x x p α=+上确定步长k α使
()(),k k k k k x p x p ϕαϕα+≤+
最后求出1k k k k x x p α+=+.然而对不同的搜索方向和步长，得到各种不同的算法.
由此，先考虑如何确定k α.设从k x 出发，已经选定下山方向k p .令 ()()k k f x p αϕα=+
()()()2T
T k k k k k k x p A x p b x p ααα=++-+
()22T
T k k k k k p Ap r p x ααϕ=-+，
其中k k r b Ap =-.由一元函数极值存在的必要条件有
()220T
T k k k k f p Ap r p αα'=-=
所确定的α即为所求步长k α，即 T k k
k T
k k
r p p Ap α=. 步长确定后，即可算出
1k k k k x x p α+=+.
此时，只要0T k k r p ≠，就有
()()()()1k k k k k k x x x p x ϕϕϕαϕ+-=+-
()
2
220T k k T
T k k
k k k k T
k k
r p p Ap r p p Ap αα=-=-<
即()()1k k x x ϕϕ+<.
再考虑如何确定下山方向k p .易知负梯度方向是()x ϕ减小最快的方向，但简单分析就会发现负梯度方向只是局部最佳的下山方向，而从整体来看并非最佳.故采用新的方法寻求更好的下山方向——共轭梯度法. 下面给出共轭梯度法的具体计算过程:
给定初始向量0x ，第一步仍选用负梯度方向为下山方向，即00p r =，于是有
00
010001000
,,T T r r x x p r b Ax p Ap αα==+=-.
对以后各步，例如第k+1步（k ≥1）,下山方向不再取k r ，而是在过点由向量k r 和1k p -所张成的二维平面
21{|,,}k k k x x x r p R πξηξη-==++∈
内找出使函数ϕ下降最快的方向作为新的下山方向k p .考虑ϕ在2π上的限制：
()1,()k k k x r p ψξηϕξη-=++
11()()T k k k k k k x r p A x r p ξηξη--=++++
12()T k k k b x r p ξη--++. 计算ψ关于,ξη的偏导得： ()()11112,
2,T T T k k k k k k
T T k k k k r Ar r Ap r r r Ap p Ap ψξηξ
ψξηη
----∂=+-∂∂=+∂
其中最后一式用到了10T k k r p -=，这可由k r 的定义直接验证.令 0ψψ
ξ
η
∂∂==∂∂， 即知ϕ在2π内有唯一的极小值点
001k k k x x r p ξη-=++%，
其中0ξ和0η满足 00101011,
0.
T T T k k k k k k T T
k k k k r Ar r Ap r r r Ap p Ap ξηξη----⎧+=⎨+=⎩
由于0k r ≠必有00ξ≠，所以可取
()0
1
001
k k k k p x
x r p ηξξ-=
-=+%
作为新的下山方向.显然，这是在平面2π内可得的最佳下山方向.令010
k η
βξ-=，则可
得
1
111.T k k k T k k r Ap p Ap β----=-
注：这样确定的k p 满足10T
k
k p Ap -=，即k p 与1k p -是相互共轭的. 总结上面的讨论，可得如下的计算公式：
T k k
k T
k k
r p p Ap α= ， 1k k k k x x p α+=+， 11k k r b Ax ++=-，
1T k k
k T
k k
r Ap p Ap β+=-， 11k k k k p r p β++=+. 在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化1k r +的计算公式：
11()k k k k k k k k r b Ax b A x p r Ap αα++=-=-+=-.
因为k Ap 在计算k α是已经求出，所以计算1k r +时可以不必将1k x +代入方程计算，而是从递推关系1k k k r b Ap α+=-得到.
再来简化k α和k β的计算公式.此处需要用到关系式
1110,T T T k k k k k k r r r p r p +-+=== 1,2,k =K .
从而可导出
1111
,T T k k k k r r r α+++=-， ()111T
T T k k k k k k k k k p Ap p r r p r αα+=-=
()1111T T
k k k k k k k k
r r p r r βαα--=+=.
由此可得
,T k k k T k k r r p Ap α=， 11.T k k k T k k
r r r r β++=.
从而有求解对称正定方程组的共轭梯度法算法如下：
0x =初值
00r b Ax =-；0k =
while 0k r ≠
1k k =+ if 1k = 00p r =
