Matematika | Analízis » Koltay Anita - Mátrix szeletelés elmélet

http://www.doksihu Eötvös Loránd Tudományegyetem Természettudományi Kar Mátrix szeletelés elmélet Szakdolgozat Koltay Anita Matematika B.Sc, elemz® szakirány Témavezet®: Mincsovics Miklós, tanársegéd Alkalmazott Analízis és Számításmatematikai Tanszék Budapest 2010 http://www.doksihu Tartalomjegyzék 1. Bevezetés 2 2. Példák 4 2.1 Kétpontos peremérték feladat 2.2 A Dirichlet-probléma 3. Klasszikus iterációs módszerek 4 6 9 3.1 Iterációs módszerek konvergenciája 10 3.2 A Richardson-, a Jacobi-, a Gauss-Seidel és a SOR-módszer 15 4. Nemnegatív mátrix szeletelés 21 4.1 Mátrix szeletelés 21 4.2 Prefaktorizációs iterációs módszerek 33 5. Összefoglalás 37 6. Függelék 38 7. Köszönetnyilvánítás 42 1 http://www.doksihu 1. fejezet Bevezetés Számos matematikai

probléma modellezhet® lineáris algebrai egyenletrendszerekkel, például nemlineáris egyenletrendszerek, dierenciál- és integrálegyenletek, interpolációs és optimalizációs feladatok. Egyenletrendszerekre vezethet® le továbbá sok zikai (pl mechanikai, villamosságtani) feladat, ökológiai, gazdasági elemzés. A lineáris algebrai egyenletrendszerek megoldhatók direkt és iterációs módszerekkel is. A direkt módszereknél olyan algoritmust adunk meg, mely véges sok lépésben kiszámolja a pontos megoldást, pontos kiindulási adatokat feltételezve és a kerekítési hibáktól eltekintve (Gauss-elimináció, LU- és Cholesky-felbontás). Gyakorlati feladatoknál a kapott eredmény persze hibával terhelt lesz, mivel már a kezdeti adatok sem pontosak (mérési, modellezési hibák). A direkt módszereket f®leg kis és közepes méret¶, illetve sávmátrixok és teli mátrixok esetén használják. A gyakorlatban sok nagyméret¶, ritka mátrixú

