# Gyakorlat 2019. 02. 27.

# Ennek kesobb jon a magyarazata, 
# a megjelenitett tizedesjegyek szamat szabalyozza.
#
# Fontos: ez nem a szamolas pontossaga, az 16 tizedesjegy korul van,
# de altalaban 3-ra kerekitve latjuk az eredmenyeket
output_precision(3)

# Matrixok
A = [2 3; 4 5]

# A transzponaltja
disp("\n\
A transzponaltja: transpose(A)")
transpose(A)

# inverz
disp("\n\
A inverze: A^-1 vagy inv(A)")
A^-1
inv(A)

# determinans
disp("\n\
A determinansa: det(A)")
det(A)

# kondicios szam
disp("\n\
A kondicios szama: cond(A) = ||A||*||A^-1||")
cond(A)

disp("\n\
Ha cond(A) ~ 10^k akkor kb k tizedesjegy pontossagot vesztunk \n\
linearis egyenletrendszer numerikus megoldasanal")

# 1. Pelda
disp("\n\
1. Pelda:")
A = [1 1; 1 1.001]

disp("Ax = [2; 2] megoldasa A\\[2; 2]")
A\[2; 2]

disp("Ax = [2; 2.001] megoldasa A\\[2; 2.001]")
A\[2; 2.001]

disp("\n\
Problema: kicsit valtozott a jobb oldal,\
\n\
mig a megoldas nagyot, nezzuk a kondicios szamot!")

disp("\n\
cond(A): ")
cond(A)

disp("\n\
Ez mar 4 * 10^3 ami nagynak szamit");

# 2. Pelda
disp("\n\
2. Pelda:")

disp("\n\
6.73 x - 8.97 y = 5.61 \n\
4.79 x - 6.39 y = 3.99")
A = [6.73 -8.97; 4.79 -6.39]

b = [5.61; 3.99]

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

disp("\n\
A kondicios szama: cond(A)")
cond(A)

disp("\n\
Valtoztassunk kicsit az egyutthatokon: 6.73 -> 6.72")
A = [6.72 -8.97; 4.79 -6.39]

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

disp("\n\
A kondicios szama: cond(A)")
cond(A)

disp("\n\
Valtoztassunk kicsit az egyutthatokon: -8.97 -> -8.96")
A = [6.72 -8.96; 4.79 -6.39]

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

disp("\n\
A kondicios szama: cond(A)")
cond(A)

disp("\n\
Valtoztassunk kicsit a jobb oldalon: 5.61 -> 5.60")
b = [5.60; 3.99]

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

disp("\n\
Valtoztassunk kicsit a jobb oldalon: 3.99 -> 4.00")
b = [5.60; 4.00]

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

disp("\n\
Latjuk, hogy nagyon kis valtozasok\n\
relative nagy elterest eredmenyeznek a megoldasokban.\n\
Ha ez az egyenletrendszer egy valodi problema, zajos meresek\n\
eredmenyekeppen all elo, akkor latjuk, hogy a modellel baj van.")

# 3. Pelda
disp("\n\
3. Pelda:")
A = [.001 1; 1 1]

b = [1; 2]

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

disp("\n\
A kondicios szama: cond(A)")
cond(A)

disp("\n\
Vegezzuk el a Gauss eliminaciot sorcsere nelkul!")

disp("\n\
Kibovitett matrix:")
Ab = [A(1,:) b(1,:) ; A(2,:) b(2,:)]

disp("\n\
1. Lepes: 2. sorbol az 1. sor 1000-szeresenek kivonasa")
Ab(2,:) = Ab(2,:) - 1000 * Ab(1,:)

disp("\n\
A modosult egyutthatomatrix:")
Ab(1:2,1:2)

disp("\n\
Ennek kondicios szama:")
cond(Ab(1:2,1:2))

disp("\n\
Hoppa! Egy relative jo egyenletrendszerbol az eliminacios lepessel\n\
egy nagyon rosszul viselkedot kaptunk.")

disp("\n\
Ezt demonstralhatjuk:\n\
Visszahelyettesitesnel egy osztassal kezdunk:")

disp("\n\
2. sor * (-1/999)")

Ab(2,:) = Ab(2,:) * (-1/999)

disp("\n\
Nem latunk eleg sok tizedesjegyet! output_precision(10)")
output_precision(10)

Ab

disp("\n\
\"Nem vesztunk\" sokat ha kerekitunk:")

Ab(2,3) = 1

disp("\n\
1. sor - 2. sor")
Ab(1,:) = Ab(1,:) - Ab(2,:)

