Numerikus módszerek '12

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, ami sokszor fellép optimalizálási problémák vagy parciális differenciál egyenletek megoldásakor. A mátrix előbbi tulajdonságai 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.