%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution to Problem 2 
% Sensitivity Analysis:
% When a Little Means a Lot
%
% Generate two linear systems of equations Ax=b
% where A and b are defined below.
%
% Solve each problem.
% Also, for each problem, solve  nprob  problems formed from replacing
% A by A+E, where the elements of E are normally distributed with mean 0
% and standard deviation tau.
%
% exlinsys.m     Dianne O'Leary 05/2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize tau, nprob, a counter for the figures, and
% a parameter delta used to define A.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tau = .0001;
nprob = 10000;
jfig = 0;
delta = .002;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For each example,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for iexample=1:2,

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Initialize  A and b, and compute the "exact" solution xtrue,
%          and plot the equations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   jfig = jfig + 1;
   figure(jfig)
   if (iexample == 1)
      A = [delta,1;1,1];
      b = [1;0];
      xval=[-2,0];
      yval = (b(1) - A(1,1)*xval)./A(1,2);
      plot(xval,yval,'r')
      hold on
      yval = (b(2) - A(2,1)*xval)./A(2,2);
      plot(xval,yval,'g')
   else
      A = [1+delta,delta-1;delta-1,1+delta];
      b = [2;-2];
      xval=[-5,5];
      yval = (b(1) - A(1,1)*xval)./A(1,2);
      plot(xval,yval,'r')
      hold on
      yval = (b(2) - A(2,1)*xval)./A(2,2);
      plot(xval,yval,'g')
   end
   title(sprintf('Equations for Linear System %d',iexample))
   xlabel('x_1')
   ylabel('x_2')

   xtrue = A \ b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute and display the condition number of the matrix.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   disp(sprintf('The condition number for matrix %d is %f.',iexample,cond(A)))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve perturbed versions, recording 
%       the change in x in array forward
%       and the residual b-Ax_{computed} in array backward.
% Plot the forward and backward errors.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   for i=1:nprob,
      E = tau*randn(2,2);
      x = (A+E)\b;
      forward(i,:) = (x-xtrue)';
      backward(i,:) = (b - A * x)';
   end

   jfig = jfig + 1;
   figure(jfig)
   plot(forward(:,1),forward(:,2),'*')
   title(sprintf('Changes in x for Linear System %d',iexample))
   xlabel('\delta x_1')
   ylabel('\delta x_2')

   jfig = jfig + 1;
   figure(jfig)
   plot(backward(:,1),backward(:,2),'*')
   title(sprintf('Residuals for Linear System %d',iexample))
   xlabel('r_1')
   ylabel('r_2')

end