disp("\n\
Azaz Gauss eliminacioval es egy \"ertelmes\" kerekitessel\n\
a (0, 1) megoldast kapjuk.")

disp("\n\
Idezzuk fel, hogy az eredeti megoldas az ~ (1, 1)")
A\b

disp("\n\
Altalanos cel:\n\
Probaljuk meg elkerulni, hogy az elminacios lepesek soran\n\
a kondicios szam sokat romoljon.")

disp("\n\
Igy a modositott eljaras:\n\
Gauss eliminacio reszleges foelemkivalasztassal (pivoting)")

disp("\n\
A foelem az az elem amivel az azonos oszlop tobbi elemet 0-zzuk.")

disp("\n\
Reszleges foelemkivalasztas: \n\
Az eppen szobajoheto elemek kozul a legnagyobb abszolut ertekut\n\
valasztjuk es sorcserevel a megfelelo helyre hozzuk.")

disp("\n\
Lassuk ujra a feladatot:")
A
b
Ab = [A(1,:) b(1,:) ; A(2,:) b(2,:)]

disp("\n\
Az elso oszlopban az 1 nagyobb abszoluterteku mint a 0.001,\n\
igy a masodik sor elso eleme lesz a foelem -> \n\
1. Lepes: 1. es 2. sor csereje.")

disp("\n\
Ket dolog csereje tipikusan egy 3. ideiglenes (temporary)\n\
valtozo bevezetesevel tortenik.")
temp = Ab(1,:)
Ab(1,:) = Ab(2,:)
Ab(2,:) = temp

disp("\n\
2. Lepes: 2. sorbol az elso 1/1000-szeresenek kivonasa")
Ab(2,:) = Ab(2,:) - Ab(1,:) / 1000

disp("\n\
A modosult egyutthatomatrix:")
Ab(1:2,1:2)

disp("\n\
Ennek kondicios szama:")
cond(Ab(1:2,1:2))

disp("\n\
Ez mar nem romlott!")

disp("\n\
Visszahelyettesites:")
Ab(2,:) = Ab(2,:) / 0.999

disp("\n\
Jo a kondicios szam, itt mar kerekithetunk")
Ab(2,3) = 1

disp("1. egyenletbe visszahelyettesites")
Ab(1,:) = Ab(1,:) - Ab(2,:)

disp("Igy a kerekites utan leolvashato megoldas (1,1), \n\
idezzuk fel, hogy az igazi megoldas hasonlo:")
A\b

#########################################################
disp("\n##############################################################")

# Matrixok felbontasai
disp("Matrixok LR / LU felbontasa")

# LR, LU felbontas: also es felso triangularis szorzata
#
# ha nemszingularis a matrix letezik (esetleges sorcserek utan)

# 1. pelda
disp("\n\
1. Pelda:")
A = [ 2 1; 1 2]

# A ~ LU (also * felso)
disp("\n\
P * A = L * U, \n\
ahol L also, U felso triangularis, P sorcsereket kodol.\n\n\
[L, U, P] = lu(A)")
[L, U, P] = lu(A)

disp("\n\
A = L * U")
L * U

# 2. pelda
disp("\n\
2. Pelda:")
A = [2 3; 3 5]
[L, U, P] = lu(A)

disp("Ekkor A = LU nem all fenn a P permutalo matrix sorcseret jelez")
A

disp("\n\
L * U")
L * U

disp("P * A = L * U teljesul, hiszen P * A = ")
P * A

disp("\n\
3. Pelda:")
B = [2 -1 0; -1 2 -1; 0 -1 2]

disp("P * B = L * U")
[L, U, P] = lu(B)

disp("P egysegmatrix, igy B = L * U")

disp("\n\
Ezert B * U^-1 = L")
B * U^-1

disp("\n\
Triangularis matrixok szorzata triangularis marad.")

disp("\n\
U*U*U*U*U")
U*U*U*U*U

disp("L*L*L*L*L")
L*L*L*L*L

#########################################################
disp("\n##############################################################")

# Matrixok mint diszkret dinamikai rendszerek elemei
disp("Matrixok diszkret dinamikai rendszerekben")

disp("\n\
Elkepzelt szituacio:\n\
- falu lakossaga 100 fo\n\
- kitor a zombi fertozes\n\
- hetente (idoegysegenkent) vizsgaljuk a falut\n\
- egy idoegyseg alatt:\n\
  - 80%-a az embereknek ember marad\n\
  - 20%-a az embereknek zombi lesz\n\
  - 99%-a a zombiknak zombi marad\n\
  -  1%-a a zombiknak meggyogyul (van gyogymod csak nehezkes!)\n\
- egy fertozott egyed inditja a betegseget, azaz\n\
  a legelejen van 99 emberunk es 1 zombink.\n\
  Ezt az adatot reprezentalhatjuk egy tablazatban is:\n\
             letszam  \n\
  emberek |     99   |\n\
  zombik  |      1   |\n\n\
- azaz az adat tulajdonkeppen egy vektor:\n\
  [ e ; z ] ahol e az emberek szama, z a zombik szama.")
