% Heath: Computer Problem 2.9
%-------------------------------------------------------------
% (a) Use Gaussian elimination without pivoting to solve
%      [e 1][x1] = [1+e]
%      [1 1][x2]   [ 2 ]
% for e = 10^(-2k), k = 1,...,10.  Compare with the true
% solution x* = [1 1]'.
%-------------------------------------------------------------
disp('Part (a)')
xs = [1 1]';      % true solution

for k=1:10;
   e = 10^(-2*k);
   A = [e 1; 1 1]; b = [1+e 2]';
   L = [1 0; 1/e 1];
   U = [e 1; 0 1-1/e];
   
   x = U\(L\b);
   fprintf('e = %8.1e   error = %8.1e  c(A) = %8.1e  c(L) = %8.1e  c(U) = %8.1e\n',...
      e,norm(x-xs),cond(A),cond(L),cond(U))
end;

%-------------------------------------------------------------
% (b) Now include one iteration of iterative refinement.
%-------------------------------------------------------------
disp('Part (b): with iterative refinement')

for k=1:10;
   e = 10^(-2*k);
   A = [e 1; 1 1]; b = [1+e 2]';
   L = [1 0; 1/e 1];
   U = [e 1; 0 1-1/e];
   
   x = U\(L\b);
   
   % iterative refinement (without extended precision)
   
   r = b - A*x;
   w = U\(L\r);
   x = x + w;
   fprintf('e = %8.1e   error = %8.1e\n',e,norm(x-xs))
end;