[分享]共轭斜量法matlab程序

风花雪月 · 发表于 2005-11-22 20:47

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

共轭斜量法数学表达式：
　<v:shapetype> <v:stroke joinstyle="miter"></v:stroke><v:formulas><v:f eqn="if lineDrawn pixelLineWidth 0"></v:f><v:f eqn="sum @0 1 0"></v:f><v:f eqn="sum 0 0 @1"></v:f><v:f eqn="prod @2 1 2"></v:f><v:f eqn="prod @3 21600 pixelWidth"></v:f><v:f eqn="prod @3 21600 pixelHeight"></v:f><v:f eqn="sum @0 0 1"></v:f><v:f eqn="prod @6 1 2"></v:f><v:f eqn="prod @7 21600 pixelWidth"></v:f><v:f eqn="sum @8 21600 0"></v:f><v:f eqn="prod @7 21600 pixelHeight"></v:f><v:f eqn="sum @10 21600 0"></v:f></v:formulas><v:path connecttype="rect" gradientshapeok="t" extrusionok="f"></v:path><lock aspectratio="t" v:ext="edit"></lock></v:shapetype><v:shape><v:imagedata></v:imagedata></v:shape>任给初始向量x(0),
r(0)=b-Ax(0),p(0)=r(0)
αk=(r(k),r(k))/(AP(k),P(k)) x(k+1)=x(k)+αkP(k)
r(k+1)=r(k)-αkAP(k) βk=(r(k+1),r(k+1))/(r(k),r(k))
P(k+1)=r(k+1)+βkP(k) (k=0,1,2,…..)
%%matlab程序
%%共轭斜量法，程序中变量alpha是αk，beta是βk
%%x为返回值，a为系数矩阵b列向量
%%函数中添加了计时函数tic,toc也可以不要
%%最多迭代次数为50次
function x=cg(a,b)
[m,n]=size(b);%判断输入的b是行向量还是列向量，如果是行向量，转化成列向量
if m<n
 b=b';
end
x=zeros(length(a),1);%初始迭代值
tic%迭代开计时
r0=b-a*x;
p=r0;
r=r0;
alpha=dot(r0,r0)/dot(a*p,p);
for k=1:50%设置最大迭代次数
 x=x+alpha*p;
 r=r-alpha*a*p;
 if dot(r,r)<0.0001
 fprintf('迭代误差为:%e\n',dot(r,r));
 fprintf('迭代次数为:%g\n',k);
 t=toc;
 fprintf('所用时间为:%g\n',t);
 break;
 end
 beta=dot(r,r)/dot(r0,r0);
 p=r+beta*p;
 alpha=dot(r,r)/dot(a*p,p);
 r0=r;
 if k>49%超过最大次数，显示警告信息
 disp('已超过最大迭代次数');
 t=toc;
 fprintf('所用时间为:%g\n',t);
 break;
 end
end

账号		自动登录	找回密码
密码			我要加入