[99; 1]

disp("\n\
Egy idoegyseg utan igy alakul a populacio: \n\
e -> 0.8  * e + 0.01 * z\n\
z -> 0.99 * z + 0.2  * e")

disp("\n\
Ez vegyuk eszre, hogy nem mas mint egy baloldali szorzas\n\
az [e; z] vektorra a kovetkezo matrixszal:")
A = [0.8 0.01; 0.2 0.99]

disp("Ez a matrix kodolja a populacio evoluciojat!")

disp("\n\
Nezzuk, hogy alakul a falu sorsa:")
populacioKezdeti = [99; 1]

disp("\n\
1 het utan: A * populacioKezdeti")
A * populacioKezdeti

disp("\n\
2 het utan: A * A * populacioKezdeti")
A * ans

disp("\n\
3 het utan: A * A * A * populacioKezdeti")
A * ans

disp("\n\
Altalanosan n het utan A^n * populacioKezdeti")

disp("\n\
Tehat, 10 het utan:")
A^10 * populacioKezdeti

disp("\n\
Tehat, 20 het utan:")
A^20 * populacioKezdeti

disp("\n\
Tehat, 30 het utan:")
A^30 * populacioKezdeti

disp("\n\
Tehat, 40 het utan:")
A^40 * populacioKezdeti

disp("\n\
Tehat, 50 het utan:")
A^50 * populacioKezdeti

disp("\n\
Tehat, 100 het utan:")
A^100 * populacioKezdeti

disp("\n\
Ugy tunik nem all be lenyegi valtozas, stabilizalodik az atlagosan\n\
-  4.7 ember\n\
- 95.2 zombi")

disp("\n\
Hogy tudnank errol meggyozodni, hogy tenyleg igy van?")

disp("\n\
Matrixok sajatertekfelbontasa:\n\
[V, L] = eig(A)")
[V, L] = eig(A)

disp("\n\
V: az A matrix sajatvektorai\n\
L: az A matrix sajatertekei az atloban")

disp("\n\
1. sajatertek L(1,1)")
L(1,1)

disp("\n\
Ehhez tartozo sajatvektor V elso oszlopa:")
V(:,1)

disp("\n\
L(1,1) * 1. oszlop")
L(1,1) * V(:,1)

disp("\n\
A * 1. oszlop")
A * V(:,1)

disp("\n\
2. sajatertek L(2,2)")
L(2,2)

disp("\n\
Ehhez tartozo sajatvektor V masodik oszlopa:")
V(:,2)

disp("\n\
L(2,2) * 2. oszlop")
L(2,2) * V(:,2)

disp("\n\
A * 2. oszlop")
A * V(:,2)

disp("\n\
az A matrix sajatertekfelbontasa: A = V * L * V^-1")
V * L * V^-1

disp("\n\
Igy A^n = (V * L * V^-1) * (V * L * V^-1) * ... * (V * L * V^-1) = \n\
V * L * (V^-1 * V) * L * (V^-1 * V) * ... * L * V^-1 = \n\
V * L * L * ... * L * V^-1 = \n\
V * L^n * V^-1")

disp("\n\
Azaz eleg lenne tudni, hogy L n-edik hatvanya hogyan alakul.")

disp("\n\
L")
L

disp("\n\
L^2")
L^2

disp("\n\
L^3")
L^3

disp("\n\
diagonalis matrixnal az atloban levo elemek hatvanyozodnak.\n\
0.79^n -> 0, 1^n ->1, igy L^n -> [0 0; 0 1]")

disp("\n\
L^10")
L^10

disp("\n\
L^100")
L^100

disp("Tegyuk vissza a 3 tizedesjegyni megjelenitest")
output_precision(3)
L^100

disp("\n\
Igy L^n -> Ln = [0 0; 0 1]")
Ln = [0 0; 0 1]

disp("\n\
Igy a falu lakossaga \"sok\" idoegyseg utan\n\
V * Ln * V^-1 * populacioKezdeti:")
V * Ln * V^-1 * populacioKezdeti

disp("\n\
Azaz valoban jol lattuk fentebb, hogy ez az allapot stabilizalodik.")