%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% jacobi.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,iter,err] = jacobi(A,x,b,iter_max,tol)
%
% [x,iter,err] = jacobi(A,x,b,iter_max,tol);
%
% A           matrix
% x           x_0
% b           right hand side
% iter_max    maximum number of iteration
% tol         tolerance
%
% x           solution
% iter        number of iteration
% err         behaviour of the error during the convergece

iter = 0;

bn2 = norm(b);
if bn2==0.0
  bn2 = 1.0;
end

r = b-A*x;
err(1) = norm(r)/bn2;
if err(1)<=tol
  return;
end

E = tril(A,-1);    % estrae la parte triangolare inferiore da A
D = diag(diag(A)); % costruisce la matrice diagonale prendendo
                   % la diagonale principale di A
F = triu(A,1);     % estrae la parte triangolare superiore da A

M = D;
N = -E-F;

for iter = 1:iter_max
   x_old = x;
   x = M\(N*x+b); % matlab si accorge automaticament che la matrice M
                  % e` diagonale ed utilizza il metodo
                  % piu` opportuno per risolvere il sistema lineare
   err(iter) = norm(x-x_old)/norm(x);
   if err(iter)<=tol
     break;
   end
end

if iter>=iter_max
  disp('ERRORE: numero eccessivo di iterazioni in Jacobi\n')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gauss_seidel.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,iter,err] = gauss_seidel(A,x,b,iter_max,tol)
%
% [x,iter,err] = gauss_seidel(A,x,b,iter_max,tol);
%
% A           matrix
% x           x_0
% b           right hand side
% iter_max    maximum number of iteration
% tol         tolerance
%
% x           solution
% iter        number of iteration
% err         behaviour of the error during the convergece


iter = 0;

bn2 = norm(b);
if bn2==0.0
  bn2 = 1.0;
end

r = b-A*x;
err(1) = norm(r)/bn2;
if err(1)<=tol
  return;
end

E = tril(A,-1);    % estrae la parte triangolare inferiore da A
D = diag(diag(A)); % costruisce la matrice diagonale prendendo
                   % la diagonale principale di A
F = triu(A,1);     % estrae la parte triangolare superiore da A

M = E+D;
N = -F;

for iter = 1:iter_max
   x_old = x;
   x = M\(N*x+b); % matlab riconosce automaticamente che la matrice M
                  % e` traingolare inferiore ed utilizza il metodo
                  % piu` opportuno per "risolvere" il sistema lineare
   err(iter) = norm(x-x_old)/norm(x);
   if err(iter)<=tol
     break;
   end
end

if iter>=iter_max
  disp('ERRORE: numero eccessivo di iterazioni in Gauss-Seidel\n')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% comp_J_GS.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Confronto metodi di Jacobi e Gauss-Seidel

echo off
disp(sprintf('\n\n1) Jacobi non converge, Gauss-Seidel converge\n'))
pause
A = [3,0,4;7,4,2;-1,-1,-2];
b = [7;13;-4];

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = D;
N = -E-F;
B_J = inv(M)*N;
rho_J = max(abs(eig(B_J)));
disp(sprintf('rho_J=%f\n',rho_J))

M = E+D;
N = -F;
B_GS = inv(M)*N;
rho_GS = max(abs(eig(B_GS)));
disp(sprintf('rho_GS=%f\n',rho_GS))

x = zeros(3,1);
[x,iter_J,err_J] = jacobi(A,x,b,1000,1E-5);

x = zeros(3,1);
[x,iter_GS,err_GS] = gauss_seidel(A,x,b,1000,1E-5);

figure
hold off
semilogy(1:iter_J,err_J,'r-')
hold on
semilogy(1:iter_J,err_J,'ro')
semilogy(1:iter_GS,err_GS,'g-')
semilogy(1:iter_GS,err_GS,'g*')
grid on
axis([0, min(iter_J,iter_GS), 0, max(max(err_J),max(err_GS))])

%%%%%%%%%%

disp(sprintf('\n\n2) Jacobi converge, Gauss-Seidel non converge\n'))
pause
A = [-3,3,-6;-4,7,-8;5,7,-9];
b = [-6;-5;3];

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = D;
N = -E-F;
B_J = inv(M)*N;
rho_J = max(abs(eig(B_J)));
disp(sprintf('rho_J=%f\n',rho_J))

