# Gyakorlat 2019. 03. 13.
output_precision(5)

disp("Gyakorlat 2019. 03. 13.");

disp("\n\
Ax = b megoldasa iterativ modszerekkel")

disp("\n\
Alap otlet F(x*) = x* egyenlet megoldasa,\n\
azaz F fixpontjanak keresese")

disp("\n\
Ilyen egyenletek eseten gyakran hasznalnak iterativ modszert,\n\
pl F iteralasat:\n\
y0 - kezdeti becsles\n\
y1 = F(y0)\n\
y2 = F(y1)\n\
...\n\
y(n+1) = F(y(n))")

disp("\n\
Ez az ugynevezett fixpont iteracio.\n\
Sok esetben konvergens,\n\
azaz a sorozat tart a fixponthoz: y(n)->x*.\n\
Megjegyzes: tipikusan akkor, ha a fixpont stabil es\n\
a kezdeti becsles a fixpont stabil halmazaba esik.\n\
Ha a fixpont globalisan aszimpotikusan stabil, akkor \n\
barmilyen kezdeti becsles megfelel.")

disp("\n\
Otlet: Irjuk at az Ax = b egyenletet ugy, hogy\n\
annak megoldasa fixpontja legyen egy F(x) = x egyenletnek")

disp("\n\
Motivacio: Ax = b numerikus megoldasa varhatoan tartalmaz\n\
szamitasi hibat hiszen geppel dolgozunk.\n\
Tehat alapbol el kell fogadni valamilyen merteku\n\
pontatlansagot az eredmenyben.\n\n\
A klasszikus eliminacio felbehagyasakor nem keletkezik\n\
hasznalhato reszeredmeny, azaz\n\
az eliminacio muveletkoltsege nem sporolhato meg.\n\n\
Ha ugyis elfogadunk szamitasi hibakat, akkor\n\
olyan modszert is elfogadnank ami kevesebb muvelettel jar\n\
mint az eliminacio, de mi donthetjuk el, hogy\n\
mennyi energiat fektetunk be.\n\
(tipikusan tobb energia ~ pontosabb eredmeny)\n\n\
Erre nyujt lehetoseget egy iterativ modszer, ahol\n\
remenyeink szerint - azaz ha konvergal az iteracio - \n\
tobb iteracio pontosabb eredmenyt nyujt, de\n\
ha mar ugy latjuk, hogy a sorozat \"konvergal\"\n\
akkor abba is lehet hagyni es a reszeredmenyt\n\
mitn kozelito megoldast lehet kezelni.\n\n\
Ezeket altalaban nagymeretu Ax = b egyenletekre\n\
szoktak alkalmazni, ahol az eliminacio \"draga\"\n\
Kulonoskeppen akkor ha amugy A-ban sok a 0\n\
azaz A ugynevezett ritka (sparse) matrix.\n\n\
Ket tipikus veszelyt rejt magaban ez a modszer:\n\
1) az iteracio lehet, hogy sosem konvergal\n\
2) latszolag tevesen konvergalhat az iteracio\n\
   egy rossz eredmenyhez (s akar kesobb a johoz is)")

disp("\n\
Pelda:\n\
 4 x1 + 2 x2 + 3 x3 =   8\n\
 3 x1 - 5 x2 + 2 x3 = -14\n\
-2 x1 + 3 x2 + 8 x3 =  27")
A = [4 2 3; 3 -5 2; -2 3 8]
b = [8 ; -14; 27]

disp("\n\
Megoldasa: ")
A \ b

disp("\n\
Ket iteracios modszert tanulunk: Jacobi es Gauss-Seidel iteraciok")

disp("\n\
Mindketto hasonlokeppen alkotja meg a fixpont egyenletet:")

disp("\n\
Az 1. egyenletbol kifejezzuk x1-et:\n\
4 x1 + 2 x2 + 3 x3 =   8\n\
=>\n\
x1 = ( 8 - (2 x2 + 3 x3) ) / 4")
disp("\n\
A 2. egyenletbol kifejezzuk x2-ot:\n\
3 x1 - 5 x2 + 2 x3 = -14\n\
=>\n\
x2 = ( -14 - ( 3 x1 + 2 x3 ) ) / -5")
disp("\n\
A 3. egyenletbol kifejezzuk x3-at:\n\
-2 x1 + 3 x2 + 8 x3 =  27\n\
=>\n\
x3 = ( 27 - ( -2 x1 + 3 x2 ) ) / 8")
disp("\n\
Igy: F([x1; x2; x3] = ")
[
"( 8 - (2 x2 + 3 x3) ) / 4";
"( -14 - ( 3 x1 + 2 x3 ) ) / -5";
"( 27 - ( -2 x1 + 3 x2 ) ) / 8"
]
disp("\n\
es nyilvan az F(x*) = x* megoldasa teljesiti, hogy Ax* = b")

disp("\n\
Az F fuggveny felepitese Octaveban a komponensek fuggvenyeinek\n\
megadasaval:")

# helyi fuggvenyek definialasa
disp("\n\
x1-re, ha x2 = 0, x3 = 0")
function eredmeny = x1uj(x2, x3)
  eredmeny = ( 8 - (2 * x2 + 3 * x3) ) / 4;
endfunction
x1uj(0, 0)

disp("\n\
x2-re, ha x1 = 0, x3 = 0")
function eredmeny = x2uj(x1, x3)
  eredmeny = ( -14 - ( 3 * x1 + 2 * x3 ) ) / -5;
endfunction
x2uj(0, 0)

disp("\n\
x3-ra, ha x1 = 0, x2 = 0")
function eredmeny = x3uj(x1, x2)
  eredmeny = ( 27 - ( -2 * x1 + 3 * x2 ) ) / 8;
endfunction
x3uj(0, 0)

