% discussed Gauss elimination, including the factorization A = LU and how
% solve  Ax = b, by first solving Ly = b (by forward substitution) and then
% solving Ux = y (by backward substitution).
%
% Then discussed pivoting, i.e., the choice of a pivot row for each unknown
% The book gives an example that illustrates a possible diffculty. Here is
% a modified version of its script NoPivot:
>> type nopivot

% Script File: NoPivot
% Examines solution to
%
%        [ delta 1 ; 1 1][x1;x2] = [1+delta;2]
%
% for a sequence of diminishing delta values.
  %   clc
disp(' Delta            x(1)                   x(2)  ' )
disp('-----------------------------------------------------')
for delta = logspace(-2,-18,9)
   if exist('goodmult')
      switch goodmult
      case 1 
         A = [1 1/delta; 1 1];
         b = [(1+delta)/delta; 2];
      case 2
         A = [delta 1; delta  delta];
         b = [1+delta; 2*delta];
      end
   else 
      A = [delta 1; 1 1];
      b = [1+delta; 2];
   end
   L = [ 1 0; A(2,1)/A(1,1) 1];
   U = [ A(1,1) A(1,2) ; 0 A(2,2)-L(2,1)*A(1,2)];
   y(1) = b(1);
   y(2) = b(2) - L(2,1)*y(1);
   x(2) = y(2)/U(2,2);
   x(1) = (y(1) - U(1,2)*x(2))/U(1,1);
   disp(sprintf(' %5.0e   %20.15f  %20.15f',delta,x(1),x(2)))
end

% I run it first as it is in the book, i.e., with goodmult undefined:
>> nopivot
 Delta            x(1)                   x(2)  
-----------------------------------------------------
 1e-002      1.000000000000001     1.000000000000000
 1e-004      0.999999999999890     1.000000000000000
 1e-006      1.000000000028756     1.000000000000000
 1e-008      0.999999993922529     1.000000000000000
 1e-010      1.000000082740371     1.000000000000000
 1e-012      0.999866855977416     1.000000000000000
 1e-014      0.999200722162641     1.000000000000000
 1e-016      2.220446049250313     1.000000000000000
 1e-018      0.000000000000000     1.000000000000000

% the output looks as in the book. Next, I run it again, but now for the 
% modified linear system, in which I have multiplied the first equation 
% by 1/delta, hence have removed the `large multiplier' that the book claims
% is the problem here: 
>> goodmult = 1;
>> nopivot
 Delta            x(1)                   x(2)  
-----------------------------------------------------
 1e-002      1.000000000000000     1.000000000000000
 1e-004      1.000000000000000     1.000000000000000
 1e-006      1.000000000000000     1.000000000000000
 1e-008      1.000000000000000     1.000000000000000
 1e-010      1.000000000000000     1.000000000000000
 1e-012      0.999877929687500     1.000000000000000
 1e-014      1.000000000000000     1.000000000000000
 1e-016      2.000000000000000     1.000000000000000
 1e-018      0.000000000000000     1.000000000000000

%  It looks less perturbed, but, still, there's similar trouble by the time
%  delta = 1e-12 and less.
% However, just to drive this point home, I'll also remove the large multiplier
% by leaving the first equation as is, but multiplying the second equation by
% delta, i.e., use nopivot with ...
>> goodmult=2;
>> nopivot
Delta            x(1)                   x(2)  
-----------------------------------------------------
 1e-002      1.000000000000001     1.000000000000000
 1e-004      0.999999999999890     1.000000000000000
 1e-006      0.999999999917733     1.000000000000000
 1e-008      0.999999993922529     1.000000000000000
 1e-010      1.000000082740371     1.000000000000000
 1e-012      0.999866855977416     1.000000000000000
 1e-014      0.999200722162641     1.000000000000000
 1e-016      1.110223024625157     1.000000000000000
 1e-018      0.000000000000000     1.000000000000000

% ... and we see the very same trouble as before, even though the offending
% multiplier is now 1.
%
% Conclusion: the difficulty is NOT the large multiplier per se. 
% The difficulty comes from using the first equation as pivot equation.
% For small delta, that first equation has trouble pinpointing x_1 for given
% x_2 since it is almost parallel to the straight line y = x_2.
% See the notes on pivoting and the pictures there.
>> diary off
