%% Kilencedik gyakorlat
%% Iterációs megoldási módszerek
%
% Meg akarjuk oldani az
% Ax = b
% egyenletet.
% Alapgondolat: Legyen A=M-N ahol M invertálható, ekkor az előző egyenlet átalakítva:
% (M-N)x = b\\
% Mx - Nx = b\\
% Mx = Nx + b\\
% x = M^{-1}N x + M^{-1} b
% Tehát keressük az f(x)=M^{-1} N x + M^{-1} b függvény egy fixpontját (azaz ahol f(x)=x), ezt az alábbi iteráció segítségével tehetjük meg:
% x_{k+1}= M^{-1}N x_k + M^{-1}b
% Ha használjuk a B = M^{-1}N és f = M^{-1} b jelöléseket, akkor az előző az alábbi iterációs képlettel írható le:
% x_{k+1} = B x_k + f,
% azaz minden lépésben megszorozzuk a megoldásokat tartalmazó x_k vektort egy B mátrixszal, majd hozzáadjuk az f vektort.
% Tétel: Az iteráció konvergál, ha \rho(B)<1(azaz |\lambda_i|<1 a B mátrix minden \lambda_i sajátértékére).
% Módszerek:
% Legyen L: főátló alatti elemeket, D a főátló elemeit, illetve R a főátló feletti elemek mátrixsza. Ekkor négy módszert különböztetünk meg:
% Jacobi iteráció: M = D, \quad N=-L-R.
% Gauss-Seidel iteráció: M = D + L, \quad N=-R.
% Relaxált Jacobi (JOR): M = \frac1\omega D, \quad N=\frac{1-\omega}{\omega}D - L - R.
% Relaxált Gauss-Seidel (SOR): M = \frac1\omega D + L, \quad N=\frac{1-\omega}{\omega}D - R.
% Tétel: Ha A szigorúan diagonálisan domináns (azaz a főátlóban lévők absz. értékei nagyobbak az adott sorban lévő abszolútérték-összegénél), akkor a Jacobi és a Gauss-Seidel iteráció konvergens.
% Definíció: Egy négyzetes mátrixot M-mátrixnak nevezünk, ha főátlóján kívül nempozitív elemei vannak, és létezik egy olyan g>0 vektor, amivel Ag >0.
% Tétel: Ha A M-mátrix, akkor konvergens a JOR és az SOR \omega \in (0,1) esetén.
% Tétel: Ha A szimmetrikus és pozitív definit, akkor konvergens az SOR iteráció \omega \in (0,2) esetén.
%
% 1. feladat: Tekintsük az alábbi egyenletrendszert:
% \left[\matrix{2 & -1 \cr -1 & 2}\right] \overline{x} = \left[ \matrix{1 \cr 3} \right].
% Melyik fenti módszerrel oldható meg?
% Megoldás: Látható, hogy diagonálisan domináns, így a Jacobi és a Gauss-Seidel iteráció konvergens lesz az esetében.
% Látható, hogy a főátlón kívüli elemek nempozitívak, továbbá ha g = \left(\matrix{1 \cr 1}\right)>0, akkor Ag = \left(\matrix{1 \cr 1}\right)>0, így M-mátrix, azaz működik a JOR és az SOR is \omega \in (0,1) esetén.
% A mátrix szimmetrikus és mivel a_{11}=2>0 és det(A)=3>0, így pozitív definit is, azaz az SOR iteráció is működni fog \omega \in (0,2) esetén.
%
% 2. feladat: (a) Alkalmazzuk a fenti egyenletrendszerre a JOR iterációt! Írjuk fel az iteráció lépéseit, majd programozzuk is le x_0 = ( \matrix{0 &  0})^T  kezdeti vektorral!
% Megoldás: Esetünkben D = \left(\matrix{2 & 0 \cr 0 & 2}\right),  L = \left(\matrix{0 & 0 \cr -1 & 0}\right),  R = \left(\matrix{0 & -1 \cr 0 & 0}\right),  így M = \frac1\omega \left(\matrix{2 & 0 \cr 0 & 2}\right) és
%  N = \frac{1-\omega}{\omega}\left(\matrix{2 & 0 \cr 0 & 2}\right) + \left(\matrix{0 & 1 \cr 1 & 0}\right)= \left(\matrix{2 \frac{1-\omega}{\omega} & 1 \cr 1 & 2 \frac{1-\omega}{\omega}}\right).
% Ezek alapján a B mátrix:
% B = M^{-1} N = \left( \matrix{2/\omega & 0 \cr 0 & 2/\omega}\right)^{-1} \left(\matrix{2 \frac{1-\omega}{\omega} & 1 \cr 1 & 2 \frac{1-\omega}{\omega}}\right) =
% = \left( \matrix{\omega/2 & 0 \cr 0 & \omega/2}\right) \left(\matrix{2 \frac{1-\omega}{\omega} & 1 \cr 1 & 2 \frac{1-\omega}{\omega}}\right) = \left( \matrix{1-\omega & \omega/2 \cr \omega/2 & 1-\omega}\right)
% Magához az iterációhoz még fel kell írnunk az f vektort:
% f = M^{-1} b = \left( \matrix{2/\omega & 0 \cr 0 & 2/\omega}\right)^{-1} \left( \matrix{1 \cr 3}\right) = \left( \matrix{\omega/2 & 0 \cr 0 & \omega/2}\right) \left( \matrix{1 \cr 3}\right) = \left( \matrix{\omega/2 \cr 3 \omega/2}\right)
% Most már elkezdhetjük a kód megírását!