egyenletrendszer keletkezik, ezeknél a direkt módszerek hátránya a nagy tárigény, valamint nagy együtthatómátrixoknál az inverz kiszámítása is nehéz és költséges, ezért ilyen típusú mátrixokra iterációs módszert alkalmaznak. Az iterációs módszerek egy (általunk megadott) kiindulási vektorból egy vektorsorozatot állítanak el®, mely az egyenletrendszer megoldásához konvergál. A gyakorlatban az iterációs lépések számára meg kell adni egy leállási kritériumot, melynek meghatározása újabb problémát vet fel. (Az egyik lehetséges módszer, hogy akkor állítjuk meg az iterációt, amikor az egymást követ® vektorok különbsége már kisebb egy általunk megadott ε-nál.) Ezeknél a módszereknél két fontos szempont van: - a vektorsorozat elemeit könnyen, "olcsón" tudjuk kiszámolni (hiszen pont ez az el®nye a direkt módszerekkel szemben), - az iteráció gyorsan konvergáljon a megoldáshoz (minél gyorsabb a konvergencia,

annál 2 http://www.doksihu kevesebb lépésre lesz szükség). Az iterációs módszereknél az egyenletrendszer együtthatómátrixát olyan mátrixok összegére (ill. különbségére) bontjuk, melyekkel könnyebben és gyorsabban tudunk számolni, mint az eredeti mátrixszal. A klasszikus iterációs módszereknél olyan egyszer¶ felbontásokat vizsgálunk, melyeknél az iteráció gyorsan és olcsón számolható, vagyis kicsi a m¶veletigény, és megnézzük, hogy mely speciális mátrixokra lesznek ezek a felbontások konvergensek. A nemnegatív mátrix szeletelés elmélet (nonnegative splitting theory ) a monoton mátrixú lineáris egyenletrendszerek iterációs megoldási módszereit és azok konvergenciájának összehasonlítását mutatja be. 3 http://www.doksihu 2. fejezet Példák Ebben a fejezetben két olyan matematikai problémát mutatok be, melynek megoldását lineáris algebrai egyenletrendszerrel modelleztem. Az egyenletrendszereket direkt

módszerrel meg is oldottam (Matlab program segítségével), és ábrákkal szemléltetem a pontos és a numerikus megoldást, valamint a hibafüggvényt. 2.1 Kétpontos peremérték feladat Tekintsük a következ® másodrend¶, önadjungált általános dierenciálegyenletet − d2 y(x) + σ(x)y(x) = f (x), dx2 a < x < b, (2.1) a kétpontos peremfeltétellel y(a) = α, y(b) = β, ahol α és β valós konstansok adottak, továbbá f (x) és σ(x) adott valós, folytonos függvények x ∈ [a,b] intervallumon és σ(x) ≥ 0. A feladat folytonos megoldása helyett keressük az alábbi diszkretizációs feladat megoldását. Deniálunk egy h = (b−a) lépésközt és az xj = a + jh, 0 ≤ j ≤ N + 1 rácspontokat. N +1 Bevezetve az yi ≡ y(xi ) jelölést, yi körül negyedrendben Taylor-sorba fejtjük annak szomszédos pontjaiban a függvényt (persze ehhez fel kell tennünk, hogy y(x) legalább négyszer folytonosan dierenciálható): yi±1 = yi ± hyi0 + h2 00

h3 000 h4 IV y ± yi + yi (xi + θi± h), 2! i 3! 4! 4 0 < |θi± | < 1. http://www.doksihu −yi00 = 2yi − yi−1 − yi+1 h2 IV + y (xi + θi h), h2 12 0 < |θi | < 1. Visszahelyettesítve ezt a 2.1 egyenletbe kapjuk, hogy 2yi − yi−1 − yi+1 h2 IV y (xi + θi h), + σ y = f − i i i h2 12 (2.2) 1 ≤ i ≤ N. Deniáljuk az A mátrixot és az y, k, τ vektorokat a következ®képpen:    1   A= 2 h      y  1   y2    y =  . , .  .    yN 2 + σ1 h2 −1 0 2 + σ2 h2 −1 −1 . . −1 0     k=   f1 + α h2 f2 . . fN + β h2     ,        ,    −1 2 + σN h2  y IV (x1 + θ1 h)  IV −h2   y (x2 + θ2 h) τ=  . 12   y IV (xN + θN h)     .   Ekkor a 2.2 egyenletünk felírható az alábbi mátrixos alakban Ay = k + τ (y), ahol τ (y) az O(h2 ) rend¶

hibatag. Tekintsünk egy konkrét példát és annak megoldását: Legyen f (x) = (π 2 + 1)sin(xπ), σ(x) = 1, és a [0, 1] intervallumon h = 111 lépésközzel keressük a diszkretizációs feladat megoldását. Matlab programmal az Ay = k egyenletrendszert megoldva, másodrendben közelítjük y(x)-et Az így kapott numerikus eredmény a pontos megoldással (y(x) = sin(xπ)), valamint a hibával együtt a 2.1 ábrán látható (a 3. kép a tridiagonális A mátrix nemnulla elemeinek elhelyezkedését szemlélteti) A hiba 7 · 10−3 ≈ 10−2 . 5 http://www.doksihu 2.1 ábra 2.2 A Dirichlet-probléma Keressük azt a zárt egységnégyzeten értelmezett u(x, y) függvényt (ill. annak közelítését), mely kielégíti az alábbi egyenletet az egységnégyzet belsejében −  ∂ 2 u(x, y) ∂x2 + ∂ 2 u(x, y)  = g(x, y), ∂y 2 0 < x, y < 1, illetve a peremfeltételt u(x, y) = 0, ha (x, y) ∈ Γ, ahol Γ a négyzet határa, g(x, y) pedig (0,1)×(0,1)-en

értelmezett, adott függvény. 1 Deniálunk egy rácshálót az egységnégyzeten h = n+1 lépésközzel. 6 http://www.doksihu A függvény értékét egy bels® (x0 , y0 ) pontjában (másodrendben) közelítjük a szomszédos pontokban ((x0 ± h, y0 ), (x0 , y0 ± h)) felvett értékekkel, u(x, y) Taylor-sorba fejtésének segítségével, feltéve, hogy a függvény megfelel®en sima, azaz kell®en sokszor folytonosan dierenciálható (ld. Függelék): 1 h2 u(x0 , y0 ) ≈ [u(x0 + h, y0 ) + u(x0 − h, y0 ) + u(x0 , y0 + h) + u(x0 , y0 − h)] + g(x0 , y0 ). 4 4 Most egy olyan Au = k egyenletrendszert kapunk, melyben A 1 tridiag[−I, tridiag[−1, 4, −1], −I] h2 alakú blokk-tridiagonális mátrix (n méret¶ blokkokkal), u vektorban u(x, y) bels® ponA= tokban felvett ismeretlen értékei szerepelnek, k-ban pedig a g(x0 , y0 ) értékek, szintén a bels® (x0 , y0 ) rácspontokban. 1. példa Legyen f (x, y) = 2π 2 sin(xπ)sin(yπ), n = 20 esetén a hiba 10−3

nagyságú (ekkor a pontos megoldás: u(x, y) = sin(xπ) · sin(yπ)): 2.2 ábra 7 http://www.doksihu Ha kisebb lépésközzel (n = 500) közelítjük a függvényt, pontosabb eredményt kapunk (10−6 nagyságrend¶ lesz csak a hiba): 2.3 ábra A diszkretizációs feladatnál minél kisebb lépésközt választunk, annál nagyobb lesz a lineáris algebrai egyenletrendszer együtthatómátrixa és annál nehezebben számolható ki a megoldás direkt módszerekkel. A 23 ábrán látható, hogy már n = 500 esetén is 2.5·105 ×25·105 nagyságú a mátrix Ha háromdimenziós tartományban kell megoldanunk egy, a Dirichlet-problémához hasonló egyenletet, akkor a mátrix mérete n2 -szeresére n®. A gyakorlatban még pontosabb közelítés esetén ennél jóval nagyobb mátrixokkal is kell számolnunk, ekkor a tárhely és a m¶veletigény is jelent®sen megn®, és nem biztos, hogy célszer¶ direkt módszereket alkalmaznunk ilyenkor. 8 http://www.doksihu 3. fejezet

Klasszikus iterációs módszerek Tekintsük az Ax = k lineáris algebrai egyenletrendszert, ahol A egy adott n × n-es valós, reguláris mátrix, x és k pedig n komponens¶ oszlopvektorok (k adott). Keressük az egyenletrendszer x = A−1 k megoldását. Az A mátrixnak egy A = M − N felbontását, ahol M reguláris, egyszer¶ felbontásnak nevezzük, és a fenti egyenletbe behelyettesítve egy általános iterációs módszert kapunk: Mx(m+1) = Nx(m) + k, x(m+1) = M−1 Nx(m) + M−1 k, m = 0, 1, 2., m = 0, 1, 2. (3.1) (3.2) (A gyakorlatban az iterációs lépéseknél 3.1 egyenletrendszert oldják meg direkt módszerekkel, azonban ennek költsége arányos M−1 kiszámításával, ezért mi a továbbiakban csak a 3.2 egyenletrendszert vizsgáljuk) Minél közelebb van A-hoz az M mátrix, a 3.2 iteráció annál gyorsabban konvergál xhez, ugyanakkor annál költségesebb is kiszámolni a vektorsorozat elemeit Ha például M = A, akkor csupán egyetlen lépésre van

szükség a megoldáshoz, viszont M−1 = A−1 -t kell "kiszámolnunk", miközben pont ennek a nehézsége indolkolja az iterációs módszer alkalmazását. A cél tehát az, hogy egy adott mátrixnál olyan felbontásokat találjunk, melyekben M−1 lényegesen könnyebben számolható mint A−1 , és sok iterációs lépés esetén is kisebb legyen a m¶veletigény, mint ha az eredeti egyenletrendszert oldanánk meg. Ugyanakkor teljesülnie kell a konvergenciának is, és persze annak gyorsasága sem elhanyagolható 9 http://www.doksihu 3.1 Iterációs módszerek konvergenciá ja Az iterációs módszerek egy kezdeti vektorból kiindulva lépésr®l-lépésre közelítik a lineáris algebrai egyenletrendszerek megoldását. A gyakorlatban legtöbbször a 0 vektorról indítják az iterációt. A kezdeti vektor ugyan nem befolyásolja egy módszer konvergenciáját, viszont, f®leg rosszul kondícionált egyenletrendszereknél, lelassíthatja azt Térjünk vissza a

3.2 iterációs módszerhez: x(m+1) = Tx(m) + g, m = 0, 1, 2., (3.3) ahol T = M−1 N az iterációs mátrix és g = M−1 k. A módszer konvergens, ha lim x(m) = x m∞ ∀ x(0) kezdeti vektor esetén, azaz lim kεε(m) k = 0, m∞ ahol ε (m) = x(m) − x az m. hibavektor Az f (x) = Tx + g függvény kontrakció (ld. Függelék), ha kTk < 1, mivel kTx1 + g − Tx2 − gk = kT(x1 − x2 )k ≤ kTkkx1 − x2 k. A Banach-féle xponttételb®l (ld. Függelék) kapjuk, hogy f (x)-nek ekkor egyértelm¶en létezik xpontja és az iteráció tetsz®leges kezdeti vektorból indulva ehhez a xponthoz tart. Tehát a kTk < 1 feltétel biztosítja a 3.3 módszer konvergenciáját Másrészt a 3.3 egyenletb®l kapjuk, hogy ε (m) = Tεε(m−1) = . = Tmε (0) , m = 0, 1, 2., vagyis ε (m) 0 ∀εε(0) -ra ⇐⇒ lim Tm = 0. m∞ Egy mátrix konvergenciájában alapvet® szerepe van a spektrálsugárnak, ezért ezt kicsit részletesebben tárgyaljuk. 10 http://www.doksihu A

spektrálsugár 3.1 Deníció Egy A = (ai,j ) n × n-es komplex mátrix spektrálsugarának nevezzük a sajátértékei közül abszolút értékben maximálisat: ρ(A) ≡ max |λi |. 1≤i≤n 3.1 Tétel 1 Bármely vektornorma által indukált mátrixnormára (ld Függelék): ρ(A) ≤ kAk. 2. Ha A szimmetrikus mátrix, akkor kAk2 = ρ(A) 3. Ha A diagonizálható mátrix, akkor létezik olyan normája, melyre kAk = ρ(A) 4. Minden A mátrixhoz tetsz®leges  > 0 esetén létezik olyan norma, melyre ρ(A) ≤ kAk ≤ ρ(A) + . 3.1 Bizonyítás 1Az A mátrix |λ| = ρ(A) sajátértékéhez tartozó v 6= 0 sajátvektorra Av = λv, |λ|kvk = kAvk ≤ kAkkvk, p p 2. kAk2 = ρ(AT A) = ρ(A2 ) = ρ(A) vagyis ρ(A) = |λ| ≤ kAk. 3. Ha S reguláris mátrix és k · k0 tetsz®leges vektornorma, akkor kxk := kSxk0 egy másik vektornorma, és a megfelel® mátrixnormákra teljesül a következ® egyenl®ség: kAk = kSAS−1 k0 . Mivel A diagonizálható, létezik S reguláris

mátrix, melyre SAS−1 diagonális és f®átlójában éppen A sajátértékei vannak. Ha a ∞-normát választjuk k · k0 -nak, akkor az el®z® egyenl®ségb®l kapjuk, hogy kAk = ρ(A) 4. Bizonyítása hasonló 3-hoz, de itt J = SAS−1 az A-nak az a Jordan-féle kanonikus alakja, melynek f®átlójában A sajátértékei állnak, a diagonális fölött pedig -ok (a többi elem 0). Az kxk-normát továbbra is kxk := kSxk∞ -ként deniálva kapjuk, hogy kAk = kJk∞ ≤ ρ(A) + . Most visszatérünk arra, hogy mikor konvergens egy mátrix. 3.2 Tétel Legyen A egy n × n-es valós mátrix lim Am = 0 ⇐⇒ ρ(A) < 1. m∞ 3.2 Bizonyítás kAm k ≥ ρ(Am ) = ρ(A)m , ezért ha ρ(A) ≥ 1, akkor Am nem konvergál 0-hoz. Ugyanakkor, ha ρ(A) < 1, akkor bármely ρ̄ ∈ (ρ(A), 1)-re A-nak létezik olyan normája, melyre kAk ≤ ρ̄, ezért kAm k ≤ kAkm ≤ ρ̄m 0. 11 http://www.doksihu 3.1 Következmény A 33 iterációs módszer akkor és csak akkor

konvergens, ha az iterációs mátrix spektrálsugara, ρ(T) < 1. A ρ(T) < 1 feltétel azonban nem biztosítja a numerikus konvergenciát, az iterácó a számítógépen túlcsordulás miatt leállhat, mivel a Tm 0 konvergencia általában nem monoton. Ezért numerikus szempontból, a számítógépen a kTk < 1 a biztos feltétel A továbbiakban a különböz® iterációs módszerek konvergenciájának sebességét szeretnénk összehasonlítani. 3.2 Deníció Legyen A egy konvergens n × n-es komplex mátrix, R∞ (A) ≡ −lnρ(A)-t a konvergencia aszimptotikus sebességének nevezzük. 3.3 Tétel (Gelfand-formula) Bármely mátrixnormára 1 ρ(A) = lim kAm k m . m∞ Egy mátrix spektrálsugara a gyakorlatban sokszor nehezen számolható ki, ezért bevezetünk egy könnyebben kezelhet® fogalmat is az iterációk összehasonlításához. 3.3 Deníció Legyen A és B két n × n-es komplex mátrix Ha van olyan pozitív valós m, melyre kAm k < 1, akkor R(Am )

≡ −ln[(kAm k)1/m ] = −lnkAm k m az A mátrix konvergenciájának átlagos sebessége m iteráció esetén. Ha R(Am ) < R(Bm ), akkor azt mondjuk, hogy B gyorsabban iterál m lépés esetén mint A. R(Am ) jelentése a következ®: Iterációnként a kezdeti hiba átlagos csökkenése (normában): σ≡  kε(m) k  m1 kε(0) k . Ha kAm k < 1, akkor 1 m σ ≤ (kAm k) m = e−R(A ) , azaz m ))−1 σ (R(A 12 1 ≤ , e http://www.doksihu tehát Nm ≡ (R(Am ))−1 megadja, hogy hány iterációs lépés kell ahhoz, hogy a hiba a kezdeti hiba e-adnyi részére csökkenjen. El®fordulhat azonban, hogy egy ε vektorra kBmε k < kAmε k, noha kAm k < kBm k. Példa: Legyenek A, B és ε a következ®k: " A= 1 2 0 0 1 2 # " , B= 3 4 0 0 1 100 # " ε= , 0 # . 1 Ekkor ε kezdeti hiba esetén " Amε = # 0 " , ( 12 )m illetve Bmε = 0 1 m ( 100 ) # , vagyis a B mátrixú iterációnál ε m kisebb lesz, noha

deníció szerint A gyorsabb m lépésben. Továbbá az is el®fordulhat, hogy m1 iteráció esetén az A mátrix gyorsabb a B mátrixnál, de m2 iterációnál már B a gyorsabb. Példa: Legyenek A és B mátrixok a következ®k: " A= 1 4 4 0 1 4 # " , B= 1 4 0 0 1 2 # . Ekkor m1 = 6 iteráció esetén 0.023 = kA6 k > kB6 k = 0015, azonban m2 = 7 iterációnál 6.8 · 10−3 = kA7 k < kB7 k = 78 · 10−3 3.1 Megjegyzés Van azonban olyan mátrixosztály, melyre teljesül kAm k monotonitása Ha A1 és A2 hermitikus (vagy normális) mátrixok (ld. Függelék), akkor m kAm i k = (ρ(Ai )) , ezért, ha ρ(A1 ) < ρ(A2 ) < 1, akkor ∀ m = 1, 2, .-re m kAm 1 k < kA2 k < 1, Az aszimptotikus és az átlagos sebesség kapcsolatát mutatja a következ® tétel. 3.4 Tétel Legyen A egy konvergens n × n-es komplex mátrix Ekkor lim R(Am ) = R∞ (A). m∞ 13 http://www.doksihu 3.3 Bizonyítás R(Am ) és R∞ (A) deníciójából, valamint

a Gelfand-formulából adódik 3.2 Következmény Mivel kAm k ≥ (ρ(A))m , ha A egy konvergens mátrix, akkor R∞ (A) ≥ R(Am ), minden olyan pozitív egész m-re, melyre kAm k < 1. El®fordulhat, hogy R(Am ) nagyon lassan konvergál R∞ (A)-hoz, és ilyenkor az aszimptotikus sebesség félrevezet® információt adhat nekünk. Tekintsük a következ® példát: Legyen A mátrix " # A= 0.99 4 0 0.99 alakú. Ekkor R∞ (A) = 0.01005 és N∞ = 995, ami azt sugallja nekünk, hogy 100 iterációs lépés már elég ahhoz, hogy a kezdeti hibát e-adnyi részére csökkentsük. Azonban kAm k kezdetben növekszik, és m < 805 esetén kAm k ≥ 1. Ahhoz, hogy kAm k ≤ 1e legyen, legalább 918 iterációra van szükség. 14 http://www.doksihu 3.2 A Richardson-, a Jacobi-, a Gauss-Seidel és a SORmódszer A Richardson-iteráció Az egyik legkorábbi és legkisebb m¶veletigény¶ módszer a Richardson-iteráció, hiszen itt még inverzet sem kell számolnunk. Az

Ax = k egyenlet mindkét oldalához hozzáadva (I − A)x-et, a kapott módszer a következ®: x(m+1) = (I − A)x(m) + k, m = 0, 1, 2. (Valójában tekinthet® úgy, hogy az A = M − N felbontásban M = I és N = I − A.) Az iteráció konvergens, ha ρ(I − A) < 1. A módszert javíthatjuk, ha bevezetünk egy paramétert: x(m+1) = (I − pA)x(m) + pk, m = 0, 1, 2., és ekkor p megfelel® megválasztásával biztosíthatjuk, s®t gyorsíthatjuk is a konvergenciát. 3.1 Állítás A Richardson-iteráció szimmetrikus, pozitív denit mátrixokra konvergens, ha p < 2 λmax , ahol λmax az A mátrix maximális sajátértéke. Tekintsük az A = M−N mátrix néhány olyan klasszikus, egyszer¶ felbontását, amikor M kiszámítása (speciális alakjának köszönhet®en) könny¶. A továbbiakban mindig feltesszük, hogy A reguláris mátrix, s®t, A minden diagonális eleme nemnulla. −1 A Jacobi-módszer Vegyük A mátrixnak azt a felbontását, amikor M = D =

diag(a1,1 , a2,2 , . , an,n ) és N = E + F, ahol E és F szigorú alsó ill. fels® háromszögmátrixok, melyek elemei A megfelel® elemeinek ellentettjei: A = D − E − F. Behelyettesítve ezt az Ax = k egyenletbe kapjuk, hogy Dx = (E + F)x + k, 15 http://www.doksihu vagyis (m+1) ai,i xi =− n X (m) ai,j xj + ki , i = 1, ., n, m = 0, 1, 2., j=1 j6=i ahol x(m) i -mel jelöljük az m-edik iterációs vektor i. komponensét Mivel feltettük, hogy A mátrix diagonális elemei nemnullák, ezért (m+1) xi =− n  X ai,j  j=1 j6=i ai,i (m) xj + ki , ai,i i = 1, ., n, m = 0, 1, 2. Mátrixos alakban pedig az el®z®ekkel ekvivalens egyenletek Dx(m+1) = (E + F)x(m) + k, m = 0, 1, 2., illetve x(m+1) = D−1 (E + F)x(m) + D−1 k, (3.4) m = 0, 1, 2., ahol B = D−1 (E + F)-et Jacobi-mátrixnak nevezzük. Mivel M itt diagonális mátrix, inverze könnyen kiszámolható, így x(m+1) gyorsan el®állítható x(m) -b®l. A Gauss-Seidel módszer A Gauss-Seidel

módszer egyik el®nye a Jacobival szemben, hogy x(m+1) kiszámításához nem kell eltárolni x(m) összes elemét, minden komponensnek csak egy iterációját tároljuk egyszerre. (m+1) ai,i xi =− i−1 X (m+1) ai,j xj j=1 − n X (m) ai,j xj + ki , i = 1, ., n, m = 0, 1, 2. j=i+1 Tehát az A = D − E − F felbontásnál most M = D − E és N = F, azaz az iteráció: (D − E)x(m+1) = Fx(m) + k, m = 0, 1, 2. Mivel a D − E mátrix reguláris, x(m+1) = (D − E)−1 Fx(m) + (D − E)−1 k, ahol (D − E)−1 F a Gauss-Seidel mátrix. 16 m = 0, 1, 2., (3.5) http://www.doksihu Ennél a módszernél M−1 = (D − E)−1 kiszámolása már kicsit nehezebb, mint a Jacobinál, hiszen itt egy alsó háromszögmátrixról van szó, azonban, mint kés®bb látni fogjuk, a vektorsorozat gyorsabban konvergál x-hez. A SOR-módszer A SOR-módszer deniálásához el®ször bevezetünk egy segédvektort: (m+1) ai,i x̃i =− i−1 X (m+1) ai,j xj j=1 − n X

(m) ai,j xj + ki , i = 1, ., n, m = 0, 1, 2. (3.6) j=i+1 Az (m + 1). iteráció i komponensét egy súlyozott átlagként kapjuk meg: (m+1) xi (m) = (1 − ω)xi (m+1) + ωx̃i (3.7) , ahol ω ≥ 0 az iterációs paraméter (súly). Az 36 és 37 egyenletekb®l (m+1) ai,i xi (m) = ai,i xi n i−1 o n X X (m) (m) (m+1) , − ai,j xj + ki − ai,i xi +ω − ai,j xj j=i+1 j=1 illetve (D − ωE)x(m+1) = {(1 − ω)D + ωF}x(m) + ωk, m = 0, 1, 2., továbbá, mivel (D − ωE) reguláris mátrix, x(m+1) = (I − ωL)−1 {(1 − ω)I + ωU}x(m) + ω(I − ωL)−1 D−1 k, ahol L ≡ D−1 E, U ≡ D−1 F , (I − ωL)−1 {(1 − ω)I + ωU} pedig a SOR-módszer mátrixa. A Gauss-Seidel az ω = 1 speciális esete a SOR-módszernek. Azonban az ω paramétert jól megválasztva gyorsabb konvergenciát is elérhetünk a Gauss-Seidelnél. A továbbiakban megvizsgáljuk, hogy a fent deniált iterációs módszerek milyen mátrixokra lesznek konvergensek, és

összehasonlítjuk a konvergenciájuk gyorsaságát. 17 http://www.doksihu 3.5 Tétel (A Stein-Rosenberg tétel) Legyen a Jacobi-mátrix B ≡ L + U alakú nemnegatív, n × n-es mátrix 0 diagonális elemekkel, G = (I−L)−1 U pedig a Gauss-Seidel mátrix. Ekkor a következ® állítások közül mindig pontosan egy teljesül: 1. ρ(B) = ρ(G) = 0 2. 0 < ρ(G) < ρ(B) < 1 3. 1 = ρ(B) = ρ(G) 4. 1 < ρ(B) < ρ(G) Tehát a Jacobi- és a Gauss-Seidel mátrixok mindegyike egyszerre konvergens vagy divergens adott A = D − L − U mátrix esetén. 3.3 Következmény Ha egy nemnegatív Jacobi-mátrix spektrálsugarára 0 < ρ(B) < 1, akkor R∞ (G) > R∞ (B). 3.6 Tétel Legyen A egy szigorúan vagy irreducibilisen diagonálisan domináns (ld Függelék) n × n-es mátrix Ekkor a megfelel® Jacobi- és Gauss-Seidel mátrix is konvergens, tehát a 3.4 és 35 iterációs módszerek az Ax = k egyenletrendszerre konvergensek lesznek bármely x(0) kezdeti

vektor esetén. 3.4 Bizonyítás A 66 denícióból és 62 tételb®l tudjuk, hogy A reguláris és diagonális elemei nemnullák. A Jacobi-mátrix bi,j elemei a következ®képpen írhatók fel A komponenseinek segítségével: ( bi,j = 0, a i=j − ai,j , i 6= j i,i ) . A 6.6 denícióból következik, hogy szigorúan diagonálisan domináns A esetén n X i = 1, ., n-re, |bi,j | < 1 j=1 irreducibilisen diagonálisan domináns A mátrix esetén pedig n X i = 1, ., n-re, |bi,j | ≤ 1 j=1 és legalább egy i-re szigorú egyenl®tlenség van. Ekkor |B| = (|bi,j |) mátrixra teljesül, hogy ρ(|B|) < 1. A 65 tételt alkalmazva kapjuk, hogy ρ(B) ≤ ρ(|B|) < 1, vagyis a Jacobi-módszer konvergens, és a Stein-Rosenberg tételb®l következik, hogy ekkor a Gauss-Seidel módszer is konvergens. 18 http://www.doksihu Szükséges és elégséges feltételek a SOR-módszer (és Gauss-Seidel módszer) konvergenciájához 3.7 Tétel (Kahan tétele) Legyen B ≡ L+U

egy n×n-es komplex mátrix 0 diagonális elemekkel, ahol L és U szigorú alsó és fels® háromszögmátrixok. Ha Lω ≡ (I − ωL)−1 {ωU + (1 − ω)I}, akkor minden valós ω -ra ρ(Lω ) ≥ |ω − 1|, ahol egyenl®ség csak akkor teljesül, ha Lω minden valós sajátértéke |ω − 1|. Kahan tételéb®l következik, hogy a SOR-módszer csak 0 < ω < 2 esetén lehet konvergens. Nézzük azt az esetet, amikor az Ax = k egyenletrendszerben A egy n × n-es hermitikus mátrix és a következ® alakban írható fel A = D − E − E∗ , ahol D és E n×n-es mátrixok és D pozitív denit hermitikus mátrix. Továbbá feltesszük, hogy D − ωE reguláris minden 0 ≤ ω ≤ 2 esetén. Ekkor az Ax = k egyenletrendszert átírva kapjuk a következ® iterációs módszert: (D − ωE)x(m+1) = {ωE∗ + (1 − ω)D}x(m) + ωk, m = 0, 1, 2. Bevezetve az L ≡ D−1 E és U ≡ D−1 E∗ jelöléseket kapjuk, hogy x(m+1) = Lω x(m) + ω(D − ωE)−1 k, m = 0, 1,

2., ahol Lω = (I − ωL)−1 {ωU + (1 − ω)I}, hasonlóan a SOR-módszerhez, azonban itt L és U nem szigorú alsó és fels® háromszögmátrixok. Erre az iterációs módszerre vonatkozik a következ® tétel: 3.8 Tétel (Ostrowski tétele) Legyen A = D − E − E∗ egy n × n-es hermitikus mátrix, ahol D pozitív denit hermitikus mátrix és D − ωE reguláris minden 0 ≤ ω ≤ 2 esetén. Ekkor ρ(Lω ) < 1 akkor és csak akkor, ha A pozitív denit és 0 < ω < 2. 19 http://www.doksihu 3.4 Következmény Legyen A = D − E − E∗ egy n × n-es hermitikus mátrix, ahol D pozitív denit hermitikus mátrix és D − E reguláris. Ekkor a Gauss-Seidel módszer konvergens akkor és csak akkor, ha A pozitív denit. Az ebben a fejezetben leírtak, az itt nem szerepl® bizonyításokkal együtt, b®vebben megtalálhatóak Richard S. Varga: Matrix Iterative Analysis [2] cím¶ könyvének 3 fejezetében 20 http://www.doksihu 4. fejezet Nemnegatív

mátrix szeletelés A nemnegatív mátrix szeletelés elmélete a Perron-Frobenius tételeken alapul, és a Richard S. Vargától származó reguláris felbontásnak egy általánosítása, melyet kés®bb tovább fejlesztettek a prefaktorizációs módszerek (AGA, EWA) kialakításánál. Az el®z® fejezetben tárgyalt módszerekt®l eltér®en, most az egyenletrendszer együtthatómátrixának felbontását nem határozzuk meg egyértelm¶en, azonban néhány, nemnegativitásra vonatkozó kikötést teszünk A, M és N mátrixokra, melyekkel biztosítjuk az iteráció konvergenciáját. 4.1 Mátrix szeletelés A továbbiakban A = (ai,j ) ≥ 0 mindig nemnegativitást fog jelenteni, vagyis minden ai,j ≥ 0 és legalább egy ai,j > 0. A = (ai,j ) > 0 pozitív mátrix, ha minden ai,j > 0 A Perron-Frobenius elmélet sok összehasonlító tételt ad a nemnegatív mátrixok sajátértékeire és sajátvektoraira vonatkozóan, ezek közül a legfontosabbak: 4.1 Tétel (I) Ha

A ≥ 0 mátrix, akkor 1. A-nak van a spektrálsugarával egyenl®, nemnegatív valós sajátértéke 2. ρ(A) > 0-hoz tartozik egy x ≥ 0 sajátvektor 3. ha az A mátrix bármelyik elemét növeljük, ρ(A) nem csökken 4.2 Tétel (II) Ha A ≥ 0 irreducibilis mátrix, akkor 1. A-nak van a spektrálsugarával egyenl®, pozitív valós sajátértéke 2. ρ(A) > 0-hoz tartozik egy x > 0 sajátvektor 21 http://www.doksihu 3. ha az A mátrix bármelyik elemét növeljük, ρ(A) is n® 4. ρ(A) az A mátrixnak egyszeres sajátértéke 4.1 Deníció Az A mátrixnak A = M−N konvergens felbontása, ha A és M reguláris mátrixok és ρ(M−1 N) < 1. 4.3 Tétel Legyen A-nak A = M − N egy felbontása Ha A és M reguláris mátrixok, akkor M−1 NA−1 = A−1 NM−1 , (4.1) továbbá M−1 N és A−1 N, illetve NA−1 és NM−1 kommutáló mátrixok. 4.1 Bizonyítás A = M − N-b®l kapjuk, hogy M−1 = (A + N)−1 = A−1 (I + NA−1 )−1 = (I + A−1 N)−1

A−1 , tehát A−1 = M−1 + M−1 NA−1 = M−1 + A−1 NM−1 , amib®l következik a tétel. A 41 egyenletet balról és jobbról megszorozva N-nel kapjuk, hogy az M−1 N és A−1 N, illetve NA−1 és NM−1 mátrixok kommutálnak. 4.2 Deníció Egy A mátrix monoton, ha reguláris és A−1 ≥ 0 Deniáljuk egy monoton A mátrixnak néhány felbontását: 4.3 Deníció Az A mátrixnak az A = M − N reguláris felbontása, ha M reguláris, M−1 ≥ 0 és N ≥ 0. 4.4 Deníció Az A mátrixnak az A = M − N nemnegatív felbontása, ha M reguláris, M−1 ≥ 0 , M−1 N ≥ 0 és NM−1 ≥ 0. 4.5 Deníció Az A mátrixnak az A = M − N gyengén nemnegatív felbontása, ha M reguláris, M−1 ≥ 0 és M−1 N ≥ 0 vagy NM−1 ≥ 0. A reguláris felbontás denícióját Varga vezette be. A 44 deníció ekvivalens a más szerz®k által gyakran használt gyengén reguláris felbontással. Ez utóbbit úgy is deniálják, hogy M−1 ≥ 0 és M−1 N ≥ 0

(az NM−1 ≥ 0 feltétel nélkül), azonban ebben az esetben további feltevések szükségesek az összehasonlító tételekben. 22 http://www.doksihu Triviálisan adódik az el®bbi deníciókból az alábbi következmény. 4.1 Következmény Egy A mátrix minden reguláris felbontása egyben nemnegatív felbontása is A-nak és minden nemnegatív felbontása egyben gyengén nemnegatív felbontása is. Vagyis az alábbi, gyengén nemnegatív felbontásokra kimondott tételek mindegyike érvényes a reguláris felbontásokra is. 4.4 Tétel Legyen az A mátrixnak A = M − N egy gyengén nemnegatív felbontása Ha A−1 ≥ 0, akkor 1. A−1 ≥ M−1 , 2. ρ(M−1 N) = ρ(NM−1 ) < 1, 3. ha M−1 N ≥ 0, akkor A−1 N ≥ M−1 N, illetve ha NM−1 ≥ 0, akkor NA−1 ≥ NM−1 , 4. ρ(M−1 N) = ρ(A−1 N) . 1 + ρ(A−1 N) (4.2) 5. Fordítva is igaz, ha ρ(M−1 N) < 1, akkor A−1 ≥ 0 4.2 Bizonyítás (1) A 43 tételb®l kapjuk, hogy A−1 = M−1 + M−1

NA−1 = M−1 + A−1 NM−1 , és mivel feltettük, hogy M−1 N ≥ 0 vagy NM−1 ≥ 0, ezért M−1 NA−1 = A−1 NM−1 ≥ 0, vagyis A−1 ≥ M−1 . (2) Tegyük fel, hogy M−1 N ≥ 0, ekkor A−1 = M−1 + M−1 NA−1 = M−1 + M−1 N(M−1 + M−1 NA−1 ) = = [I + M−1 N]M−1 + (M−1 N)2 A−1 = [I + M−1 N + (M−1 N)2 ]M−1 + (M−1 N)3 A−1 = = [I + M−1 N + (M−1 N)2 + . + (M−1 N)k−1 ]M−1 + (M−1 N)k A−1 Mivel A−1 , M−1 és M−1 N nemnegatív mátrixok, az I + M−1 N + (M−1 N)2 + . 23 (4.3) http://www.doksihu sorozat konvergens, és a 6.6 tételb®l ρ(M−1 N) < 1 Hasonlóan kapjuk az NM−1 ≥ 0 esetre, hogy I + NM−1 + (NM−1 )2 + . konvergens. Felhasználva a 64 tételt, ρ(M−1 N) = ρ(NM−1 ) < 1 (3) A 4.3 egyenletben k ∞ esetén (M−1 N)k 0, és a 66 tételb®l I + M−1 N + (M−1 N)2 + . = (I − M−1 N)−1 ≥ I ≥ 0 Továbbá, mivel A−1 = (I − M−1 N)−1 M−1 , (4.4) A−1 N = (I − M−1

N)−1 M−1 N ≥ M−1 N. (4.5) rögtön adódik, hogy A másik eset (NM−1 ≥ 0) hasonlóan levezethet®. (4) Mivel M−1 N és A−1 N kommutáló mátrixok, azonosak a sajátvektoraik, vagyis M−1 Nvi = λi vi , A−1 Nvi = τi vi minden vi sajátvektorra. A 45 egyenletb®l (I − M−1 N)−1 M−1 Nvi = τ vi , azaz M−1 Nvi = τi vi . 1 + τi Ebb®l a következ® összefüggést kapjuk a sajátértékekre: λi = τi . 1 + τi Az M−1 N ≥ 0 és A−1 N ≥ 0 feltevésekb®l, valamint az I. Perron-Frobenius tételb®l a megfelel® spektrálsugarak kapcsolata: ρ(M−1 N) = ρ(A−1 N) . 1 + ρ(A−1 N) (5) M−1 N ≥ 0 és ρ(M−1 N) < 1 esetén a 6.6 tételb®l adódik, hogy (I − M−1 N)−1 ≥ 0, és 4.4 egyenletb®l közvetlenül kapjuk, hogy ekkor A−1 ≥ 0 ( Az NM−1 ≥ 0 és ρ(NM−1 ) < 1 esetre ugyanígy belátható.) 24 http://www.doksihu 4.2 Következmény Egy monoton A mátrix minden gyengén nemnegatív felbontása konvergens

felbontás, és fordítva, ha A-nak létezik gyengén nemnegatív konvergens felbontása, akkor A−1 ≥ 0. A továbbiakban egy mátrix különböz®, gyengén nemnegatív felbontásai által kapott iterációk gyorsaságát hasonlítjuk össze, melyhez bevezetünk egy segédtételt. 4.5 Tétel Legyen A = M1 − N1 = M2 − N2 az A mátrix két gyengén nemnegatív felbontása és A−1 ≥ 0. Ekkor, ha a következ® egyenl®tlenségek bármelyike teljesül −1 1. A−1 N2 ≥ A−1 N1 ≥ 0 (vagy M−1 2 N2 ≥ M1 N1 ≥ 0) −1 2. A−1 N2 ≥ N1 A−1 ≥ 0 (vagy M−1 2 N2 ≥ N1 M1 ≥ 0) −1 3. N2 A−1 ≥ N1 A−1 ≥ 0 (vagy N2 M−1 2 ≥ N1 M1 ≥ 0) −1 4. N2 A−1 ≥ A−1 N1 ≥ 0 (vagy N2 M−1 2 ≥ M1 N1 ≥ 0) , akkor −1 ρ(M−1 1 N1 ) ≤ ρ(M2 N2 ) < 1. (4.6) 4.3 Bizonyítás A 64 és 65 tételb®l, valamint a 44 tétel 4 állításából következik 4.6 A zárójelekben megadott feltételek esetén pedig a 64 és a 65 tételb®l rögtön adódik az

állítás. El®ször hasonlítsuk össze a felbontásokat az N mátrix szerint: 4.6 Tétel Legyen A = M1 − N1 = M2 − N2 az A mátrix két azonos típusú, gyengén −1 −1 nemnegatív felbontása, azaz vagy M−1 ≥ 0 és 1 N1 ≥ 0 és M2 N2 ≥ 0 vagy N1 M1 −1 N2 M−1 ≥ 0. Ha N2 ≥ N1 , akkor 2 ≥ 0, továbbá legyen A −1 ρ(M−1 1 N1 ) ≤ ρ(M2 N2 ) < 1. (4.7) 4.4 Bizonyítás Az N2 ≥ N1 egyenl®tlenségb®l következik, hogy vagy A−1 N2 ≥ A−1 N1 ≥ 0 vagy N2 A−1 ≥ N1 A−1 ≥ 0, és a 4.5 tételb®l adódik 47 4.1 Megjegyzés Evidens, hogy az el®bbi két tétel mindegyike érvényes A nemnegatív felbontásaira is. Az iterációs módszereknél nem mindig tudjuk az N mátrixokat összehasonlítani (vagyis sem N1 ≥ N2 , sem N2 ≥ N1 nem teljesül, kivéve pl. a Jacobi- és a Gauss-Seidel módszert), az M−1 -eket azonban gyakrabban Azt várhatjuk, hogy minél közelebb van M az A-hoz, annál gyorsabb lesz az iteráció. A

következ®kben megnézzük, hogy M−1 hogyan befolyásolja ρ(M−1 N)-et. 25 http://www.doksihu 4.7 Tétel Legyen A = M1 − N1 = M2 − N2 az A mátrix két nemnegatív felbontása és −1 A−1 ≥ 0. Ha M−1 1 ≥ M2 , akkor −1 ρ(M−1 1 N1 ) ≤ ρ(M2 N2 ) < 1. (4.8) 4.5 Bizonyítás A Perron-Frobenius, a 44 és a 64 tételekb®l λ1 = ρ(N1 M−1 1 ) = −1 = ρ(M−1 1 N1 ) < 1 és λ2 = ρ(M2 N2 ) < 1, a hozzájuk tartozó sajátvektorok, x1 és x2 pedig nemnegatívak. Ezért N1 M−1 1 x1 = λ1 x1 ≥ 0 (4.9) és T xT2 M−1 2 N2 = λ2 x2 ≥ 0. (4.10) Szorozzuk meg a 4.9 egyenletet balról, a 410 egyenletet pedig jobbról A−1 -zel: −1 A−1 N1 M−1 1 x 1 = λ1 A x 1 , (4.11) −1 xT2 M−1 = λ2 xT2 A−1 . 2 N2 A (4.12) −1 −1 −1 −1 A−1 = M−1 = M−1 2 + M 2 N2 A 1 + A N1 M 1 , (4.13) Tudjuk, hogy melyet átírva kapjuk, hogy −1 −1 −1 M−1 − A−1 N1 M−1 2 N2 A 1 = M1 − M2 ≥ 0, (4.14) −1 M−1 ≥ A−1 N1

M−1 2 N2 A 1 ≥ 0, (4.15) vagyis és a nemnegativitás deníciójából következik, hogy ekkor −1 −1 −1 −1 M−1 2 N2 A x1 ≥ A N1 M1 x1 = λ1 A x1 . (4.16) Szorozzuk meg a 4.16 egyenletet balról xT2 -tal, a 412 egyenletet pedig jobbról x1 -gyel: −1 T −1 xT2 M−1 2 N2 A x 1 ≥ λ1 x2 A x1 , (4.17) −1 T −1 xT2 M−1 2 N2 A x 1 = λ2 x2 A x 1 , (4.18) λ1 xT2 A−1 x1 ≤ λ2 xT2 A−1 x1 . (4.19) melyekb®l adódik, hogy Mivel xT2 A−1 x1 > 0, ebb®l következik, hogy λ1 ≤ λ2 , amivel beláttuk az állítást. 26 http://www.doksihu −1 4.2 Megjegyzés Az el®bbi tételnél er®sebb állítás is igaz: ha M−1 1 > M2 , akkor 4.8-ban szigorú egyenl®tlenség teljesül A bizonyítás hasonlóan levezethet® 4.1 Állítás Legyen A = M1 − N1 = M2 − N2 az A mátrix két gyengén nemnegatív −1 felbontása. Ha N2 ≥ N1 , akkor M−1 1 ≥ M2 ≥ 0. −1 Az állítás megfordítása azonban nem igaz, ezért az M−1 1 ≥ M2

feltétel gyengébb megkötés az N2 ≥ N1 -nél. Gyengén nemnegatív felbontások esetén M−1 -ek összehasonlítására az alábbi tételek vonatkoznak: 4.8 Tétel Legyen A = M1 − N1 = M2 − N2 az A monoton mátrix két különböz® típusú −1 −1 gyengén nemnegatív felbontása, azaz vagy M−1 1 N1 ≥ 0 és N2 M2 ≥ 0, vagy N1 M1 ≥ 0 −1 −1 és M−1 2 N2 ≥ 0. Ha M1 ≥ M2 , akkor −1 ρ(M−1 1 N1 ) ≤ ρ(M2 N2 ) < 1. (4.20) −1 4.6 Bizonyítás Tekintsük azt az esetet, amikor M−1 1 N1 ≥ 0 és N2 M2 ≥ 0. Ekkor T y1T M−1 1 N1 = λ1 y1 (4.21) N2 M−1 2 y 2 = λ2 y2 . (4.22) és A Perron-Frobenius, 4.2 és a 64 tételb®l kapjuk, hogy λ1 = ρ(M−1 1 N1 ) < 1 és λ2 = −1 −1 ρ(M2 N2 ) = ρ(N2 M2 ) < 1, valamint a megfelel® sajátvektorok, y1 ill. y2 nemnegatívak Megszorozva a 4.21 egyenletet jobbról A−1 -zel, a 422 egyenletet pedig balról A−1 -zel azt kapjuk, hogy −1 = λ1 y1T A−1 (4.23) y1T M−1 1 N1 A és −1

A−1 N2 M−1 2 y 2 = λ2 A y 2 . (4.24) −1 A 4.13 és 41 egyenletekb®l, valamint az M−1 1 ≥ M2 feltevésb®l következik, hogy −1 −1 ≥0 A−1 N2 M−1 2 ≥ M 1 N1 A (4.25) T −1 −1 = λ1 y1T A−1 . y1T A−1 N2 M−1 2 ≥ y1 M1 N1 A (4.26) és 27 http://www.doksihu Megszorozva 4.26 egyenletet jobbról y2 -vel, 424 egyenletet pedig balról y1T -tal, azt kapjuk, hogy T −1 y1T A−1 N2 M−1 (4.27) 2 y2 ≥ λ1 y1 A y2 és (4.28) T −1 y1T A−1 N2 M−1 2 y 2 = λ2 y1 A y2 . −1 Mivel y1T A−1 y2 > 0, λ1 ≤ λ2 , amivel beláttuk az állítást. Az N1 M−1 1 ≥ 0 és M2 N2 ≥ 0 eset ugyanígy bizonyítható. −1 4.3 Megjegyzés Az el®bbi tételnél er®sebb állítás is igaz: ha A−1 > 0 és M−1 1 > M2 , −1 akkor ρ(M−1 1 N1 ) < ρ(M2 N2 ) < 1. Az el®bbi tétel természetesen arra az esetre is igaz, amikor A egy nemnegatív és egy gyengén nemnegatív felbontását hasonlítjuk össze. Azonban két azonos típusú

gyengén nemnegatív felbontás esetén az M−1 -ek összehasonlítása nem ad információt a konvergenciák sebességér®l. Ezt szemlélteti a következ® példa. Tekintsük az alábbi monoton A mátrixot:  1  A=  0 −1 −1 1 0 0   −1  , 2  2 2 1   melyre A−1 =   1 2 1 ,   1 1 1 és annak néhány felbontását: 28 (4.29) 1 1 29 1 −1  − 14 23 4 4  − 23   3  0  −1 19 23 9 23 −1 4 5 2 5 −1 0 1 2  0  −1 − 56 5 6 1 3 −1 1   1 1  − 12  M   45 23    − 20 23  0 12 5   − 65   0  − 32   3 0  −1   2 0 0 0  0 0 0  0 0 0  9 23 0 9 23 0 3 23 1 − 23 0    0 −1 −1   5 5  2 0 25 5     0 0 −1   2  0 0 1    −4 −4  23 23  0    − 1 −1 0   12  6 1 1 0

6 3  N     3 4 2 3     9 5 4 5 2 5 2 3 9 5 9 5 1 5 1 2 4 5 4 5 3 5 2 3 2 2 1 2 3          1 2 1     2 3 2 2 1 1 2     1 2 1     2 3 2 2 1    1 2 1     M−1 0 0 0   0 0 0 M−1 N  0 0 2 3 0 0 2 3 0   0 0 0     −2 −3 −1   5 5 5  7 3 1 5 5     −1 −1 −1   3 2 3  2 2 1 3 3  2 3 0    −1 −1 −1   3 3 3   0 0 0  0 0 0  1 5 1 5 0 0   0 0   0 1 5 1 5        0 0 0    0 16 16     0 0 0    0 0 31       −1 −1 −1   0 0 0   3  2 4   1 1 1 2 1 0 3 2 12 6  NM−1 − 2 5 √ 1 6 1 3 0 ρ(M−1 N) http://www.doksihu 30 −1 1 2 8  0  −2 1 2  −2 2 7  0  −2 1

2  −2 2 6  0  −2 1 2  − 32 1 5  0  −2 2  M      −2   4 0  −2   5 0  −2   4 0  −1   2 0 1 1 1 1 1 −1 1 1 −1 1 0 1 1    −1   2 0  −1   3          0     −1   2 0        0   0 −1 0 1 −1   0  0    0  −1    0  −1    0  − 12  N 5 8 1 8 1 8 3 4 1 4 1 4 5 6 1 3 1 3 1 2 1 2 1 0 1 2 1 2 5 8 5 8 1 8 2 3 2 3 1 6 4 3 4 3 1 3 M−1 1 4 1 4 1 4 1 4 1 4 1 4 1 3 1 3 1 3 2 3 2 3 2 3 1 2 0 0   1 2 3 8     0     1 4 1 2     0     1 6 1 2     0     2 3 3 4 1 2 1 2  3 4  0     0   0 0 1 2 3 8   0   0 0 1 2 1 3 0 0 2 3 1 6 1 6 1 2 1

2     5 8 1 8 1 8   0  0  1 8 5 8 1 8 1 8 5 8 1 8 1 6 2 3 1 6 1 6 1 6  1 2       0   0 1 8 1 8 5 8 1 2   0   0  0   0 − 13 0 M−1 N   0  0        0 0 0       0 32 13   NM−1 3 4 3 4 2 3 1 2 ρ(M−1 N) http://www.doksihu 1 0 31 1 0 0     12  0 1 0   0 0 2     −1   2 0 0 0 0 1   0 0 0 N   0 0 0  0 1 0 1 1 0 1 1   1 2 1 2 1 2 1 2 1 2 1 2 1 0 0   0 1  0 0  3 2 1 2 1 2 M−1 0 0 0  0 0 0  0 1 0     0 0 0    1 1 12             0 0 0     1 1 1 2 2   NM−1 0 1 0     0 0 1    1 0 0  1 0 0     0 1 0    1 0 0 2  0 1 0     0 0 1   2  1 0 0     

   0 0 1   0 1 0   1 0 1       2 2  1 1 0 0 0 0 0 0 0 2 2     0 0 0    1 0 0       0 0 0   −1      0 1 0 2 0    11  0 1 0   −1 0 2   1 1 −1 10   0  −1 1 9  0  −1 1  M 0     0 0   0 0 0 0  0   0 0 0 1 0 1 2 1 2 1 2 1 2 1 2 1 2 0 1 0     0 0 1    1 0 0 2     0 0 1    0 12 0        0  0  M−1 N 3 q q 1 2 1 2 1 2 1 2 ρ(M−1 N) http://www.doksihu http://www.doksihu A táblázatban meggyelhet®, hogy az els® négy felbontás azonos típusú (M−1 N ≥ 0), míg az ötödik a másik típusú (NM−1 ≥ 0) gyengén nemnegatív felbontás, a következ® három nemnegatív, az utolsó négy pedig reguláris mátrix szeletelés. Látható, hogy azonos −1 −1 −1 −1 −1

felbontásoknál M−1 1 ≥ M3 > 0 és ρ(M1 N1 ) < ρ(M3 N3 ), ugyanakkor M2 ≥ M3 > −1 0 és ρ(M−1 2 N2 ) > ρ(M3 N3 ). A korábbi tételekre is példák a fenti szeletelések: −1 −1 −1 - M−1 1,2,3,4,5 > M6 és ρ(M1,2,3,4,5 N1,2,3,4,5 ) < ρ(M6 N6 ) egy nemnegatív és egy gyengén nemnegatív felbontásra és szigorú egyenl®tlenségre, −1 −1 −1 -M−1 1,2,3 > M5 és ρ(M1,2,3 N1,2,3 ) < ρ(M5 N5 ) a különböz® típusú gyengén nemnegatív felbontásokra és szigorú egyenl®tlenségre, −1 −1 −1 -M−1 7 ≥ M8 és ρ(M7 N7 ) ≤ ρ(M8 N8 ) két nemnegatív felbontásra, −1 −1 −1 -M−1 6 > M7,8 és ρ(M6 N6 ) < ρ(M7,8 N7,8 ) két nemnegatív felbontásra és szigorú egyenl®tlenségre. Azonos típusú gyengén nemnegatív felbontások is összehasonlíthatók M−1 alapján, de csak akkor, ha feltesszük, hogy A és M1 vagy M2 szimmetrikus mátrix. 4.9 Tétel Legyen A = M1 −N1 = M2 −N2 egy szimmetrikus és

monoton A mátrix két −1 gyengén nemnegatív felbontása. Ha M−1 1 ≥ M2 , és M1 vagy M2 szimmetrikus, akkor −1 ρ(M−1 1 N1 ) ≤ ρ(M2 N2 ) < 1. Természetesen még sok összehasonlító tétel létezik gyengén nemnegatív és ezáltal nemnegatív és reguláris felbontásokra is, melyek például felhasználják A, M, N transzponáltjait is, és a különböz® mátrixok ill. mátrixszorzatok közötti relációt vizsgálják (Mindezek, valamint az itt nem szerepl® bizonyítások is részletesen megtalálhatók Woznicki: Nonnegative Splitting Theory [1] cím¶ jegyzetének 3. fejezetében) 32 http://www.doksihu 4.2 Prefaktorizációs iterációs módszerek A prefaktorizációs iterációs módszerekben a nemnegatív mátrix szeletelés elméletének eredményeit fejlesztették tovább. Vegyük ismét az Ax = k lineáris algebrai egyenletrendszert. Ha A = M − N az A-nak egy felbontása, akkor az iterációs módszerünk x(m+1) = M−1 Nx(m) + M−1 k, m = 0,

1, 2, . (4.30) alakú, és konvergens, ha ρ(M−1 N) < 1. M-et úgy is tekinthetjük, mint A egy közelítését, melyet általában reguláris mátrixok szorzataként állítanak el® úgy, hogy könnyen kiszámolható és invertálható legyen. M-et A prefaktorizációjának is nevezik, N = M − A-t pedig a különbségmátrixnak (residual matrix ). Ha N a nullmátrix, akkor az egyenletrendszer megoldását közvetlenül a Gausseliminációs direkt módszerrel kapjuk meg (strict factorization) Ha N nem a nullmátrix, akkor pedig a megoldást iterálva érjük el (partial factorization). Deniáljuk egy reguláris A mátrix felbontását a következ®képpen: A = D − E − F, (4.31) ahol D egy reguláris és diagonális, E, F pedig szigorú alsó ill. fels® háromszögmátrixok Bevezetünk további szigorú alsó (H) és fels® (Q) háromszögmátrixokat, melyek segítségével a faktorizációs módszer M mátrixa így írható fel: M = [I − (E + H)K−1 ]K[I − K−1

(F + Q)], (4.32) ahol K = D − diag{(E + H)K−1 (F + Q)} egy reguláris, diagonális mátrix, N pedig N = odiag{(E + H)K−1 (F + Q)} − H − Q. (4.33) A fenti egyenletekben diag{C} olyan diagonális mátrixot jelöl, melynek a f®átlóbeli elemei megegyeznek C diagonális elemeivel, a többi eleme pedig 0, és odiag{C} = C-diag{C}. Ekkor a 4.30 módszer iterációs mátrixa: F ≡ M−1 N = [I − K−1 (F + Q)]−1 K−1 [I − (E + H)K−1 ]−1 N. (4.34) Mivel I − (E + H)K−1 és I − K−1 (F + Q) reguláris alsó ill. fels® háromszögmátrixok, a módszer könnyen végrehajtható egy ún. two-sweep, el®re-hátrafelé számoló algoritmussal 33 http://www.doksihu Megszorozva a 4.30 egyenletet balról I − K−1 (F + Q)-val, majd ezt átrendezve kapjuk, hogy x(m+1) = K−1 {(F + Q)x(m+1) + [I − (E + H)K−1 ]−1 [Nx(m) + k]}. Bevezetve a következ® jelölést c(m+1) ≡ [I − (E + H)K−1 ]−1 [Nx(m) + k], és megszorozva ezt balról I − (E +

H)K−1 -zel, az eredeti módszerb®l egy két lépéses iterációt kapunk: ) c(m+1) = (E + H)K−1 c(m+1) + Nx(m) + k, x(m+1) = K−1 [(F + Q)x(m+1) + c(m+1) ], m = 0, 1, 2. (4.35) Mivel (E + H)K−1 és K−1 (F + Q) szigorú alsó ill. fels® háromszögmátrixok, a c(m+1) vektor komponenseit index szerinti növekv® sorrendben (forward elimination sweep), x(m+1) elemeit pedig index szerinti csökken® sorrendben (backward substitution sweep) számolja az algoritmus. A 4.35 egyenletrendszer a prefaktorizációs módszereknek egy általános alakja, melyb®l H és Q mátrixok megválasztásával deniálhatjuk a konkrét iterációkat. Nézzük, hogy a már korábban megismert speciális módszereket hogyan kapjuk meg a megfelel® mátrixok behelyettesítésével, és hogyan deniáljuk az EWA és AGA iterációkat. 1. A Jacobi-módszer H = −E, Q = −F, KJ = D MJ = D, NJ = E + F −1 B = MJ NJ = D−1 (E + F). 2. A Gauss-Seidel módszer I. típusú H = −E, Q = 0, KG = D MG

