% template.m
% This is a template for electronic submission of homework
% assignments for the course MATH 3Q03; use of this template is 
% obligatory
% Email the complete file to the instructor before the
% due date; before submission make sure that the code executes
% correctly in MATLAB

% Written by Bartosz Protas, September 2016

function hw2_sol

close all, clear all

% Student information
  
disp('  ------------   Student Information  ------------  ')  

% Please enter your information here
Name = 'Johny';
Surname = 'Good';
ID = 123456789;

fprintf('  Student:  %s %s (ID: %d) \n', Name, Surname, ID);



disp('  ------------------------------------------------  ')  
disp('  ')
disp('  ')
disp('  ')
disp('  ')


disp('  ----------   Solution to Question #1  ----------  ')  
disp('  ')
% Insert here your solution to Question #1
% (if you define your own functions, they may appear 
%  at the end of the file)

% Insert here your solution to Question #1
% (if you define your own functions, they may appear 
%  at the end of the file)

% setting up the matrix A;
A = [6  1  0  3; 5  8  4  0;  2  2 7 1;  0  1  2  4];

% setting up the RHS vectors
y1 = [ 4 1 9 2 ]';

for i=1:4
  y2(i) = y1(i)^2;
  y3(i) = y1(i)^3;
  y4(i) = 1.0 / y1(i);
  y5(i) = sqrt(y1(i));
end

% part (a) - performing the LU decompostion;
[L, U] = lu(A);
disp('         Matrix L:')
L
disp('         Matrix U:')
U

% part (b) solving the system in two steps;

z1 = Lbck(L,y1);   x1 = Ubck(U, z1);
z2 = Lbck(L,y2);   x2 = Ubck(U, z2);
z3 = Lbck(L,y3);   x3 = Ubck(U, z3);
z4 = Lbck(L,y4);   x4 = Ubck(U, z4);
z5 = Lbck(L,y5);   x5 = Ubck(U, z5);

% part (c)
disp('    '); disp('    '); 
disp('First entries x1(1),...,x5(1) of the solution vectors:')
[x1(1) x2(1) x3(1) x4(1) x5(1)]





disp('  ')
disp('  ------------------------------------------------  ')  
disp('  ')
disp('  ')
disp('  ')
disp('  ')


disp('  ----------   Solution to Question #2  ----------  ')  
disp('  ')
% Insert here your solution to Question #2
% (if you define your own functions, they may appear 
%  at the end of the file)


i = 1;
for n=25:25:500
% part (a)  
  An = hilb( n );    
% part (b)  
  [N1, N2, Ninf, Nf] = Matrix_Norms( An, n );
% part (c)  
  buf = sprintf('n=%3d  N1=%8g  N2=%8g  Ninf=%8g  Nf=%8g', n, N1, N2, ...
		Ninf, Nf);
  disp(buf);
  Data(i,1) = n;
  Data(i,2) = N1;
  Data(i,3) = N2;
  Data(i,4) = Ninf;
  Data(i,5) = Nf;
  
  i = i + 1;
  
end;
  
% plotting the results;
figure(1);  plot(Data(:,1),Data(:,2),'x-',...
		 Data(:,1),Data(:,3),'o-',...
		 Data(:,1),Data(:,4),'s-',...
		 Data(:,1),Data(:,5),'d-');
legend('1-norm','2-norm','inf-norm','Frobenius norm');
    



disp('  ')
disp('  ------------------------------------------------  ')  
disp('  ')
disp('  ')
disp('  ')
disp('  ')


disp('  ----------   Solution to Question #3  ----------  ')  
disp('  ')
% Insert here your solution to Question #3
% (if you define your own functions, they may appear 
%  at the end of the file)


N = 10;  
M = 21;
TOL = 1.0e-6;

W = linspace(0.0,2.0,M);

max_iter = 1000; % maximum number of iterations 
m = 1:1:max_iter;

A = gallery('poisson', N); % FD approximation of the Laplacian
f = zeros(N*N, 1); % right hand side
x0 = ones(N*N, 1); % starting guess


for i=1:M
  [it(i), err_gs] = sor(A, f, x0, max_iter, TOL, W(i));
end;

figure(2)
plot(W,it,'o-');
xlabel('\omega'); ylabel('number of iterations');

  


disp('  ')
disp('  ------------------------------------------------  ')  
disp('  ')
disp('  ')
disp('  ')
disp('  ')


% Please copy as many headers as necessary

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x] = Lbck( L, y )
% Parameters: 
%  L - lower triangular matrix
%  y - RHS vector
  
  n = length(y);
  for i=1:n;
    x0(i) = y(i) / L(i,i);
    for j=1:i-1
      x0(i) = x0(i) - L(i,j) * x0(j);      
    end
  end

  x = x0';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x] = Ubck( U, y )
% Parameters: 
%  U - upper triangular matrix
%  y - RHS vector
  
  n = length(y);
  for i=n:-1:1;
    x0(i) = y(i) / U(i,i);
    for j=n:-1:i+1
      x0(i) = x0(i) - U(i,j) * x0(j) / U(i,i);      
    end
  end

  x = x0';

  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [n1,n2,ninf,nf] = Matrix_Norms( A,  N )
% Calculates various matrix norms

% calculating the "1" norm;
  n1 = 0.0;
  for j=1:N
    sum = 0.0;
    for i=1:N
      sum = sum + abs( A(i,j) );
    end;
    if ( sum > n1 )
      n1 = sum;
    end;
  end;
  
% calculating the "2" norm
  n2 = sqrt( max( abs( eig( A'*A ) ) ) );
    
% calculating the "1" norm;
  ninf = 0.0;
  for i=1:N
    sum = 0.0;
    for j=1:N
      sum = sum + abs( A(i,j) );
    end;
    if ( sum > ninf )
      ninf = sum;
    end;
  end;
  
% calculating the Frobenius norm;
  nf = 0.0;
  for i=1:N
    for j=1:N
      nf = nf + A(i,j)^2;
    end;
  end;
  nf = sqrt(nf);  

  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% w - over-relaxation paramater;
function [m,er] = sor(A, f, x, max_iter, TOL, w)

er = [];
n = size(A,1);
D = diag(diag(A));
L = tril(A,-1);
U = triu(A,1);
normf = norm(f);

m = max_iter;
i = 1;
while (i <= max_iter)
  x = (D+w*L)\(-(w*U + (w-1)*D)*x + w*f);
  if (normf)
    er = [er; norm(A*x - f)/normf];
  else
    er = [er; norm(x)];
  end
  if (er(i) < TOL)
    m = i;
    i = max_iter;
  end
  i = i + 1;
end
