%% Iterační metody pro systémy lineárních rovnic
% Příklad: třídiagonální matice

n=10;
D=2*diag(ones(n,1));
E1=-diag(ones(n-1,1),-1);
A=D+E1+E1'

%
% Řešíme systém Ax=b pro b=[1,\dots,1]^T

b=ones(n,1);

% Rozklad matice

L=tril(A,-1);
U=triu(A,1);

%
% Jacobiova metoda

TJ=-inv(D)*(L+U);
gj=inv(D)*b;

% Kolik iterací je potřeba na dosažení přesnosti 10^{-6}?

x=gj; % poc. iterace
x1=TJ*x+gj;
ni=1;
while norm(x-x1)>1e-6
    x=x1;
    x1=TJ*x+gj;
    ni=ni+1;
end
x1
ni

%
% Spektrální poloměr

rj=max(abs(eig(TJ)))

%
% Gaussova-Seidelova metoda

TG=-inv(D+L)*U;
gg=inv(D+L)*b;

% Kolik iterací je potřeba na dosažení přesnosti 10^{-6}?

x=gg; % poc. iterace
x1=TG*x+gg;
ni=1;
while norm(x-x1)>1e-6
    x=x1;
    x1=TG*x+gg;
    ni=ni+1;
end
x1
ni

%
% Spektrální poloměr

rg=max(abs(eig(TG)))

%
% Relaxační metoda

om1=0.8;
TR1=inv(D+om1*L)*((1-om1)*D - om1*U);
g1=om1*inv(D+om1*L)*b;

% Spektrální poloměr

r1=max(abs(eig(TR1)))

%

om2=1.2;
TR2=inv(D+om2*L)*((1-om2)*D - om2*U);
g2=om1*inv(D+om2*L)*b;

% Spektrální poloměr

r2=max(abs(eig(TR2)))

%
% Optimální omega

om_o=2/(1+sqrt(1-rj^2))
TRo=inv(D+om_o*L)*((1-om_o)*D - om_o*U);
go=om_o*inv(D+om_o*L)*b;
ro1=max(abs(eig(TRo)))
ro2=abs(1-om_o)

%
% Kolik iterací je potřeba na dosažení přesnosti 10^{-6}?

x=go; % poc. iterace
x1=TRo*x+go;
ni=1;
while norm(x-x1)>1e-6
    x=x1;
    x1=TRo*x+go;
    ni=ni+1;
end
%x1
ni

%
% Numerický výpočet optimální hodnoty omega

n=10;
D=2*diag(ones(n,1));
E1=-diag(ones(n-1,1),-1);
A=D+E1+E1';
L=tril(A,-1);
U=triu(A,1);
f=@(om) max(abs(eig(inv(D+om*L)*((1-om)*D - om*U))));
options=optimset('Display','iter');
om_opt=fminbnd(f,0,2,options)

%