= D(I − D−1 F), NG = E −1 −1 G = MG NG = (I − D F)−1 D−1 E. II. típusú H = 0, Q = −F, KG = D −1 f G = D(I − D E), eG = F M N 34 http://www.doksihu e G = (I − D−1 E)−1 D−1 F. f −1 N Ge = M G 3. Az EWA-módszer H = Q = 0, KE = D − diag{EK−1 E F} −1 −1 ME = (I − EKE )KE (I − KE F), NE = odiag{EK−1 E F} −1 −1 −1 −1 E = M−1 E NE = (I − KE F) KE (I − EKE )NE . 4. Az AGA-módszer H = HA , Q = QA , KA = D − diag{(E + HA )K−1 A (F + QA )} −1 MA = [I − (E + HA )K−1 A ]KA [I − KA (F + QA )] NA = odiag{(E + HA )K−1 A (U + QA )} − HA − QA −1 −1 −1 −1 −1 A = M−1 A NA = [I − KA (F + QA )] KA [I − KA (E + HA )] NA . Az AGA módszer az algoritmusok egy széles osztályát reprezentálja. A konkrét iterációt HA és QA megválasztása határozza meg, melyekre kikötés, hogy HA + QA 6= 0 (kivéve a HA = −E és QA = −F eseteket), és feltétel az is, hogy NA nemnulla elemeinek elhelyezkedése nem

