% Iteracni metoda pro 3-diagonalni matice
%-----------------------------------------%
% Necht je dana soustava Ax=B, kde
% A(N,N) je 3-diagonalni matice:
% A = [a b 0 0 0 .. 0]
%     [c a b 0 0 .. 0]
%     [0 c a b 0 .. 0]
%     ...
%     [0 0 0 ...0 c a]
%
% tedy s prvky na diagonale   = a,
%              pod diagonalou = c,
%              nad diagonalou = b,

% Matice je tedy urcena pouze 4 cisly: N,a,b,c
% Pro iteracni metodu (GS/J) neni tedy potreba
% vzhledem ke strukture matice tuto vubec sestavovat:

N = 500;        % inicializace rozmeru matice (ulohy)
a = -2;
b = 1;
c = 1;

B = ones(N,1);    % vektor prave strany (napr. same 1 ...)
x_J = B;          % pocatecni odhad reseni, napr x = B
x_GS = B;


% pevne dany pocet iteraci 
% lepe -> iterovat dokud neni napr.: norm(x-x_old, 2) < tolerance

pocet_iteraci = 100000; 

for k=1:pocet_iteraci
    
    x_old_J  = x_J; % ulozime k-tou iteraci 
    x_old_GS = x_GS;
    
    % JACOBI---------------------------------
    % 1.radek
    x_J(1) = (1/a)*(B(1) - b*x_old_J(2));
    % 2. - N-1 radek
    for i=2:N-1
        x_J(i) = (1/a)*(B(i) - b*x_old_J(i+1) - c*x_old_J(i-1));
    end
    % N-ty radek
    x_J(N) = (1/a)*(B(N) - c*x_old_J(N-1));

    % GAUSS ---------------------------------
    % 1.radek
    x_GS(1) = (1/a)*(B(1) - b*x_GS(2));
    % 2. - N-1 radek
    for i=2:N-1
        x_GS(i) = (1/a)*(B(i) - b*x_GS(i+1) - c*x_GS(i-1));
    end
    % N-ty radek
    x_GS(N) = (1/a)*(B(N) - c*x_GS(N-1));
    %-----------------------------------------

    err_J(k)  = norm(x_old_J-x_J,  2);
    err_GS(k) = norm(x_old_GS-x_GS,2);

end

iter = 1:pocet_iteraci;
plot(iter, err_GS,'-b',iter,err_J,'-r')
