% Script File InverseMatrixError
%
% Plots, as a function of n, the relative error in using the matrix identity 
%          Inverse(X)   =   sum of (I-X)k, k>=0
%
%
% nTerms = 150
% A=[0.10 0.10 0.10; 0.20 0.75 0.05; 0.10 0.30 0.60]
% B=[0.90 -0.10 0.00; -0.20 0.25 -0.05; -0.10 -0.30 0.40]

   close all
   nTerms = 150;
   I=[1 0 0; 0 1 0; 0 0 1]
   A=[0.10 0.10 0.10; 0.20 0.75 0.05; 0.10 0.30 0.60]
   B=[0.90 -0.10 0.00; -0.20 0.25 -0.05; -0.10 -0.30 0.40]
   
   A_sum=zeros(3,3);
   B_sum=zeros(3,3);
   A_relerr=ones(nTerms,1);
   B_relerr=ones(nTerms,1);
   for n=0:(nTerms-1)
      A_sum=A_sum+(I-A)^n;
      B_sum=B_sum+(I-B)^n;
      A_relerr(n+1)=norm(A_sum-inv(A))/norm(inv(A));
      B_relerr(n+1)=norm(B_sum-inv(B))/norm(inv(B));
   end
   subplot(2,1,1), semilogy(0:nTerms-1,A_relerr,'r-')
   ylabel('Relative Error in Partial Sum.')
   xlabel('Order of Partial Sum.')
   title(sprintf('Matrix A'))
   subplot(2,1,2), semilogy(0:nTerms-1,B_relerr,'b-')
   ylabel('Relative Error in Partial Sum.')
   xlabel('Order of Partial Sum.')
   title(sprintf('Matrix B'))