egyezik meg az E + HA + F + QA mátrix nemnulla elemeinek elhelyezkedésével. Mivel HA és QA szigorú alsó ill. fels® háromszögmátrixok, nemnulla elemeiket közvetlenül ki tudjuk számolni az (E + HA )K−1 A (F + QA ) implicit alakból. Tehát, ha megadjuk, hogy HA -nak és QA -nak hol helyezkednek el a nemnulla elemei, akkor egy rekurzív, szimultán algoritmussal meg tudjuk határozni KA , HA , QA , NA elemeit úgy, hogy a fenti kikötések teljesüljenek. 4.4 Megjegyzés Ha azt kapjuk, hogy NA minden eleme 0 lesz, akkor a Gauss-eliminációval megegyez® direkt AGA-módszert deniálja a képlet 4.10 Tétel Legyen a Jacobi-mátrix B = D−1 (E + F) egy n × n-es nemnegatív mátrix (n > 2), melynek diagonális elemei 0-k, és ρ(B) < 1. Továbbá legyenek G = M−1 G NG , −1 E = M−1 E NE , A = MA NA a korábban deniált Gauss-Seidel, EWA- és AGA-módszerek iterációs mátrixai. Ekkor ezen mátrixok mindegyike konvergens és 0 ≤ ρ(A) ≤ ρ(E) ≤ ρ(G) ≤

