m5cg.m

求解大规模对称矩阵的CG算法

clear all;
format long;  %修改默认数据类型
tic;   %开始计时
n=1000;   %矩阵规模
e=0.0000001;  %精度
A=hilb(n);   %生成系数矩阵
K=cond(A);  %条件数
M=triu(A,-1);
M=tril(M,1);
x=ones(n,1);   %生成精确解
b=A*x;   %生成精确解对应的b
X=zeros(n,1);  %初始
ek1=sqrt((X-x)'*A*(X-x)); %开始是不会知道精确解的！
sulv=0;  % 评价指标
r=b-A*X;
ni=inv(M);
z=ni*r;
p=z;
k=1;%迭代次数
fan=norm(r);
while (k<n)&&(fan>e)
zr=r'*z;
afa=zr/(p'*(A*p));
X=X+afa*p;
ek2=sqrt((X-x)'*A*(X-x));
y(k,1)=ek2/ek1;
sulv=sulv+y(k,1);
ek1=ek2;
r=r-afa*(A*p);
z=ni*r;
fan=norm(r);
beta=r'*z/zr;
p=z+beta*p;
k=k+1;
end
en=norm(X-x);
sulv=sulv/(k-1);
B=A*X;
xl=1:k-1;
toc;
plot(xl,y),xlabel('迭代步'),ylabel('残差速率'),title('下降速率随迭代步的变化')