% Az omega paraméter értéke:
omega=0.5;

% Kezdeti vektor
x=zeros(2,1)

% A B mátrix:
B = [1-omega omega/2; omega/2 1-omega]

% Az f vektor:
f = [omega/2 3*omega/2]'

% Lépések száma
k=1000;

for i=1:k
    x = B*x+f;
end

x

%
% (b) Teszteljük, hogy mely értékek esetén milyen gyors lesz a konvergencia! Ezt úgy tegyük meg, hogy tekintsük pontosnak a megoldásunkat, ha
% \left\Vert x_k - \left( \matrix{5/3 \cr 7/3} \right)\right\Vert_{\infty}<10^{-6}.
% Mivel lehet magyarázni ezen jelenséget?
% Megoldás: Írjunk kódot, ami adott \omega értékre visszaadja a lépések számát!

omega=linspace(0.1,1,10)
k=zeros(10,1)

for i=1:10
    k(i)=lepesszam(omega(i));
end

plot(omega,k)

% A jelenség magyarázatához vizsgáljuk meg a B = \left( \matrix{1-\omega & \omega/2 \cr \omega/2 & 1-\omega}\right) mátrix sajátértékeit! Ezek \lambda_1 = 1-\omega/2 és \lambda_2 = 1-3\omega/2. Mivel \omega\in (0,1), és itt \left| 1 - \omega/2\right| > \left| 1 - 3\omega/2\right|, ezért \rho(B)=\left| 1 - \omega/2\right|, amelynek az értéke monoton csökken, amint \omega tart az 1 felé, emiatt az 1 körüli \omega értékekre lesz a leggyorsabb a konvergencia (általánosságban, minél kisebb \rho(B), annál gyorsabb a folyamat.
%
% 3. feladat: Oldjuk meg az 1. feladat egyenletét Gauss-Seidel iteráció segítségével is! Írjunk programot, amely elvégzi az iterációt!
% Megoldás: Felhasználva a 2. feladatban kiszámított mátrixokat, M = \left( \matrix{2 & 0 \cr -1 & 2} \right) és N = \left( \matrix{0 & 1 \cr 0 & 0} \right), azaz
% B = M^{-1} N = \left( \matrix{2 & 0 \cr -1 & 2} \right)^{-1} \left( \matrix{0 & 1 \cr 0 & 0} \right) =  \left( \matrix{1/2 & 0 \cr 1/4 & 1/2} \right) \left( \matrix{0 & 1 \cr 0 & 0} \right) = \left( \matrix{0 & 1/2 \cr 0 & 1/4} \right)
% és
% f = M^{-1} b = \left( \matrix{1/2 & 0 \cr 1/4 & 1/2} \right)  \left( \matrix{1  \cr 3 } \right) = \left( \matrix{1/2  \cr 7/4 } \right) .
% Azaz így az iteráció alakja:
% x_{k+1} = \left( \matrix{0 & 1/2 \cr 0 & 1/4} \right) x_k + \left( \matrix{1/2  \cr 7/4 } \right) .
% Az iteráció kódja:

% Kezdeti vektor
x=zeros(2,1)

% A B mátrix:
B = [0 1/2; 0 1/4]

% Az f vektor:
f = [1/2 7/4]'

% Lépések száma
k=1000;

for i=1:k
    x = B*x+f;
end

x

%

function [k]=lepesszam(omega)
x=zeros(2,1);

B = [1-omega omega/2; omega/2 1-omega];

f = [omega/2 3*omega/2]';

pontos = [5/3 7/3]';

k=0;
while norm(x-pontos,'Inf')>10^(-6)
    x = B*x+f;
    k=k+1;
end

end