ρ(B) < 1. 35 http://www.doksihu A korábbi táblázatban az utolsó négy reguláris felbontás közül a 12. a Jacobi-módszert, a 11. az I típusú, a 10 pedig a II típusú Gauss-Seidel módszert szemlélteti, és a 9 EWAtípusú mátrix szeletelés Látható, hogy a spektrálsugarakra valóban teljesül a tétel: −1 −1 −1 ρ(M−1 9 N9 ) = ρ(M10 N10 ) < ρ(M11 N11 ) < ρ(M12 N12 ). Az A mátrix 4.29-ben deniált felbontására csak egyetlen AGA-módszer létezik, mivel az algoritmussal azt kapjuk, hogy NA a nullmátrix, MA pedig az A mátrix lesz. A prefaktorizációs iterációs módszerekr®l részletesebb leírás található Woznicki: Nonnegative Splitting Theory [1] cím¶ jegyzetének 4. fejezetében 36 http://www.doksihu 5. fejezet Összefoglalás A mátrix szeletelés elmélet a lineáris algebrai egyenletrendszerek konvergens iterációs megoldási módszereinek összehasonlítását teszi lehet®vé. A második fejezetben az algebrai

egyenletrendszerekkel modellezhet® számos matematikai probléma közül mutattam be kett®t, és azok numerikus megoldását ábrákkal szemléltettem. A harmadik fejezetben az iterációs módszerek konvergenciájának fogalmával és összehasonlításával foglalkoztam, majd néhány klasszikus iterációs módszert mutattam be. A negyedik fejezetben a nemnegatív mátrix szeletelés elméletének alapjait és annak egy alkalmazását (prefaktorizációs módszerek) ismertettem. 37 http://www.doksihu 6. fejezet Függelék Taylor-sorfejtés: Tekintsük az u(x, y) kétváltozós függvényt. Ha u(x, y) legalább (k + 1)-szer folytonosan dierenciálható, akkor egy (x0 , y0 ) pont körül Taylor-sorba fejthet®: u(x, y) = u(x0 , y0 ) + (x − x0 ) ∂u ∂u (x − x0 )2 ∂ 2 u (x0 , y0 ) + (y − y0 ) (x0 , y0 ) + (x0 , y0 )+ ∂x ∂y 2! ∂x2 (y − y0 )2 ∂ 2 u (x − x0 )(y − y0 ) ∂ 2 u (x , y ) + (x0 , y0 ) + . 0 0 2! ∂y 2 2! ∂x∂y X (x − x0 )k1 (y − y0

)k2 ∂ k+1 u X (x − x0 )k1 (y − y0 )k2 ∂ k u (x , y )+ (θ), + 0 0 k1 !k2 ! ∂ k1 x∂ k2 y k1 !k2 ! ∂ k1 x∂ k2 y k ,k ≥0 k ,k ≥0 + 1 2 k1 +k2 =k+1 1 2 k1 +k2 =k ahol X k1 ,k2 ≥0 k1 +k2 =k+1 (x − x0 )k1 (y − y0 )k2 ∂ k+1 u (θ) k1 !k2 ! ∂ k1 x∂ k2 y a Lagrange-féle maradéktag. A Dirichlet-problémában a rácsháló egy bels® pontja, (x0 , y0 ) körül Taylor-sorba fejtjük a függvényt annak szomszédos pontjaiban: u(x0 ± h, y0 ) = u(x0 , y0 ) ± h ∂u h2 ∂ 2 u h3 ∂ 3 u h4 ∂ 4 u (x0 , y0 ) + (x , y ) ± (x , y ) + (θ1 ) 0 0 0 0 ∂x 2 ∂x2 3! ∂x3 4! ∂x4 u(x0 , y0 ± h) = u(x0 , y0 ) ± h ∂u h2 ∂ 2 u h3 ∂ 3 u h4 ∂ 4 u (x0 , y0 ) + (x , y ) ± (x , y ) + (θ2 ) 0 0 0 0 ∂y 2 ∂y 2 3! ∂y 3 4! ∂y 4 38 http://www.doksihu Ezt a 4 egyenletet összeadva a függvényt másodrendben közelítjük (x0 , y0 ) pontban: 1 [u(x0 + h, y0 ) + u(x0 − h, y0 ) + u(x0 , y0 + h) + u(x0 , y0 − h) − 4u(x0 , y0 )] = h2 = h ∂ 2u

