%% Metoda Jacobi pentru rezolvarea sistemelor de ecuatii liniare
% Date de intrare: 
    %Matricea coeficientilor
    %Vectorul termenilor liberi
    %Aproximatia initiala
    %Toleranta de calcul
% Date de iesire: 
    %Solutia sistemului
    %Numarul necesar de iteratii
    %Proba solutiei
    %Suplimetar se pot vizualiza iteratiile in variabila Iteratii in care 
        % pe prima coloana este numarul iteratiei iar pe urmatoarele n
        % coloane sunt valorile necunoscutelor in iteratia curenta.

clearvars
clc

%% Date de intrare

% Matricea coeficientilor A
A= [2 1 -1
    1 2 1
    -1 1 -10]; 

% Vectorul termenilor liberi
B = [4 4 4]'; 
n=size(B,1);

% Toleranta
tol=1e-6;

%Aproximatia initiala
%Varianta 1: zerouri
%x=zeros(n,1); 

%Varianta 2: vectorul g (b(i)/a(i,i))
for i=1:n
    g(i)=B(i)/A(i,i);
end
x=g';


%% Programul principal
for i=1:n
  suma_rand=0;
  for j=1:n
    suma_rand=suma_rand+abs(A(i,j));
  end
  if suma_rand-abs(A(i,i)) > abs(A(i,i))  % verificam daca matricea este diagonal dominanta
    disp('OBSERVATIE: matricea nu este diagonal dominanta');
    break;
    %disp(suma_rand-abs(A(i,i)));
    %disp(A(i,i));
  end
 end
  
normVal=Inf; 
itr=0;

while normVal>tol
    xold=x;
    
    for i=1:n
      sigma=0;
        for j=1:n
            if j~=i
                sigma=sigma+A(i,j)*x(j);
            end
        end
      x(i)=(1/A(i,i))*(B(i)-sigma);
    end
    itr=itr+1;
    normVal=norm(xold-x);
    Iteratii(itr,1)=itr;
        for k=1:n 
         Iteratii(itr,k+1)=x(k);
        end
end

%% Afisare rezultate

for j=1:n 
    if x(j)==inf
      disp ('   ');
      disp('Solutia nu s-a putut calcula');
      fprintf ( 'Numarul curent de iteratii este %1.0f  \n',itr);
    return
    end
    
end

disp('Solutia sistemului de ecuatii este:')
for k=1:n 
fprintf ( '   x(%1.0f) = %2.10f  \n',k,x(k))
end

disp ('   ');
fprintf ( 'Solutia s-a determinat dupa %1.0f iteratii  \n',itr);

%proba solutiei; practic se face operatia A*X-B
disp ('   ');
disp('Proba solutiei A*X-B=0:') 

for k1=1:n 
fprintf ( '  Ecuatia %1.0f: %2.20f  \n', k1, A(k1,:)*x-B(k1));
end