M = E+D;
N = -F;
B_GS = inv(M)*N;
rho_GS = max(abs(eig(B_GS)));
disp(sprintf('rho_GS=%f\n',rho_GS))

x = zeros(3,1);
[x,iter_J,err_J] = jacobi(A,x,b,1000,1E-5);

x = zeros(3,1);
[x,iter_GS,err_GS] = gauss_seidel(A,x,b,1000,1E-5);

hold off
figure
semilogy(1:iter_J,err_J,'r-')
hold on
semilogy(1:iter_J,err_J,'ro')
semilogy(1:iter_GS,err_GS,'g-')
semilogy(1:iter_GS,err_GS,'g*')
grid on
axis([0, min(iter_J,iter_GS), 0, max(max(err_J),max(err_GS))])

%%%%%%%%%%

disp(sprintf('\n\n3) Jacobi converge piu lentamente di Gauss-Seidel\n'))
pause
A = [4,1,1;2,-9,0;0,-8,-6];
b = [6;-7;-14];

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = D;
N = -E-F;
B_J = inv(M)*N;
rho_J = max(abs(eig(B_J)));
disp(sprintf('rho_J=%f\n',rho_J))

M = E+D;
N = -F;
B_GS = inv(M)*N;
rho_GS = max(abs(eig(B_GS)));
disp(sprintf('rho_GS=%f\n',rho_GS))

x = zeros(3,1);
[x,iter_J,err_J] = jacobi(A,x,b,1000,1E-5);

x = zeros(3,1);
[x,iter_GS,err_GS] = gauss_seidel(A,x,b,1000,1E-5);

hold off
figure
semilogy(1:iter_J,err_J,'r-')
hold on
semilogy(1:iter_J,err_J,'ro')
semilogy(1:iter_GS,err_GS,'g-')
semilogy(1:iter_GS,err_GS,'g*')
grid on
axis([0, max(iter_J,iter_GS), 0, max(max(err_J),max(err_GS))])

%%%%%%%%%%

disp(sprintf('\n\n4) Jacobi converge piu velocemente di Gauss-Seidel\n'))
pause
A = [7,6,9;4,5,-4;-7,-3,8];
b = [22;5;-2];

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = D;
N = -E-F;
B_J = inv(M)*N;
rho_J = max(abs(eig(B_J)));
disp(sprintf('rho_J=%f\n',rho_J))

M = E+D;
N = -F;
B_GS = inv(M)*N;
rho_GS = max(abs(eig(B_GS)));
disp(sprintf('rho_GS=%f\n',rho_GS))

x = zeros(3,1);
[x,iter_J,err_J] = jacobi(A,x,b,1000,1E-5);

x = zeros(3,1);
[x,iter_GS,err_GS] = gauss_seidel(A,x,b,1000,1E-5);