∂x (x0 , y0 ) + 2 i h2 h ∂ 4 u i ∂ 2u ∂ 4u (x , y ) + (θ ) + (θ ) 0 0 1 2 ∂y 2 12 ∂x4 ∂y 4 A másod- és annál magasabb fokú tagokat elhagyva, továbbá felhasználva, hogy ∂ 2u ∂ 2u (x , y ) + (x0 , y0 ) = −g(x0 , y0 ), 0 0 ∂x2 ∂y 2 kapjuk, hogy h2 1 u(x0 , y0 ) ≈ [u(x0 + h, y0 ) + u(x0 − h, y0 ) + u(x0 , y0 + h) + u(x0 , y0 − h)] + g(x0 , y0 ). 4 4 6.1 Deníció Legyen (X, ρ) egy metrikus tér ρ metrikával Egy f:XX függvény kontrakció, ha létezik olyan q ∈ [0, 1), melyre ρ(f (x1 ), f (x2 )) ≤ q·ρ(x1 , x2 ) minden x1 , x2 ∈ X re 6.1 Tétel (Banach-féle xponttétel) Ha (X, ρ) egy teljes metrikus tér és f:XX egy kontrakció, akkor f-nek létezik egyértelm¶ xpontja, és az x(m+1) = f (x(m) ) iteráció (m=0,1,2.) tetsz®leges x(0) esetén ehhez a xponthoz konvergál. Mátrixok típusai 6.2 Deníció Egy n × n-es komplex elem¶ A mátrix adjungáltjának nevezzük azt az A∗ mátrixot, melyre a∗i,j = āj,i . 6.3

