%Arnoldi方法
n=1000  %方程階數(shù)
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;
r0=b-A*X0;
beita=norm(r0);
V(:,1)=r0/beita;
for m=1:1000
    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
     e1=e;
     e1(1)=1;
     ym=H(1:m,1:m)\(beita*e1(1:m));
     if norm(ym)<1.0e5
        if H(m+1,m)==0
          X=X0+V(:,1:m)*ym;
          break;
         end
        if H(m+1,m)~=0
            m
            em=e;
            em(m)=1;
           nr(m)=H(m+1,m)*norm(em(1:m)'*ym);
          if  norm(nr(m))/norm(b)<=precision_control
              X=X0+V(:,1:m)*ym
              break;
          end
        end
     end
 end
nr(1:m)
plot(nr(1:m),'-ro');
xlabel('迭代步數(shù)')
ylabel('殘差')