Numerikus módszerek

Differenciálegyenlet-megoldó módszerek

A differenciálegyenlet-megoldó numerikus módszereknek alapvetően két különböző csoportjük van, az explicit és implicit módszerek. Mindkét esetben arra keressük a választ, hogyha ismert a rendszer állapota egy $t_n\,$ időpontban, valamint tetszőleges $t\,$-re elő tudjuk állítani a deriváltakat, akkor mi a lehető legjobb becslés a rendszer jellemzőire (változóira) egy $t_n + \Delta t\,$ időpontban.

Explicit módszerek

Az explicit módszerek egy zártan kiértékelhető (értsd: behelyettesítesz és kész) alakot adnak az $t_n + \Delta t\,$-beli értékekre egy, vagy néhány korábbi pontban ismert állapotokból, valamint az ott ismert deriváltakból.

Euler-módszer

A legegyszerűbb egylépéses módszer. Az y(x=0)=y0 kezdőfeltétellel megadott, y’=f(x,y) diffegyenlet megoldása esetén az Euler lépés alakja (Taylor-sorfejtés első két tagja): $y_(n+1)=y_n+h*f(x_n,y_n)$

Hibája: Taylor-sorfejtést tovább írjuk, a különbség $O(h^2)$ lesz. Használata nem javasolt, mert a hibák hamar felösszegződnek, a megoldás „felrobban”. Ennek elkerülésére érdemes lehet használni az implicit Euler-módszert: $y_(n+1)=y_n+h*f(x_{n+1},y_{n+1})$. Ez nagy h értékekre is stabil marad.

Runge-Kutta módszer

Miért használjuk? Mert sokkal pontosabb, mint az Euler-módszer.

Másodrendű RK (vagy midpoint method - középponti módszer)

$k_1=h\,f(x_n,y_n)\,$

$k_2=h\,f \left ( x_n+\frac{1}{2} h,y_n+\frac{1}{2} k_1 \right ) \,$

$y_{n+1}=y_n+k_2+O(h^3)\,$

Ez a módszer tehát harmadrendig pontos. Általánosan az n-ed rendű RK-nak $O(h^{n+1})$ hibája van.

Negyedrendű RK

$k_1=h\,f(x_n,y_n)\,$

$k_2=h\,f \left ( x_n+\frac{1}{2} h,y_n+\frac{1}{2} k_1 \right ) \,$

$k_3=h\,f \left ( x_n+\frac{1}{2} h,y_n+\frac{1}{2} k_2 \right ) \,$

$k_4=h\,f(x_n+h,y_n+k_3)\,$

$y_{n+1} = y_n + \frac{1}{6} k_1 + \frac{1}{3} k_2 + \frac{1}{3} k_3 + \frac{1}{6} k_4 + O(h^5)$

A negyedrendű módszerben négyszer kell kiértékelni az f függvényt, míg az Euler-módszernél egyszer kellett. Ezért ennek a használata akkor gazdaságos, ha ugyanakkora pontosság mellett legalább négyszer akkora lehet a lépésköz.

Adaptív RK

Egy differenciálegyenlet megoldása során lehetnek gyorsan és lassan változó szakaszok a függvényben. A lassan változó szakaszok integrálása során nagyobb lépéseket is tehetünk a hiba növekedése nélkül. Ennek a megoldására szolgál az adaptív Runge-Kutta módszer. Alapötlete, hogy egy lépést tegyünk meg egyszer teljesen (2h-val, $y_1$), egyszer pedig két fél lépésben (h-val, $y_2$). Mindegyik lépés 4 függvény kiértékelést igényel (3*4), de ebből kettő megegyezik, így 11 kiértékelés szükséges a 2*4 helyett, ami a két fél lépésből jönne össze.

$y(x+2h) = y_1 + (2h)^5 \Phi + O(h^6)\,$

$y(x+2h) = y_1 + 2(h^5) \Phi + O(h^6)\,$

A kettő közti különbség:

$\Delta = y_2 - y_1\,$

A különbség $h^5$-nel skálázik. Ha egy $h_1$ lépés eredménye $\Delta_1$, és mi $\Delta_0$ hibát akarunk elérni, akkor $h_0$ lépést kell tennünk, ami: $h_0 = h_1 \left | \frac{\Delta_0}{\Delta_1} \right | ^{0.2}$

