% computes solution of P16 on A#7 in Math 694

N = 100;
A = [P16map(N,[1 0 0]) P16map(N,[0 1 0]) P16map(N,[0 0 1])];
close all % close figures that appear
% size(A) = 100 x 3

disp('answer to (c) is row 50 of A:')
A(50,:)

h = 1/N;
xj = (0:h:1-h)';
b = 0.3*ones(N,1) + 0.05*sin(pi*(xj+h));
% size(b) = 100 x 1

[Q,R] = qr(A,0);
z = Q' * b;       % 3 x 1
y = Q * Q' * b;   % = P b; 100 x 1

disp('answer to (d) is vector c of coefficients:')
c = R \ z   % called "x" in general analysis


% show solution to least squares problem as a heat distribution
plot(xj,b,'o-',xj,y,'*-')
legend('b','y = A c = P b')


disp('additional information related to conditioning of least squares problem:')
theta = acos(norm(y)/norm(b));   % theta =  0.016265
eta = norm(A)*norm(c)/norm(y);   % eta   =  1.0904
kappa = cond(A);                 % kappa =  58.265   = sigma_1 / sigma_3
disp(['[theta eta kappa] = ' num2str([theta eta kappa])])

% see table in theorem 18.1 in Trefethen & Bau
table = [1/cos(theta),      kappa / (eta * cos(theta));
         kappa/cos(theta),  kappa + (kappa * kappa * tan(theta)) / eta] 
%table =
%     1.0001    53.4401
%    58.2725   108.9059