共轭梯度法实验报告
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
数值代数实验报告
一、实验名称:用共轭梯度法解线性方程组。
二、实验目的:进一步熟悉理解掌握共轭梯度法解法思路,提高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 =