Ha $h_1>h_0$, akkor meg kell ismételni a számolást egy kisebb lépéssel, ha pedig $h_1<h_0$, akkor a következő lépésben használhatjuk $h_0$-t lépésként.

Implicit módszerek

Az implicit módszerek az eddigiekkel szemben más logikán alapszanak: tegyük fel, hogy a $t_n + \Delta t\,$-pontban vagyunk, és ezzel fejezzük ki az itteni deriváltakat, és arra keressük a választ, hogy milyen értékeire az egyenletrendszernek leszünk konzisztensek a korábbi, $t_n\,$-beli értékekkel. Ez tehát egy bonyolultabb feladat, tipikusan egy egyenletrendszert kell megoldanunk minden lépésben, ami jobb esetben lineáris.

Miért érdekesek ezek a módszerek? Azért, mert sokkal stabilabbak (esetleg feltétel nélkül stabilak!), mint az explicit módszerek, ahol nem választhatunk akármekkora lépést, mert akkor a megoldás hasnzálhatatlan lesz (elszáll).

Egyenletrendszerek megoldása

Lineáris egyenletrendszernek nevezzük az olyan egyenletrendszert, melyben minden változó az első hatványon szerepel. Ezeket felírhatjuk mátrixformában is: $A\vec{x} = \vec{b}$. Az A mátrixban tároljuk az egyenletrendszerben szereplő együtthatókat, $\vec{x}$ vektor az ismeretlenek vektora, $\vec{b}$ vektorban pedig az egyenletek jobb oldalainak értékei szerepelnek. Célunk meghatározni az $\vec{x}$-t, ehhez beszorozzuk jobbról az előző egyenletet A^-1-zel: $\vec{x} = A^{-1}\vec{b}$. A feladatunk tehát az A mátrix inverzét meghatározni.

Gauss-Jordan elimináció

A Gauss-elimináció során az A mátrix és a $\vec{b}$ soraival egyszerre végzünk transzformációkat (másképp fogalmazva: az egyenleteket alakítjuk át). A megengedett transzformációk:

• mátrix két sorának felcserélése (két egyenlet cseréje)
• mátrix sorának számmal való szorzása
• a mátrix egy sorához egy másik skalárszorosának hozzáadása

A cél, hogy ilyen transzformációk sorozatának alkalmazásával A mátrixból egy felső háromszög mátrixok csináljunk. Ezután alulról fölfelé haladva behelyettesítéssel megkapjuk sorban $\vec{x}$ elemeit. Példa

LU dekompozíció

Legyen A egy kvadratikus mátrix, erre az LU felbontás: A = LU. Az L és U mátrixok alsó illetve felső háromszög mátrixok (főátlóban is vannak elemek). Példa

Alkalmazása

1. Determináns számolás: det(A) = det(L)*det(U)

A háromszög mátrixok determinánsa a főátlóban lévő elemek szorzata, ezért ez nagyon egyszerűen számolható

1. Mátrix inverz: $A^{-1} = U^{-1}L^{-1}$

Háromszög mátrix inverze könnyen számolható, ezért jó ez a felbontás.

Singular Value Decomposition - Szinguláris érték dekompozíció

Sokszor előfordul, hogy a mátrix, amit invertálni akarunk, közel van a szingulárishoz (det A ≈ 0), és sem a Gauss-elimináció, se az LU felbontás nem jár sikerrel.

Adott egy A mátrixunk, ami MxN méretű. Ezt a mátrixot felírhatjuk UWV^T formában, ahol:

• U: MxN, oszlop-ortogonális mátrix
• W: NxN, diagonális, a szinguláris értékeket tartalmazza
• V^T: NxN, ortogonális mátrix

A módszer segítségével nem-négyzetes mátrixokat is tudunk invertálni (pszeudoinverz) a következő módon:

$A^{-1} = V\,[diag(1/w_j)]\,U^T$

Ami problémát okozhat, az a szinguláris értékek reciproka, ha az nulla, vagy nagyon közel van nullához. A megoldás, hogy ebben az esetben az 1/w_j-t nullának vesszük.

Optimalizációs módszerek - szélsőérték keresés

Newton-iteráció

Bár eredetileg a Newton-iterációt gyökkeresésre szokás alkalmazni, azonban ha ismert, vagy előállítható a függvény deriváltja, aminek a szélsőértékét keressük, akkor a derivált függvény zérushelye éppen a függvény szélsőértékét adja. A Newton-iteráció alakja szélsőérték keresésre:

