数值代数实验报告
一、实验名称：用共轭梯度法解线性方程组。
二、实验目的：进一步熟悉理解掌握共轭梯度法解法思路，提高matlab编程能力。
三、实验要求：已知线性方程矩阵，应用共轭梯度法在相关软件编程求解线性方程组的解。
四、实验原理：
1．共轭梯度法：
考虑线性方程组
b?Ax的求解问题，其中A是给定的n阶对称正定矩阵，b是给定的n维向量，是待求解x的n维向量.为此,定义二次泛函
TT?xbx2Ax(x).
?的极小值点. 求方程组的解，等价于求二次泛函设A对称正定定理1 ,)(xb?Ax?的极小值点问题表明，求解线性方程组问题就转化为求二次泛函. 定理1)x(x，确定求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量0?ppxpx?x?找一个点而方向为一个下山方向 ，沿着经过点的直线00000?px?x?， 0001?有使得对所有实数
px?xp， 00000?xxp，沿着出发，使再确定一个下山的方向达到极小.然后从即在这条直线上)x(111直
pxxp?x?xx达到极线得在小：即再跨出一步，找到 使1112111p?x?px?.
11111重复此步骤，得到一串
,p,p,,,p,， 和202110?xpp，再在直线.一般情况下，先在称点找下山方向为搜索方向，为步长kkkk px?x? 上确定步长使kkk,x?x?p?p kkkkk?pxx.
然而对不同的搜索方向和步长，得到各种不同的算法最后求出.kkk1?k
?.设从出发，已经选定下山方向.由此，先考虑如何确定令 pxkkk? pfxkkTT px?xpx?2pbAkkkkkkT2T ，xpAp?2rpkkkkkr?b?Ap.由一元函数极值存在的必要条件有 其中kkTT 0p?2r?2?pfApkkkk，即即为所求步长所确定的 kTpr?kk. ? kTpApkk 步长确定后，即可算出?px?x?.
kk?1kkT0r?p ，就有此时，只要kk?xpxxx? 2kkkk?1kkTprkk2TTrppAp?2?0 kkkkkkTAppkkxx?.
即kk?1?p减小最快的方向，但简单分.易知负梯度方向是再考虑如何确定下山方向)(xk析就会发现负梯度方向只是局部最佳的下山方向，而从整体来看并非最佳.故采用新的方法寻求更好的下山方向——共轭梯度法.
下面给出共轭梯度法的具体计算过程:
xp?r，于是有给定初始向量 ，第一步仍选用负梯度方向为下山方向，即000Trr00.
Axbp,r,xx? 0011000TpAp00rrp和而是在过点由向量）,下山方向不再取，对以后各步，例如第k+1步（k?1kkk?1所张成的二维平面
Rp,?x?}r?,{?x|x 1k2kk?p上的限制：在.内找出使函数考虑下降最快的方向作为新的下山方向 k2p)r(x, 1kk?kT?p)rxA(x?(rp) 1kkkk?1kk?Tp)r(x?2?b?. 1kkk 的偏导得：关于计算,TTT,rrApArr?2r kkkk?1kk ?TT,prAp?Ap?2Trr0p?，这可由的定义直接验证.令 其中最后一式用到了 1k?1kk?k?11kk?k， ?0 内有唯一的极小值点即知在2p?xxr ，1?kk0k0
满足 和其中00TTTrr,r?rApAr?k01k0?kkkk? TTpAp?rAp?0.?1?k?11?00kkk
00r?，所以可取必有由于 0k?10?r?xpxp kkkk?1 ?000?，则可显然，这是在平面令内可得的最佳下山方向..作为新的下山方向 21?k?0 得TApr?1kk? . 1k?TApp1k?1k?Tppp满足. 与注：这样确定的，即是相互共轭的pAp?0k1k?k?1kk总结上面的讨论，可得如下的计算公式：
Tprkkp?x?x ， ， ? kkkk?1kTpApkkr?b?Ax， 1k?1?kTAprkk1?pr?p?. ， kk?11kk?kTpApkk在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称r的计算公式： .首先来简化的计算公式1kApr)p?b?Ax?b?A(xr?.
kk?1kkk?1kkk?rApx代入方程计算，而是在计算时可以不必将是已经求出，所以计算因为1kk?k1kApb?r?得到.
从递推关系kk?1k的计算公式.此处需要用到关系式和 再来简化kkTTTpprrr?r0, .
2,1,k?k?k1?1k1kkk?从而可导出
1TTr,rr， 1?k1?1?kk?11kTTT?prprr?pAp? k1kkkkkk11kkTT?rrprr? . kkk?1k?k1kkk 由此可得TTrrrr1k?k1?kk.
， ., kkTTpAprrkkkk从而有求解对称正定方程组的共轭梯度法算法如下：
x初值 ?0r?b?Ax； 0k?00r?0 while k 1kkif 1?kp?r 00．
else
TT? rrr?r2?1k?1?kkk?2k?2?prp 2k?1k1kk?2end
TT? Aprp?r1?1kk?k1?k?k?11?p?x?x 1kkk?k11?Apr?r? 1kk?k1k?1?end
xx? k. 注：该算法每迭代一次仅需要使用系数矩阵做一次矩阵向量积运算A? rp 定理2具有如下基本性质： 由共轭梯度法得到的向量组和iiT0rp? ，（1）;0?i?j?kjiT0r?r ，（2）， j?i;0?i,kj?jiT0App? ，） ，（3ji?;kj?0?i,ji?1)k?A,r,{r,,r}?span{p,p}?,(span ）4，（0kk00 其中k?,Ar{rA,r}r(A,,k?1)?span, ， 0000 子空间.下面给出共轭梯度法全局最优性定理：通常称之为Krylovx 用共轭梯度法计算得到的近似解满足定理3 k)x?x?(Ax,?minr,xk: 00k 或 ，),k(A,r:?x?xminxx?x?x?00*k*AA T?x)k(A,r,的解，其中是方程组是由所定义的Krylov子空间. ，Ax?xxb?Ax*0A?,r,k?1),pA(,,pr,r表明，向量组定理2的正Krylov和分别是子空间00k0kx因此，步便可得到方程组的解.交基和共轭正交基.由此可知，共轭梯度法最多n*.
理论上来讲，共轭梯度法是直接法 r之间的正交性很快损失，以致于其有限使然而实际使用时，由于误差的出现，k很大，以至于迭此外，在实际应用共轭梯度法时，由于一般步终止性已不再成立.nnO因此，实际上将共轭梯度法代次所耗费的计算时间就已经使用户无法接受了.r是否已经很小及迭代次数是否已经达到最大允而且通常是作为一种迭代法使用，kk来终止迭代.许的迭代次数从而得到解对称正定线性方程组的实用共轭梯度法，max其算法如下：
初值x?
T? r?r;?Ax0;k?r?b? kandkb while
max2 1k?k? if 1k? rp? else
? ppr?; end
T pxp?x;?Ap;Tr;r?r;r?
end
?，这不仅可以产生向量算法中，系数矩阵的作用仅仅是用来由已知向量pAp?A充分利用的稀疏性，而且对某些提供矩阵较为困难而由已知向量产生向量pAA?又十分方便的应用问题是十分有益的。
Ap?2．算例 x 43 1 102 121运用共轭梯度法求解下述方程，并解释你所观察到的结果 x27-1 2 ?-3 19 2 14 -5?x 2 -17 333 2 3 12 -1?17x解：已知共轭梯度法的MATLAB程序代码如下所示： 415 13 -5 - 12 4-xfunction[x,n]=conjgrad(A,b,x0) 5%采用共轭梯度法求线性方程组Ax=b的解
%线性方程组的系数矩阵：A
%线性方程组中的常数向量：b
%迭代初始向量：x0
%线性方程组的解：x
%求出所需精度的解实际的迭代步数：n
r1=b-A*x0;
p=r1;
n=0;
for i=1:rank(A)
if(dot(p,A*p)<1.0e-50)
break;
end
alpha=dot(r1,r1)/dot(p,A*p);
x=x0+alpha*p;
r2=r1-alpha*A*p;
if(r2<1.0e-50)
break;
end
belta=dot(r2,r2)/dot(r1,r1);
p=r2+belta*p;
n=n+1;
end
用共轭梯度法求解，在MATLAB命令窗口中输入求解程序：
>> A=[10 1 2 3 4 2 4 2 1 1;1 9 -1 2 -3 3 -3 -1 9 2;2 -1 7 3 -5 4 -5 7 -1
3;3 2 3 12 -1 5 -1 3 2 4;4 -3 -5 -1 15 1 3 2 3 2;2 3 4 5 1 10 4 5 3 4;4 -3
-5 -1 3 4 9 1 -3 2;2 -1 7 3 2 5 1 7 -1 1;1 9 -1 2 -3 3 -3 -1 12 10;1 2 3 4
2 4 2 1 10 15];
>> b=[12;-27;14;-17;12;12;-27;14;-17;12];
>> x0=[0;0;0;0;0;0;0;0;0;0];
>> [x,n]=conjgrad(A,b,x0)
输出的计算结果为：
x = 8.5669
-7.1981
-7.3479
-13.9388
1.1303
11.5188
-26.8966
-2.2388
0.0829
13.2786
输出的迭代次数为
n = 10
共轭梯度法理论上只要迭代n步，就能得出方程组的解，但是由于存在计算误差， 步时得到准确解。n即分数向小数转化时存在舍入误差，很难保证在第
3． 总结
本文首先给出最小二乘问题的定义，随后从盲人下山法开始讨论了共轭梯度法的具体推导过程及其相关性质与算法.继而重点给出正则化方法的实用共轭梯度算法并举例进行检验.最后，需要说明虽然正则化方法是求一般线性方程组，b?Axnm?且列满秩的最小二乘解的一种方法且简单易行，但是也有许多不nm?RA?AT形成时运算量大，中某些信息会丢失；当病足之处，如时一般无解；AAnm?AA2T态时其收敛性速度由于很大变得非常之慢等，故为了避免正则化方?))A(A(A22法的缺点，还可运用残量极小化方法或残量正交化方法等更好的方法来解决此类问题.
在实际的科学与工程问题中，常常将问题归结为一个线性方程组的求解问题，而求解线性方程组的数值解法大体上可分为直接法和迭代法两大类.直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法.因此，直接法又称为精确法.迭代法则是采取逐次逼近的方法，亦即从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列,其极限才是方程组的精确解，只经过有限次运算得不到精确解.当线性方程组的系数矩阵为对称正定矩阵时，我们常用共轭梯度法（或简称CG法）求解，目前有关的方法与理论已经相当成熟，并且已成为求解大型稀疏线性方程组最受欢迎的一类方法.