gmres.m

实现兰乔斯迭代

%---tongyue200612009
%GMres方法
%此m文件直接在malab命令窗命令窗口中输入函数名即可，无需调用函数
n=1000  %方程阶数
m=5;
H=zeros(n);
V=zeros(n);
e=zeros(n,1);
B=randn(n);
D=diag(1:1:n);
A=B*D*inv(B);
b=randn(n,1);
X0=ones(n,1);
nr=ones(n,1);
precision_control=1.0e-6;
for k=1:n
    if k==1
       X0=ones(n,1);
     else
       X0=xm;
     end
       r0=b-A*X0;
       beita=norm(r0);
       V(:,1)=r0/beita;
    for j=1:m
        w=A*V(:,j);
          for i=1:j
          H(i,j)=V(:,i)'*w;
          w=w-H(i,j)*V(:,i);
          end
        H(j+1,j)=norm(w);
        if H(j+1,j)~=0
          V(:,j+1)=w/H(j+1,j);
         end
     end
     em=e;
     em(m)=1;
     HW=[H(1:m,1:m);H(m+1,m)*em(1:m)'];
     e1=e;
     e1(1)=1;
     ym=HW'*HW\(HW'*beita*e1(1:(m+1)));
     xm=X0+V(:,1:m)*ym;
     nr(k)=norm(b-A*xm);
     if  norm(nr(k))/norm(b)<=precision_control
        X=X0+V(:,1:m)*ym
        break;
     end      
end
nr(1:k)
plot(nr(1:k));
xlabel('迭代步数')
ylabel('残差')