$x_{n+1} = x_n - \frac{f'(x_n)}{f''(x_n)}$

ahol ' a deriválást jelöli. Látszik, hogy a második deriváltra is szükségünk van. A Newton-módszer kiterjeszthető többváltozó esetére is, ekkor az első deriváltak vektora és a második deriváltakat tartalmazó mátrixra van szükségünk az analóg formula kiértékeléséhez.

Konjugált gradiens módszer (iteratív változat)

A konjugált gradiens módszer lineáris egyenletrendszert tud megoldani, abban az esetben, ha az együtthatómátrix szimmetrikus és pozitiv definit. Ez szemléletesen annak felel meg, hogy egy olyan sokdimenziós potenciálban keresünk minimumot, amelynek nincs nyeregfelülete, hanem egyértelmű minimuma van. Legyen az egyenletrendszer: $Ax = b\,$. Belátható, hogy ennek a megoldása a fenti A-ra vonatkozó feltételek mellett minimalizálja a következő függvényt:

$f(x) = \frac{1}{2}x^T A x - x^T b\,$.

Vezessük be a konjugáltság fogalmát két vektorra:

$v_1^T A v_2 = 0\,$.

Ha $A\,$ $N\,$ dimenziós, akkor létezik ennyi darab konjugált vektor (jelöljük $d_i$-vel), amelyek lineárisan függetlenek, tehát bázist alkotnak, ezért minden vektor kifejthető rajtuk, speciálisan a megoldás x is: $\alpha_i$ koefficiensekkel. Ezt behelyettesítve az egyenletrendszerbe kapjuk, hogy:

$\alpha_i = -\frac{d_i^T b}{d_i^T A d_i} \,$.

Tehát a cél ezek után $d_i$-k előállítása. Az ötlet az, hogy állítsuk elő sorban ezeket az irányokat, egyenként minimalizálva, egyre nagyobb altereiben A-nak. Konkrétan az állítás: a következő $x_i\,$ sorozat a megoldáshoz konvergál:

$x_{i+1} = x_i + \alpha_i d_i\,$

ahol a kifejtési együttható:

$\alpha_i = -\frac{g_i^T d_i}{d_i^T A d_i}\,$

és:

$d_{i} = -g_{i} + \frac{g_{i}^T A d_{i-1}}{d_{i-1}^T A d_{i-1}} d_{i-1}\,$, de az első lépésben egyszerűen a legmeredekebb irányba lépünk: $d_{0} = -g_{0}\,$

végül a gradiens irány egyszerűen a derivált kiértékelése:

$g_i = Ax_i + b\,$,

Szemléletesen az történik, hogy ahányszor előállítunk egy a minimum felé mutató gradiens irányt, abból levonjuk a korábbi irányok konjugáció által definiált vetületét. Tehát minden lépésben egy irányban megkeressük az abban az irányban elérhető legjobb minimumot és a későbbi lépéseket úgy választjuk meg, hogy erre az irányra merőlegesek legyenek, így nem rontják el a korábban elért rész-minimumokat. Ebből az is látszik, hogyha mind az N irányban elvégeztük ezeket a lépéseket, akkor elvileg a minimumban kell lennünk (nincs több független irány, amiben lehetne minimumot keresni). A gyakorlatban a kerekítési hibák ebben megakadályoznak.

Szimulált hőkezelés

A módszer elnevezése az anyagtudomány területéről ered, ahol hőkezelési eljárásokkal lehet változtatni egy anyag kristályosodási méretét. Számítógépes módszerként egy sokdimenziós függvény energiaminimumának megkeresésére lehet használni.

Működési struktúra:

1. Válasszunk egy tetszőleges kezdőállapotot (i)
2. Válasszuk ki az egyik szomszédot (j)
3. Ha a szomszéd energiája kisebb, átlépünk rá
4. Ha nagyobb a szomszéd energiája, az energiakülönbség, és egy globális T paraméter által meghatározott valószínűség szerint elfogadjuk a nagyobb energiájú helyet

Az elfogadási valószínűségek tehát:

• ha f(j) ≤ f(i), akkor 1
• ha f(j) > f(i), akkor $e^{\frac{f(i)-f(j)}{T}}$

A T hőmérsékletet kezdetben nagynak választva, majd folyamatosan csökkentve ezzel a módszerrel megtalálhatjuk a függvény globális minimumát.