Deníció Egy n × n-es komplex elem¶ A mátrixot hermitikus mátrixnak nevezünk, ha A = A∗ , azaz ha ai,j = āj,i . 6.4 Deníció Egy n × n-es komplex elem¶ A mátrixot normálisnak nevezünk, ha kommutál az adjungáltjával, azaz A∗ A = AA∗ . 39 http://www.doksihu 6.1 Következmény Minden hermitikus mátrix normális mátrix 6.5 Deníció Minden n ≥ 2-re, egy n × n-es komplex A mátrix reducibilis, ha létezik olyan n × n-es permutáló P mátrix, hogy " A1,1 A1,2 PAPT = 0 # , A2,2 ahol A1,1 egy r × r-es, A2,2 pedig egy (n − r) × (n − r)-es részmátrix (1 ≤ r < n). Ha nem létezik ilyen P permutáló mátrix, akkor A irreducibilis. 6.6 Deníció Egy A = (ai,j ) n × n-es komplex mátrix diagonálisan domináns, ha |ai,i | ≥ n X ∀ 1 ≤ i ≤ n-re. |ai,j | j=1 j6=i Szigorúan diagonálisan domináns A, ha |ai,i | > n P ∀ 1 ≤ i ≤ n-re. Irreducibilisen |ai,j | j=1 j6=i diagonálisan domináns A, ha irreducibilis