hold off
figure
semilogy(1:iter_J,err_J,'r-')
hold on
semilogy(1:iter_J,err_J,'ro')
semilogy(1:iter_GS,err_GS,'g-')
semilogy(1:iter_GS,err_GS,'g*')
grid on
axis([0, max(iter_J,iter_GS), 0, max(max(err_J),max(err_GS))])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conv_J.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Velocita` di convergenza del metodo di Jacobi

echo off

A = [4 1 0 ; 2 5 1 ; -1 2 4];
b = [1;0;3];

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = D;
N = -E-F;
B_J = inv(M)*N;
rho_J = max(abs(eig(B_J)));
disp(sprintf('rho_J_1=%f\n',rho_J))

x = zeros(3,1);
[x,iter_J_1,err_J_1] = jacobi(A,x,b,1000,1E-5);

%%%

A = [4 1 0 ; 2 5 1 ; -1 2 4] + 2*eye(3);

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = D;
N = -E-F;
B_J = inv(M)*N;
rho_J = max(abs(eig(B_J)));
disp(sprintf('rho_J_2=%f\n',rho_J))

x = zeros(3,1);
[x,iter_J_2,err_J_2] = jacobi(A,x,b,1000,1E-5);

%%%

A = [4 1 0 ; 2 5 1 ; -1 2 4] + 20*eye(3);

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = D;
N = -E-F;
B_J = inv(M)*N;
rho_J = max(abs(eig(B_J)));
disp(sprintf('rho_J_3=%f\n',rho_J))

x = zeros(3,1);
[x,iter_J_3,err_J_3] = jacobi(A,x,b,1000,1E-5);


figure
hold off
semilogy(1:iter_J_1,err_J_1,'ro-')
hold on
semilogy(1:iter_J_2,err_J_2,'go-')
semilogy(1:iter_J_3,err_J_3,'bo-')

grid on
axis([0, min(iter_J_1,iter_J_2), 0, max(max(err_J_1),max(err_J_2))])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conv_GS.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrici non a dominanza diagonale, Gauss-Seidel

echo off
disp(sprintf('1) Matrice non a dominanza diagonale per cui Gauss-Seidel non converge\n'))
pause
A = [3 -5 47 20; 11 16 17 10; 56 22 11 -18; 17 66 -12 7];
b = [18;26;34;82];

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = E+D;
N = -F;
B_GS = inv(M)*N;
rho_GS = max(abs(eig(B_GS)));
disp(sprintf('rho_GS=%f\n',rho_GS))

x = ones(4,1);
[x,iter_GS,err_GS] = gauss_seidel(A,x,b,1000,1E-5);


figure
semilogy(1:iter_GS,err_GS,'g-')
grid on

%%%%%%%%%%

disp(sprintf('2) Matrice non a dominanza diagonale (ultima riga) per cui Gauss-Seidel converge\n'))
pause
A = [56 22 11 -18; 17 66 -12 7; 3 -5 47 20; 11 16 17 10];
b = [34;82;18;26];

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = E+D;
N = -F;
B_GS = inv(M)*N;
rho_GS = max(abs(eig(B_GS)));
disp(sprintf('rho_GS=%f\n',rho_GS))

x = ones(4,1);
[x,iter_GS,err_GS] = gauss_seidel(A,x,b,1000,1E-5);


figure
semilogy(1:iter_GS,err_GS,'g-')
grid on

%%%%%%%%%%

disp(sprintf('3) Matrice non a dominanza diagonale (ultima riga) per cui Gauss-Seidel converge\n'))
pause
A = [56 22 11 -18; 17 66 -12 7; 3 -5 47 20; 11 16 17 13];
b = [34;82;18;26];

E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = E+D;
N = -F;
B_GS = inv(M)*N;
rho_GS = max(abs(eig(B_GS)));
disp(sprintf('rho_GS=%f\n',rho_GS))

x = ones(4,1);
[x,iter_GS,err_GS] = gauss_seidel(A,x,b,1000,1E-5);


figure
semilogy(1:iter_GS,err_GS,'g-')
grid on


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% iter_trid.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Confronta i diversi metodi iterativi per matrici tridiagonali

echo off
sprintf('Matrici tridiagonali\n')
pause

n = 10;

d  = ones(n,1);
lu = ones(n-1,1);
A  = -diag(lu,-1)+2*diag(d)-diag(lu,1);
b  = A*d;



E = tril(A,-1);
D = diag(diag(A));
F = triu(A,1);

M = D;
N = -E-F;
B_J = inv(M)*N;
rho_J = max(abs(eig(B_J)));
sprintf('rho_J=%f\n',rho_J)

M = E+D;
N = -F;
B_GS = inv(M)*N;
rho_GS = max(abs(eig(B_GS)));
sprintf('rho_GS=%f\n',rho_GS)
sprintf('rho_J^2=%f\n',rho_J^2)

x = zeros(n,1);
[x,iter_J,err_J] = jacobi(A,x,b,1000,1E-5);

x = zeros(n,1);
[x,iter_GS,err_GS] = gauss_seidel(A,x,b,1000,1E-5);


figure
hold off
semilogy(1:iter_J,err_J,'ro-')
hold on
semilogy(1:iter_GS,err_GS,'g*-')
grid on
axis([0, max([iter_J,iter_GS]),...
       0, max([max(err_J),max(err_GS)])])