%% Iterativ megoldók az Ax = b egyenletrendszerre (Jacobi, Gauss-Seidel; relaxáltjaik: JOR, SOR)
% Meg akarjuk oldani az
%  Ax = b
% egyenletet, azonban A mérete túl nagy ahhoz, hogy a hagyományos, n^3-ös direkt módszereket használhassuk.
% Motiváció: parcdiff



% Tétel (Banach fixpont): Ha X teljes metrikus tér, f: X \to X
% kontrakció, akkor f-nek pontosan egy fixpontja van, továbbá tetszőleges x_0 esetén az x_n = f(x_{n-1})
%  képlettel definiált sorozat tart ehhez a fixponthoz.
% Következmény: Ha f: \mathbb{R}^n \to \mathbb{R}^n egy kontrakció az adott normában, azaz Lipschitz-folytonos és a megfelelő 0 \leq L < 1 konstans mellett minden x,y pontra \| f(x) - f(y) \| \leq L \| x -y\| teljesül, akkor f-nek van fixpontja és tetszőleges kezdőpontból indulva, az f ismételt alkalmazásával tartani tudunk ehhez.
% Példa: f(x) = \frac{1}{2}x melynek fixpontja a 0 és valóban bármely x_0-ból indulva, f-et n-szer alkalmazva a 2^{-n} x_0 érték a nullához fog tartani.
% Példa: Véges dimenzióban az affin függvények f(x) = Ax +b alakúak valamely A mátrixra és b vektorra. Ekkor f(x) - f(y) = Ax - Ay = A(x - y), továbbá adott norma és az indukált mátrixnorma esetén \| f(x) - f(y) \| \leq \| A \| \| x -y \|, tehát adott normában ez az affin függvény kontrakció, ha a megfelelő A mátrix indukált normája kisebb, mint egy.
% Tétel: Véges dimenzióban a normák ekvivalensek.
% Következmény: Véges dimenzióban ha egy sorozat konvergens, akkor minden normában konvergens; ha van olyan norma, amiben nem konvergens, akkor a többi normában sem lehet konvergens.
% Következmény: Egy affin függvény ismételt alkalmazásával kapott iteráció konvergens, ha van olyan indukált norma, ami szerinti normája a függvény lineáris részét leiró A mátrixnak kisebb, mint 1.
% Állitás: Jelölje \rho(A) az A mátrix legnagyobb abszolútértékű sajátértékének abszolútértékét, azaz legyen \rho a spektrálsugár. Ekkor tetszőleges indukált mátrixnormára és k\geq1egész számra \rho(A)^k \leq \| A^k\|, hiszen ha \lambda sajátérték v normált sajátvektorral, akkor \|A^k \| = \sup_{\|u\| = 1} \| A^ku \| \geq \| A^k v\| = |\lambda^k| \| v\| = |\lambda^k|.
% Tétel: \rho(A) = \inf\{\| A\| : \| \cdot \| \text{ indukált mátrixnorma}\}
% Következmény: Indukált mátrixnormák segitségével tetszőlegesen meg lehet közeliteni a spektrálsugár értékét.
% Tétel: Egy affin függvény ismételt alkalmazásával kapott iteráció konvergens, ha a függvény lineáris részét leiró A mátrix spektrálsugara kisebb, mint 1.
% Biz.: Ha ez kisebb mint egy, akkor valamely indukált normája is kisebb mint egy, tehát kontrakció, tehát konvergens.
%
%% Hogyan készitünk egy ilyen fixpont-iterációt I.
% Legegyszerűbb megközelités (egyszerű- vagy Richardson-iteráció)
% Ax = b\\
% 0 = b - Ax \\
% x = x - Ax + b \\
% \\
% f(x) = (I-A)x +b
% Egy gond ezzel, hogy sokszor az (I-A) mátrix spektrálsugara még nem elég kicsi. Ezen segithetünk egy \omega paraméter bevezetésével:
% Ax = b\\
% 0 = \omega(b - Ax) \\
% x = x - \omega Ax + \omega b \\
% \\
% f_{\omega}(x) = (I-\omega A)x +\omega b
% Itt \rho(I-\omega A) < 1 pontosan akkor teljesül, ha A \lambda sajátértékeire | 1 - \omega \lambda | < 1. A konvergencia akkor a leggyorsabb, ha ez a spektrálsugár minél kisebb. Például ha az A mátrix szimmetrikus és pozitiv definit akkor az optimális választás \omega-ra: \frac{2}{\lambda_{\min} + \lambda_{\max}}.
% 0. Feladat: Miért ez az optimális választás?

omegas = linspace(0, 3, 100);
lambdas = [0.5, .6, .75, 1.5]

for lambda=lambdas
    plot(omegas, abs(1-lambda*omegas));
    hold on;
end

legend(".5", ".6", ".75", "1.5")
xlabel("\omega")

hold off;

%% Hogyan készitünk egy ilyen fixpont-iterációt II.
% A fentihez hasonló átalakitás, ha A = M - N felbontással élünk, majd ezzel számolva
% Ax = b\\
% (M-N)x = b\\
% Mx = Nx + b\\
% x = M^{-1}N x + M^{-1} b\\
% \\
% f(x) = M^{-1}N x + M^{-1} b
% Itt tehát az iterációs mátrix B = M^{-1}N és ennek a spektrálsugarát már hatékonyabban tudjuk befolyásolni az M és N alkalmas megválasztásával.
% Módszerek:
% Legyen A = L + D + Uegy felbontása az A mátrixnak rendre szigorú alsó-háromszög, diagonális, és szigorú felső-háromszög mátrixokra. Ekkor tekintsük az alábbi nevezetes módszerek formuláit:
% Jacobi iteráció: M = D, \quad N=-L-U
% .
% Gauss-Seidel iteráció: M = D + L, \quad N=-U.
% Relaxált Jacobi (JOR): M = \frac1\omega D, \quad N=\frac{1-\omega}{\omega}D - L - U.
% Relaxált Gauss-Seidel (SOR): M = \frac1\omega D + L, \quad N=\frac{1-\omega}{\omega}D - U.