és |ai,i | > n P |ai,j | legalább egy i-re. j=1 j6=i 6.2 Tétel Legyen A = (ai,j ) egy n × n-es szigorúan vagy irreducibilisen diagonálisan domináns komplex mátrix. Ekkor A reguláris 6.7 Deníció Egy A n×n-es komplex mátrix konvergens (0-hoz tart), ha az A, A2 , A3 sorozat a 0 mátrixhoz tart. Vektor- és mátrixnormák 6.8 Deníció Legyen x egy Cn -beli vektor Ekkor k · k : Cn R norma Cn -en, ha teljesülnek az alábbiak: 1. kxk ≥ 0, és kxk = 0 ⇐⇒ x = 0 2. kλxk = |λ| · kxk minden x ∈ Cn és λ ∈ C esetén 3. kx + yk ≤ kxk + kyk minden x, y ∈ Cn -re 6.9 Deníció Legyen x egy komplex n-elem¶ oszlopvektor Ekkor x Euklideszi-normája: ∗ kxk2 ≡ (x x) 1/2 = n X j=1 40 2 |xj | 1/2 . http://www.doksihu 6.10 Deníció Legyen x egy komplex n-elem¶ oszlopvektor Ekkor x ∞-normája: kxk∞ ≡ max |xi |. 1≤i≤n 6.11 Deníció Legyen A egy n×n-es komplex mátrix, ekkor A-nak az kxk vektornorma által indukált

mátrixnormája: kAk = sup x6=0 kAxk . kxk 6.3 Tétel Ha A egy n × n-es mátrix, x egy tetsz®leges n-elem¶ vektor, akkor a fent deniált normákra teljesül az alábbi egyenl®tlenség: kAxk ≤ kAk kxk. Spektrálsugárra vonatkozó tételek 6.4 Tétel Legyen A és B két n × n-es mátrix Ekkor ρ(AB) = ρ(BA) 6.5 Tétel Legyen A és B két n × n-es mátrix, melyekre 0 ≤ |B| ≤ A Ekkor ρ(B) ≤ ρ(A). 6.6 Tétel Ha egy G mátrixra ρ(G) < 1, akkor I − G reguláris és (I − G)−1 = I + G + G2 + ., ahol az egyenlet jobb oldalán a sorozat konvergens. Fordítva is igaz: ha a jobb oldali sorozat konvergens, akkor ρ(G) < 1. 41 http://www.doksihu 7. fejezet Köszönetnyilvánítás Köszönöm témavezet®mnek, Mincsovics Miklósnak, hogy segített megismerkedeni ezzel a témával, rámutatott az összefüggésekre, valamint a megfelel® szakirodalom ajánlásával, tanácsaival és ötleteivel hozzájárult a dolgozat elkészüléséhez. Továbbá

köszönöm Chmelik Gábornak a Latex használatában nyújtott sok segítséget. 42 http://www.doksihu Nyilatkozat A szakdolgozat szerz®jeként fegyelmi felel®sségem tudatában kijelentem, hogy a dolgozatom önálló munkám eredménye, saját szellemi termékem, abban a hivatkozások és idézések standard szabályait következetesen alkalmaztam, mások által írt részeket a megfelel® idézés nélkül nem használtam fel. 43 http://www.doksihu Irodalomjegyzék [1] : Nonnegative Splitting Theory, Japan J. Indust, Appl.Math,11, 289-342, 1994 Zbigniew I. Woznicki [2] Richard S. Varga : Matrix Iterative Analysis, Prentice-Hall, New Jersey, 1962 [3] Stoyan Gisbert, Takó Galina [4] Douglas [5] Gilbert Strang : Numerikus módszerek I., Typotex, Budapest, 2002 (http://www.tankonyvtarhu/konyvek/numerikus-modszerek-1/numerikusmodszerek-1-081029-9) : A Concise Introduction to Numerical Analysis, 2001 (http://www.imaumnedu/ arnold/59700-01/nabookpdf) N. Arnold :

Linear Algebra and its Applications, Thomson Learning, 3rd edition, 1988 44