%
% Application Programming in MatLab-environment
% Autumn 2000
% 
% Exercise 2
%

% ----- Problem 1 ----- 

% Let us first clear all variables from memory
% and close all figures.
clear
close all
disp('Problem 1')

% We create five pairs (a,b) of 4x4 matrices with random entries.
a = rand(4,4,5);
b = rand(4,4,5);

% We compute their inverse matrices into the arrays ai and bi,
% and then we check the validity of the matrix formulas.
ai = []; bi = [];
for i=1:5
  ai(:,:,i) = inv(a(:,:,i));
  bi(:,:,i) = inv(b(:,:,i));
  %
  disp((a(:,:,i)*b(:,:,i))' - b(:,:,i)'*a(:,:,i)');
  disp(ai(:,:,i)' - inv(a(:,:,i)'));
  disp(bi(:,:,i)' - inv(b(:,:,i)'));
  %
  disp('The results should be (almost) zero. Press any key.');
  pause;
end


% ----- Problem 2 ----- 
disp('Problem 2')

% Define matrix A
a = [2 1 ; 1 3]
%a = [1 1 ; 1 2]
%a = [2 sqrt(6) ; sqrt(6) 3];
%a = [1 2 ; 2 -2];
%a = [1 2 ; 3 4];

% Compute eigenvalues and -vectors of A
[u,d] = eig(a);

% Check the validity of the formula AU = UD
disp(a*u - u*d)
% Check the validity of the diagonalisation formula
disp(u'*a*u - d)

% Illustrate the basis vectors e1 and e2 and their
% transformations U'*e1 and U'*e2 using the quiver-function

% Remember to set the scale to zero! Otherwise automatic scaling
% is used, which causes division by zero in this special case.
%
% xo and yo include x- and y-coordinates for the starting points 
% of the vectors. Vectors to be drawn are joined to obtain the
% matrix V. Please, pay attention to the use of hold on and hold off.
xo = zeros(1,2); yo = xo;
e1 = [1 ; 0]; e2 = [0 ; 1];
u1 = u'*e1;   u2 = u'*e2;
V  = [e1 e2];
o1 = quiver(xo,yo,V(1,:),V(2,:),0,'b'); hold on
V  = [u1 u2];
o2 = quiver(xo,yo,V(1,:),V(2,:),0,'r'); hold off
title('Basis vectors and their transformations');
xlabel('x-coordinate'); ylabel('y-coordinate');
legend([o1(1) ; o2(1)],'basis vectors','transformed vectors');
clear e1 e2 o1 o2 u1 u2

disp('Press any key to continue.')
pause

% ----- Problem 3 ----- 
disp('Problem 3')

% Construct two-dimensional grid in the rectangle [-4,4]x[-4,4]
xp = -4:0.2:4;
[x,y] = meshgrid(xp);

% We calculate the dimensions of the grid and reshape the 
% matrices x and y into vector form. The two two vectors 
% are joined to obtain the matrix xy.
[n1,n2] = size(x);
n = n1*n2;
xy = [reshape(x,1,n); reshape(y,1,n)];

% Computation of the values of the quadratic form at the grid points.
% At the coordinate  crd = [xc ; yc], the value of the quadratic form
% is given by  crd'*a*crd
%
% First way to compute the desired values:
% zxy = diag(xy' * a * xy);
% - this is not recommended, because it leads to unnecessary 
%   computational operations and memory consumption

% Second way is to use a for-loop:
% for i=1:n, zxy(i)=xy(:,i)'*a*xy(:,i); end
% - this is also a bad choice, because for-loops are slow
%   in MatLab

% Third way is a recommended. It computes the values
% by using operator .* and the function sum
zxy = sum(xy .* (a*xy));

% Then we reshape the computed values, which are in the vector zxy,
% back into matrix form and draw the figures using mesh and contour.
z = reshape(zxy,n1,n2);
mesh(x,y,z)
title('Surface plot of the quadratic form');
xlabel('x-coordinate'); 
ylabel('y-coordinate');
zlabel('value of the quadratic form');

% Draw the contour plot and add the eigenvectors into the figure
% provided they are well-defined.
figure
contour(x,y,z,20); hold on
title('Contour plot of the quadratic form');
xlabel('x-coordinate'); 
ylabel('y-coordinate');

eps = 1.e-6;
if all(diag(d)>eps),
  u1 = u(:,1)./sqrt(d(1,1));
  u2 = u(:,2)./sqrt(d(2,2));
  U  = [u1 u2];
  quiver(xo,yo,U(1,:),U(2,:),0,'r');

  % we add one contour line which corresponds to the 
  % value of the quadratic form at the end points of the 
  % eigenvectors
  v = u1'*a*u1;
  fprintf(1,'\nThe value of the quadratic form at the end points: %5.2f\n',v);
  contour(x,y,z,[v v]); hold off
else
  disp('Not possible to draw the vectors, because');
  disp('the matrix A is not positive definite.');
end

disp('Press any key to continue.');
pause
close all

% Values of the quadratic form with the square root of A
% are computed and illustrated in similar way.
asq = u * sqrt(d) * u'
% alternative way:  asq = sqrtm(a)
zxy = sum(xy .* (asq*xy));
z = reshape(zxy,n1,n2);

mesh(x,y,z)
title('Surface plot of the quadratic form');
xlabel('x-coordinate'); 
ylabel('y-coordinate');
zlabel('value of the quadratic form');

figure
contour(x,y,z,20); hold on
title('Contour plot of the quadratic form');
xlabel('x-coordinate'); 
ylabel('y-coordinate');

if all(diag(d)>eps),
  u1 = u(:,1)./sqrt(sqrt(d(1,1)));
  u2 = u(:,2)./sqrt(sqrt(d(2,2)));
  U  = [u1 u2];
  quiver(xo,yo,U(1,:),U(2,:),0,'r');
  v = u1'*a*u1;
  fprintf(1,'\nThe value of the quadratic form at the end points: %5.2f\n',v);
  contour(x,y,z,[v v]); hold off
else
  disp('Not possible to draw the vectors, because');
  disp('the matrix A is not positive definite.');
end

clear xy zxy u1 u2
disp('Press any key to continue.');
pause
close all

% ----- Problem 4 ----- 
disp('Problem 4');

% Singular value decomposition (with 'economy size' option, see help)
b = [ 1 2 ; 3 4 ; 5 6 ; 7 8 ]
[u,s,v] = svd(b,0);
disp('Error in the formula:')
disp(b - u*s*v');

disp('Press any key to continue.');
pause

% Eigensystem of b'*b
[w,e] = eig(b'*b);

% Singular values are square roots of the eigenvalues of b'*b.
% Possibly in different order.
disp('Singular values of B and square roots of the eigenvalues of B^T*B:')
disp([diag(s) sqrt(diag(e))]);

disp('Error in matrix representation:')
disp(b'*b - v * s'*s * v');
disp('Diagonal of S^T*S and eigenvalues of B^T*B')
disp([diag(s'*s) diag(e)]);