disp("\n\
====================\n\
Jacobi iteracio:\n\
- egy kezdeti becsles\n\
- majd F kiertekelesenek ismetlese (iteralasa) a becslesen\n\n\
[x1uj; x2uj; x3uj] = F([x1; x2; x3] = \n\
[\n\
( 8 - (2 x2 + 3 x3) ) / 4;\n\
( -14 - ( 3 x1 + 2 x3 ) ) / -5;\n\
( 27 - ( -2 x1 + 3 x2 ) ) / 8\n\
]\n\n\
Elonyok:\n\
- parhuzamosithato, minden egyenletet lehet kulon szamolni\n\
- tobbszor konvergal mint a kesobb bemutatatando Gauss-Seidel\n\n\
Hatranyok:\n\
- lassabban konvergal mint a Gauss-Seidel\n\
- tobb tarhelyet igenyel hisz egy iteracios lepesben minden\n\
  egyenlet megoldasahoz ismerni kell az elozo lepes eredmenyet.")

disp("\n\
Jacobi iteracio bemutatasa: Ax = b")
A
b
disp("\n\
Kezdeti becsles:")
x = [0; 0; 0]
xuj = [
x1uj(x(2), x(3));
x2uj(x(1), x(3));
x3uj(x(1), x(2))
]

disp("\n\
Figyeljuk meg, hogy egyszerre kell x-et es xuj-at tarolni,\n\
azonban xuj komponenseit egymastol fuggetlenul is tudjuk szamolni")

disp("\n\
Iteracio: x = xuj es ujra:")
x = xuj
xuj = [
x1uj(x(2), x(3));
x2uj(x(1), x(3));
x3uj(x(1), x(2))
]
x = xuj;

disp("\n\
Tovabbi 10 lepes utan: (koztes eredmenyeket nem irjuk ki")
# az eredmeny nem irodik ki ha ; zarja a sort -> probald meg levenni a ;-t
# for ciklus: az i valtozo ertekei 1-tol 10-ig valtoznak es minden egyes ertekre
# vegrehajtodik a ciklus teste
for i = [1:10]
  xuj = [
    x1uj(x(2), x(3));
    x2uj(x(1), x(3));
    x3uj(x(1), x(2))
  ];
  x = xuj;
endfor
x

disp("\n\
Tovabbi 10 lepes utan: (koztes eredmenyeket nem irjuk ki")
# az eredmeny nem irodik ki ha ; zarja a sort -> probald meg levenni a ;-t
for i = [1:10]
  xuj = [
    x1uj(x(2), x(3));
    x2uj(x(1), x(3));
    x3uj(x(1), x(2))
  ];
  x = xuj;
endfor
x

disp("\n\
Tovabbi 10 lepes:")
for i = [1:10]
  xuj = [
    x1uj(x(2), x(3));
    x2uj(x(1), x(3));
    x3uj(x(1), x(2))
  ];
  x = xuj
endfor

disp("\n\
\"Latszik\", hogy -1, 3, 2 stabilizalodik\n\n\
Ez is a megoldas:")
A\b


disp("\n\
====================\n\
Gauss-Seidel iteracio:\n\
- egy kezdeti becsles\n\
- majd F kiertekelesenek ismetlese (iteralasa) a becslesen\n\
  a koztes eredmenyek felhasznalasaval.\n\n\
x1uj = ( 8 - (2 x2 + 3 x3) ) / 4\n\
x2uj = ( -14 - ( 3 x1uj + 2 x3 ) ) / -5\n\
x3uj = ( 27 - ( -2 x1uj + 3 x2uj ) ) / 8\n\n\
Elonyok:\n\
- gyorsabban konvergal mint a Jacobi\n\
- kevesebb tarhelyet igenyel hisz egy iteracios lepesben\n\
  eleg az azon lepes reszeredmenyeit hasznalni\n\n\
Hatranyok:\n\
- nem parhuzamosithato, minden egyenletnek\n\
  szuksege van az elozoek reszeredmenyeire\n\
- elofordul, hogy nem konvergal mikor a Jacobi igen.")

disp("\n\
Gauss-Seidel iteracio bemutatasa: Ax = b")
A
b
disp("\n\
Kezdeti becsles:")
x = [0; 0; 0]

disp("\n\
1. iteracio, 1. egyenlet")
x(1) = x1uj(x(2), x(3))

disp("\n\
1. iteracio, 2. egyenlet")
x(2) = x2uj(x(1), x(3))

disp("\n\
1. iteracio, 3. egyenlet")
x(3) = x3uj(x(1), x(2))

disp("\n\
Figyeljuk meg, hogy egyszerre elegendo x-et tarolni es valtoztatni,\n\
azonban a komponenseket nem tudjuk egymastol fuggetlenul szamolni")

disp("\n\
Iteracio: ujra:")
disp("\n\
2. iteracio, 1. egyenlet")
x(1) = x1uj(x(2), x(3))

disp("\n\
2. iteracio, 2. egyenlet")
x(2) = x2uj(x(1), x(3))

disp("\n\
2. iteracio, 3. egyenlet")
x(3) = x3uj(x(1), x(2))

disp("\n\
Tovabbi 5 iteracio:")
for i = [1:5]
  x(1) = x1uj(x(2), x(3));
  x(2) = x2uj(x(1), x(3));
  x(3) = x3uj(x(1), x(2));
  # i + 2 mert fentebb megcsinaltunk mar ket iteraciot
  disp("iteracio sorszama:"), disp(i + 2)
  disp("a megfelelo becsles a megoldasra:")
  x
endfor

disp("\n\
Figyeljuk meg, hogy sokkal gyorsabban kerultek az iteraltak az\n\
1, -3, 2 megoldas kozelebe mint a Jacobi iteracio soran.")