%(L+D+U)x = b
%Dx = b - (L+U)x
%
% x = inv(D)*b - inv(D)*(L+U)x

% 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.
% Észrevehetjük, hogy ebben a példában egy ismerős tridiagonális mátrixcsalád egy eleméről van szó. Ezek ritka mátrixok (tehát nagyon sok bennük a 0 elem), igy a velük való számolás lényegesen gyorsabb tud lenni.

n = 50;
e = ones(n,1);
A = spdiags([-e 2*e -e], [-1 0 1], n, n);
spy(A)
full(A)

% láthatjuk, hogy a sajátértékek mind pozitivak
all(eig(A) > 0)
% a mátrix persze szimmetrikus is

%
% 2. feladat: (a) Alkalmazzuk a fenti egyenletrendszerre a JOR iterációt! Írjuk fel az iterációhoz használható mátrixokat és vektorokat, 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),  U = \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 az iterációs 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)
% Az affin eltolást megadó vektor:
% 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 eltolás vektora:
r = [omega/2;
     3*omega/2];

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

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

x

[2, -1;
-1,  2]*x - [1;3]

%
% (b) Milyen \omega választással gyorsitható a konvergencia?
% A B = \left( \matrix{1-\omega & \omega/2 \cr \omega/2 & 1-\omega}\right) mátrix sajátértékei \lambda_1 = 1-\omega/2 és \lambda_2 = 1-3\omega/2. Mivel \omega\in (0,1), itt

omegas = linspace(0, 2, 100);
plot(omegas, abs(1-omegas/2), omegas, abs(1-1.5*omegas))
legend(".5", "1.5")

% azaz \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.
% Konklúzió: Itt nem érdemes relaxálni.
%
% 3. feladat: Irjunk általános függvényt a fenti A = M -N felbontással kapott iterációkhoz. Emlék: f(x) = M^{-1}N x + M^{-1} b.

A = [2, -1;
    -1,  2]

b = [1;3];

[x, diff_a, diff_r, step_num] = iteracio(diag(diag(A)), diag(diag(A)) - A, b, zeros(2,1), 1e-3, 1e-4, 1000)

% mi lehet a spektrálsugár ennél a választásnál?

eig(inv(diag(diag(A)))*(diag(diag(A)) - A))

%

tril(magic(4), -1)

%

[x, diff_a, diff_r, step_num] = jacobi_iteracio(A, b, zeros([2,1]), 1e-3, 1e-3, 100)
[x, diff_a, diff_r, step_num] = sor_iteracio(A, b, zeros([2,1]), 0.5, 1e-3, 1e-3, 100)

%

function [x, diff_a, diff_r, step_num] = iteracio(M, N, b, x0, atol, rtol, max_steps)
    r = M\b;
    % erdemes-e kihozni ezt ide vagy sem (?)
    B = M \ N;
    
    x = x0;
    for step_num=1:max_steps
        % ha x_ - x nagyon közel esnek egymáshoz, akkor megállunk
        % x_uj = x_regi + d
        % d = x_uj - x_regi
        
        % ne hasznaljunk inv(M)-et
        d = B*x + r - x;
        
        % Ha d kicsi, azaz x_n es x_{n+1} kozel vannak, 
        % azaz x_n es f(x_n) kozel vannak,
        % azaz kozel vagyunk a fixponthoz.
        diff_a = norm(d, "inf");
        diff_r = diff_a / norm(x, "inf");
        
        x = x + d;
        
        % ha az elteres pici, akkor megallunk
        if diff_a < atol || diff_r < rtol
            return
        end
    end
end

% 4. feladat: Implementáljuk a 4 releváns módszert az előző feladat megoldásának felhasználásával.


function [x, diff_a, diff_r, step_num] = jor_iteracio(A, b, x0, omega, atol, rtol, max_steps)
    assert(abs(omega) > 1e-32)
    
    M = diag(diag(A))/omega;
    N = M - A;
    [x, diff_a, diff_r, step_num] = iteracio(M, N, b, x0, atol, rtol, max_steps);
end

function [x, diff_a, diff_r, step_num] = jacobi_iteracio(A, b, x0, atol, rtol, max_steps)
    [x, diff_a, diff_r, step_num] = jor_iteracio(A, b, x0, 1, atol, rtol, max_steps);
end

function [x, diff_a, diff_r, step_num] = sor_iteracio(A, b, x0, omega, atol, rtol, max_steps)
    assert(abs(omega) > 1e-32)
    
    M = diag(diag(A))/omega + tril(A, -1);
    N = M - A;
    [x, diff_a, diff_r, step_num] = iteracio(M, N, b, x0, atol, rtol, max_steps);
end

function [x, diff_a, diff_r, step_num] = gs_iteracio(A, b, x0, atol, rtol, max_steps)
    [x, diff_a, diff_r, step_num] = sor_iteracio(A, b, x0, 1, atol, rtol, max_steps);
end

% Szorgalmi feladat: Csináljuk végig a 2. feladatot JOR helyett SOR iterációval.

%