%% Annotated Matlab transcript -- 4/26
%% CS 323 - Doug DeCarlo

% condition numbers in matlab

% if the 5 in this matrix was 4.9, the matrix would be singular
A = [7 10; 5 7];

cond(A, inf)
cond(A, 2)
cond(A, 1)

% we don't get exactly the identity matrix back...
A * inv(A)

% here it's only invertible because 4.9 isn't represented exactly
A = [7 10; 4.9 7];
inv(A) * A

A = [7 10; 4.9000000000001 7];
inv(A) * A
cond(A,inf)

% --- Jacobi iteration example

% initial guess
x = [0 0 0];
% right hand side
b = [10 19 0];
% actual answer
xans = [1 2 -1];
for i = 1:20
  xn = [1/9*(b(1)-x(2)-x(3))  1/10*(b(2)-2*x(1)-3*x(3))  1/11*(b(3)-3*x(1)-4*x(2))];
  error(i) = norm(xans - xn,inf);
  ratio(i) = norm(xans-xn,inf)/norm(xans-x,inf);
  x = xn;
end

% we're pretty close to ans
xn

% plot error
plot(1:length(error),error)
semilogy(1:length(error),error,'.-')

% ratio converges
plot(1:length(error),ratio)

% --- Gauss Seidel iteration

% initial guess
x1 = 0;
x2 = 0;
x3 = 0;
% right hand side
b = [10 19 0];
% actual answer
xans = [1 2 -1];
for i = 1:20
  x = [x1 x2 x3];
  x1 = 1/9*(b(1)-x2-x3);
  x2 = 1/10*(b(2)-2*x1-3*x3);
  x3 = 1/11*(b(3)-3*x1-4*x2);
  xn = [x1 x2 x3];

  error(i) = norm(xans - xn,inf);
  ratio(i) = norm(xans-xn,inf)/norm(xans-x,inf);
end

% (you get divide by zero errors here because
% it acually converges and the ratio is 0/0)

% plot error
plot(1:length(error),error)
semilogy(1:length(error),error,'.-')

% ratio doesn't actually converge for this...
plot(1:length(error),ratio)