篇一：共轭梯度实验报告
共轭梯度法求解线性方程组
一.共轭梯度法
共轭梯度法是求解正定系数矩阵线性方程组
Ax?b(A为对称正定矩阵)1
的一种迭代方法，该方法的好处是可以用于求解病态问题，大大减少了计算误差。 求解方程组（1）的本质是求二次函数
TT?(x)?1xAx?bx2
的极小值点。由于极值的必要条件，在极小点x*处，函数?(x)的梯度满足
??(x)?Ax?b?03
即极小点x*是方程（1）的解，因此线性方程组的解转化为一个求极值问题。
二.共轭的几何意义
从本质上来讲，共轭的方向就是某种意义下的正交方向。当A为单位阵时，共轭就是通常意义下的正交。为了便于理解，考虑二维（n?2）情况，函数
TT??1
xx?bx3
的等值线为一族同心圆，而式（2）表示的函数?(x)的等值线是一族同心椭圆。从直观上来讲，圆的切线与切点到圆心的连线正交，那么，椭圆的切线与切点到椭圆中心的连线是共轭。图1给出了共轭的几何意义。
（a） 正交方向（b）共轭方向
图 1(转载自：www.xiaocaOfaNWen.com 小草 范 文 网:共轭梯度实验报告) 正交与共轭的几何意义
从图1可以看出，如果能找到两个共轭方向，那么沿一个方向前进，找到极小点后，再沿另一个方向继续前进，找到极小点，这个极小点就是椭圆的重心，即二次函数的极小点。
对于n维问题，只要能找到n个相互共轭的非零方向，从某一点出发，沿每个方向前进，找到一维极小点，然后再从这个点出发，沿第2个方向前进找极小点，当n个方向都完成这一动作后，就能够得到n维二次函数的极小点。
三.共轭梯度算法程序流程
1. 输入数据A,b ，置r(0)?b?Ax(0)，d(0)?r(0)，精度要求?和k?0；
2. 计算
(r(k))Tr(k)
?k?(k)T(k),x(k?1)?x(k)??kd(k),r(k?1)?r(k)??kAd(k); (d)Ad
3. 若x(k?1)?x(k)??或k?n?1,则停止计算（x(k?1)作为方程组的解）；否则，计算
(r(k?1))Tr(k?1)
(k?1)(k?1)(k)?k?1?,d?r??d; k?1(k)T(k)(r)r
4. 置k?k?1，转第2步。
四.共轭梯度算法的调用及应用
1. 程序内容
function [x,k] = cgg(A,b,ep )
b=b';
if nargin<3,ep=1e-5;end
n=length(A);
x=zeros(n,1);
r=b;d=r;k=0;
while k<n
alpha=(r'*r)/(d'*A*d);
r1=r;
s=alpha*d;
x=x+s;
r=r-A*s;
if norm(s)<ep break;end
beta=(r'*r)/(r1'*r1);
d=r+beta*d;
k=k+1;
end
其中输入项：A为待求解方程组的系数矩阵，b为方程左端常数项，ep为计算精度，输出项：x为方程组的解，k为迭代次数。直接在matlab命令框里输入命令
[x,k] = cgg(A,b,ep )
即可。
2. 实际应用
用共轭梯度法求如下解线性方程组，Ax?b，其中
?2?1???1??????12?1???0??,b??? A???????12?1???0????1??12?????
其中A的阶数n分别取100,200,400，指出计算结果是否可靠。
解：在matlab里产生相关系数矩阵，将精度设置为1e-9,调用共轭算法得出100阶，200阶，400阶的方程组解和对应的迭代次数为
x100?(?1,?1,
x200?(?1,?1,
x400?(?1,?1,,?1,?1)T,k?50,?1,?1)T,k?100 ,?1,?1)T,k?200
由计算结果可知计算结果是可靠的，用共轭梯度算法求解时，随着方程组的阶数增大，方程组的解是稳定的，而且迭代次数只为对应方程组阶数的一半，远小于方程组的阶数，计算速度很快。
篇二：共轭梯度法实验报告
数值代数实验报告
一、实验名称：用共轭梯度法解线性方程组。
二、实验目的：进一步熟悉理解掌握共轭梯度法解法思路，提高matlab编程能力。 三、实验要求：已知线性方程矩阵，应用共轭梯度法在相关软件编程求解线性方程组的解。 四、实验原理：
1．共轭梯度法：
考虑线性方程组
Ax?b
的求解问题，其中A是给定的n阶对称正定矩阵，b是给定的n维向量，x是待求解的n维向量.为此,定义二次泛函
?(x)?xTAx?2bTx.
定理1 设A对称正定,求方程组Ax?b的解，等价于求二次泛函?(x)的极小值点. 定理1表明，求解线性方程组问题就转化为求二次泛函?(x)的极小值点问题.求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量x0，确定一个下山方向p0，沿着经过点x0而方向为p0的直线x?x0??p0找一个点
x1?x0??0p0，
使得对所有实数?有
??x0??0p0????x0??p0?，
即在这条直线上x1使?(x)达到极小.然后从x1出发，再确定一个下山的方向p1，沿着直
线x?x1??p1再跨出一步，即找到?1使得??x?在x2?x1??1p1达到极小：
??x1??1p1????x1??p1?.
重复此步骤，得到一串
?0,?1,?2,
x?xk??pk上确定步长?k使
和 p0,p1,p2,
，
称pk为搜索方向，?k为步长.一般情况下，先在xk点找下山方向pk，再在直线
??xk??kpk????xk??pk?,
最后求出xk?1?xk??kpk.然而对不同的搜索方向和步长，得到各种不同的算法.
由此，先考虑如何确定?k.设从xk出发，已经选定下山方向pk.令 f??????xk??pk?
??xk??pk?A?xk??pk??2bT?xk??pk?
T
??2pkApk?2?rkTpk???xk?，
T
其中rk?b?Apk.由一元函数极值存在的必要条件有
T
f?????2?pkApk?2rkTpk?0
所确定的?即为所求步长?k，即 步长确定后，即可算出
此时，只要rkTpk?0，就有
rkTpk
. ?k?T
pkApk
xk?1?xk??kpk.
??xk?1????xk????xk??kpk????xk?
T
Apk?2?krkTpk?? ??k2pk
?rkTpk?
2
即??xk?1????xk?.
T
pkApk
?0
再考虑如何确定下山方向pk.易知负梯度方向是?(x)减小最快的方向，但简单分析就会发现负梯度方向只是局部最佳的下山方向，而从整体来看并非最佳.故采用新的方法寻求更好的下山方向——共轭梯度法. 下面给出共轭梯度法的具体计算过程:
给定初始向量x0，第一步仍选用负梯度方向为下山方向，即p0?r0，于是有
r0Tr0
?0?T,x1?x0??0p0,r1?b?Ax0.
p0Ap0
对以后各步，例如第k+1步（k?1）,下山方向不再取rk，而是在过点由向量rk和pk?1所张成的二维平面
?2?{x|x?xk??rk??pk?1,?,??R}
内找出使函数?下降最快的方向作为新的下山方向pk.考虑?在?2上的限制：
???,????(xk??rk??pk?1)
?(xk??rk??pk?1)TA(xk??rk??pk?1)
?2bT(xk??rk??pk?1). ??计算?关于?,??2??rTAr??rTAp?rTr?,
kkkk?1kk
??
??TTT?2?rAp??p?kk?1k?1Apk?1?,r其中最后一式用到了rkpk?1，这可由的定义直接验证.令 0k??
????
??0， 即知?在?2内有唯一的极小值点
??
??
x?xk??0rk??0pk?1，
其中?0和?0满足
??0rkTArk??0rkTApk?1?rkTrk,
?TT
??0rkApk?1??0pk?1Apk?1?0.
1
由于rk?0必有?0?0，所以可取
pk?
作为新的下山方向.显然，这是在平面?2内可得的最佳下山方向.令?k?1?得
?0
?x?xk??rk?
?0
p ?0k?1
?0
，则可?0
rkTApk?1
?k?1??T.
pk?1Apk?1
T
注：这样确定的pk满足pkApk?1?0，即pk与pk?1是相互共轭的.
总结上面的讨论，可得如下的计算公式：
rkTpk
，xk?1?xk??kpk， ?k?T
pkApk
rk?1?b?Axk?1，
rkT?1Apk
，pk?1?rk?1??kpk. ?k??T
pkApk
在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化rk?1的计算公式：
rk?1?b?Axk?1?b?A(xk??kpk)?rk??kApk.
因为Apk在计算?k是已经求出，所以计算rk?1时可以不必将xk?1代入方程计算，而是从递推关系rk?1?b??kApk得到.
再来简化?k和?k的计算公式.此处需要用到关系式
rkTrk?1?rkTpk?1?rkT?1pk?0, k?1,2,
从而可导出
1
rkT?1??rkT?1rk?1,， 1T?k1TT
pkApk?pk?rk?rk?1??pkrk
??k1Tk1
?rk?rk??k?1pk?1??rkTrk.
.
由此可得
?k?k
rkTrkrkT?1rk?1?k?T,，?k?T..
pkApkrkrk
从而有求解对称正定方程组的共轭梯度法算法如下： x0?初值
r0?b?Ax0；k?0 while rk?0
k?k?1
if k?1p0?r0
else
?k?2?rkT?1rk?1rkT?2rk?2pk?1?rk?1??k?2pk?2 end
T
?k?1?rkT?1rk?1pk?1Apk?1
xk?xk?1??k?1pk?1 rk?rk?1??k?1Apk?1 end
x?xk
注：该算法每迭代一次仅需要使用系数矩阵A做一次矩阵向量积运算. 定理2由共轭梯度法得到的向量组?ri?和?pi?具有如下基本性质： （1）piTrj?0，0?i?j?k; （2）riTrj?0， i?j，0?i,j?k; （3）piTApj?0， i?j，0?i,j?k; （4）span{r0,其中
,rk}?span{p0,,pk}??(A,r0,k?1)，
?(A,r0,k?1)?span{r0,Ar0,,Akr0}，
通常称之为Krylov子空间.下面给出共轭梯度法全局最优性定理： 定理3 用共轭梯度法计算得到的近似解xk满足
??xk??min???x?:x?x0??(A,r0,k)?
或
xk?x*
A
?min?x?x*
A
:x?x0??(A,r0,k)?，
其中x
A
?x*是方程组Ax?b的解，?(A,r0,k)是由所定义的Krylov子空间.
定理2表明，向量组r0,,rk和p0,,pk分别是Krylov子空间?(A,r0,k?1)的正
交基和共轭正交基.由此可知，共轭梯度法最多n步便可得到方程组的解x*.因此，理论上来讲，共轭梯度法是直接法.
然而实际使用时，由于误差的出现，使rk之间的正交性很快损失，以致于其有限步终止性已不再成立.此外，在实际应用共轭梯度法时，由于一般n很大，以至于迭代O?n?次所耗费的计算时间就已经使用户无法接受了.因此，实际上将共轭梯度法作为一种迭代法使用，而且通常是rk是否已经很小及迭代次数是否已经达到最大允许的迭代次数kmax来终止迭代.从而得到解对称正定线性方程组的实用共轭梯度法，其算法如下：
x?初值
k?0;r?b?Ax;??rTr
while
??b
2
?and?k?k
max
?
k?k?1
if k?1p?relse
????;p?r??pend
??Ap;???pT?;x?x??p r?r???;???;??rTr end
算法中，系数矩阵A的作用仅仅是用来由已知向量p产生向量??Ap，这不仅可以充分利用A的稀疏性，而且对某些提供矩阵A较为困难而由已知向量p产生向量
??Ap又十分方便的应用问题是十分有益的。
2．算例
?10 1 2 3 4??x1??12??1 9 -1 2 -3??x???27?运用共轭梯度法求解下述方程，并解释你所观察到的结果 ???2????2 -1 73-5??x3???14? ??????3 2312 -1?17x解：已知共轭梯度法的???4??MATLAB?程序代码如下所示：
??12???4 -3-5 -1 15?????x5??function[x,n]=conjgrad(A,b,x0) %采用共轭梯度法求线性方程组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);
篇三：运筹学实验报告(F-R共轭梯度法、Wolfe简约梯度法)
一、实验目的:
1、掌握求解无约束最优化问题的 F-R共轭梯度法，以及约束最优化问题 Wolfe简约梯度法。
2、学会用MATLAB编程求解问题，并对以上方法的计算过程和结果进行分析。
二、实验原理与步骤:1、F-R共轭梯度法
基本步骤是在点X 处选取搜索方向d, 使其与前一次
(k)
(k)
的搜索方向d
(k?1)
关于A共轭，即
d(k?1)?d(k),Ad(k?1)??0
(k)(k)(k?1)
然后从点X 出发，沿方向d求得f(X)的极小值点X ,
即
f(X(k?1))?minf(X(k)??d)
??0
(k)
如此下去, 得到序列{X}。不难求得?d
(k)
,Ad(k?1)??0的解为
(k)
X
(k?1)
?X
(k)
?b?AX(k),d(k)?(k)?d(k?1)(k?1)
?d,Ad?
(k)
d注意到的选取不唯一，我们可取
d(k)???f(X(k))??k?1d(k?1)
(k)(k?1)
由共轭的定义?d,Ad??0可得：
?r(k),Ad(k?1)?
?k?1??
?d(k?1),Ad(k?1)?
共轭梯度法的计算过程如下：
第一步：取初始向量X, 计算
?d(0)?r(0)???f(X(0))?b?AX(0)?
?r(0),Ad(0)??
??0??(0)(0)
?d,Ad??
?X(1)?X(0)??d(0)
0?
(0)
第k?1步：计算
?r(k)???f(X(k))?b?AX(k)
?(k)(k?1)
??????r,Ad
?k?1?d(k?1),Ad(k?1)???
?d(k)?r(k)??k?1d(k?1)?(k)(k)
?r,Ad???k??
??d(k),Ad(k)??(k?1)(k)(k)X?X??d?0?
2、Wolfe简约梯度法
Wolfe基本计算步骤：
第一步：取初始可行点x0∈Xl ,给定终止误差ε>0 ，令k:=0；
k
第二步：设IB是xk的m个最大分量的下标集，对矩阵A进行相应
分解 A=(Bk ,Nk)；
第三步：计算 ?f x =
k
?Bf xk
?f xk N
,然后计算简约梯度
?1k
rN=?(BkNk)T?Bf xk +?Nf xk ；
第四步：构造可行下降方向 pk. 若||pk||≤ε ，停止迭代，输
出xk。否则进行第五步。
第五步：进行有效一维搜索，求解min
f xk+tpk ,得到最优
解tk. 令xk+1=xk+ tkpk , k:=k+1, 转入第二步。
三、实验内容：
1、（运筹学P153页第20题）用F-R法求解
min(1?x1)2+2(x2?x12)2
选取初始点x0=(0,0)T, ε=10?6.
2、（运筹学P154页第25题）用Wolfe法求解以下问题： minf x1,x2 =2x12+2x22?2x1x2?4x1?6x2
s.t.x1+x2≤2
x1+5x2≤5 x1≥0,x2≥0选取初始可行点x0=(0,0)T, ε=10?6.
四、问题求解： 问题1求解：（F-R法） 程序代码如下： （1）主函数
syms x1 x2 r;
f=(1-x1)^2+2*(x2-x1^2)^2; x=[x1,x2];
df=jacobian(f,x); df=df.';
error=0.000001; x0=[0,0]';
g1=subs(df,x,x0); k=0;
while(norm(g1)>error) if k==0
d=-g1; else
bta=g1'*g1/(g0'*g0);
d=-g1+bta*d0; end
y=subs(f,x,x0+r*d); result=jintuifa(y,r);
result2=golden(y,r,result); step=result2; x0=x0+step*d;
g0=g1;g1=subs(df,x,x0); d0=d;k=k+1; end; k x0
（2）子函数
进退法确定一维搜索区间：
function result=jintuifa(y,r) t0=0; step=0.0125; t1=t0+step;
ft0=subs(y,{r},{t0}); ft1=subs(y,{r},{t1}); if(ft1<=ft0) step=2*step; t2=t1+step;
ft2=subs(y,{r},{t2}); while(ft1>ft2)
t1=t2;step=2*step; t2=t1+step;
ft1=subs(y,{r},{t1}); ft2=subs(y,{r},{t2}); end else
step=step/2;t2=t1;t1=t2-step; ft1=subs(y,{r},{t1}); while(ft1>ft0)
step=step/2;t2=t1;t1=t2-step; ft1=subs(y,{r},{t1}); end end
result=[t2];
黄金分割法进行一维搜索：
function result=golden(y,r,m)