%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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