%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution to Problem 2 and Problem 5:
% 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')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For Problem 5, solve linear systems with perturbed right-hand sides
% and compare with the result predicted by the confidence intervals 
% for the pertubation of the % solution.  
% We use kappa = 1.96 to get a 95% confidence interval on
% each component of x.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

   kappa = 1.96;
   u = A\[1;0];
   deltax1 = kappa * norm(u)*tau;
   u = A\[0;1];
   deltax2 = kappa * norm(u)*tau;
   disp(sprintf('Confidence intervals for perturbation of'))
   disp(sprintf('Right-hand side for Linear System %d',iexample))
   disp(sprintf(' %f <= x_1 <= %f', xtrue(1)-deltax1,xtrue(1)+deltax1))
   disp(sprintf(' %f <= x_2 <= %f', xtrue(2)-deltax2,xtrue(2)+deltax2))
   disp(sprintf(' %f of the x_1 values lie in the confidence interval',...
            sum( abs(forward(:,1))<=deltax1 ) /nprob))
   disp(sprintf(' %f of the x_2 values lie in the confidence interval',...
            sum( abs(forward(:,2))<=deltax2 ) /nprob))

end

figure(1)
print -depsc ls1_eqns.eps
figure(2)
print -depsc ls1_forw.eps
figure(3)
print -depsc ls1_back.eps
figure(4)
print -depsc ls2_eqns.eps
figure(5)
print -depsc ls2_forw.eps
figure(6)
print -depsc ls2_back.eps