Eötvös Loránd Tudományegyetem Természettudományi kar
Lineáris algebra és mátrixok alkalmazása a numerikus analízisben Szakdolgozat
Témavezető: Fialowski Alice egyetemi docens Algebra és Számelmélet Tanszék
Készítette: Borostyán Dóra Matematika BSc matematikai elemző
Budapest 2013
Tartalomjegyzék Bevezetés ......................................................................................................................................... 4 1.
2.
3.
Előismeretek ....................................................................................................................... 9 1.1.
Mátrixok................................................................................................................ 9
1.2.
Lineáris egyenletrendszerek megoldhatósága .................................................. 14
Direkt módszerek............................................................................................................. 15 2.1.
Gauss-elimináció ................................................................................................ 15
2.2.
LU-felbontás ....................................................................................................... 17
2.3.
Cholesky-felbontás ............................................................................................. 20
Legkisebb négyzetek módszere ...................................................................................... 22 3.1.
4.
Főbb eredmények ............................................................................................... 23
Iterációs módszerek......................................................................................................... 27 4.1.
Jacobi-iteráció ..................................................................................................... 28
4.2.
Gauss-Seidel-iteráció.......................................................................................... 30
4.3.
Richardson-iteráció............................................................................................. 33
Összefoglalás ................................................................................................................................ 36 Köszönetnyilvánítás .................................................................................................................... 37 Irodalomjegyzék .......................................................................................................................... 38 Nyilatkozat.................................................................................................................................... 39
3
Bevezetés Szakdolgozatom témája a lineáris algebra alkalmazása a numerikus analízisben. Bármennyire is idegennek hangzik, ez a témakör jelen van a hétköznapi életünkben több területen is, mint gondolnánk. Fizika, mechanika, kémia, építészet, pénzügy és közgazdaságtan területein előforduló problémákra is alkalmazzuk a lineáris algebrai egyenletrendszerek felírását. Ez a problémamegoldás magában foglalja a mátrixok alkalmazását is, hiszen az egyenletrendszereket mátrixok formájában is felírhatjuk. A lineáris egyenletrendszerek megoldásán kívül a másik fő mátrixszámítási művelet a mátrix sajátértékeinek és sajátvektorainak meghatározása. Motivációnak szeretnék ezekre a módszerekre két mindennapjainkban is előforduló modellre példát mutatni. (A példák és ábráik a 3 -ból származnak.) Elsőként egy fizikai példa: „Egy mechanikai rendszer rezgései” Egy mátrix sajátértékeinek és sajátvektorainak kiszámítása alapvető matematikai segédeszköz egy mechanikai szerkezet rezgéseinek tanulmányozásához. Ebben az esetben a sajátértékek a frekvencia négyzetes értékei, a sajátvektorok pedig a rendszer rezgésmódjai. Vegyük például egy épület rezgési frekvenciájának kiszámítását, amely fontos szerepet játszik abban, hogy meghatározzuk az épület szilárdságát egy földrengés esetén. Az egyszerűség kedvéért ezt most egy egyszerűbb modellen vizsgáljuk meg, de a módszer ugyanez bonyolultabb modellekre is. Tekintsünk egy kétemeletes épületet (földszint + két emelet), amelynek a merev mennyezeteinek tömegeit 𝑚 1, 𝑚 2, 𝑚 3-mal jelöljük. A falaknak elhanyagolható a tömege, de a rugalmasságuk egy rugóval van modellezve, amelynek a rugóállandói 𝑘1, 𝑘2 , 𝑘3 (minél nagyobb 𝑘𝑖 , annál szilárdabb a fal). A mennyezet vízszintes elmozdulását 𝑦1, 𝑦2, 𝑦3-mal jelöljük, az épület alapja pedig földhöz rögzített. Másképpen, ez a kétemeletes épület egy olyan rendszert határoz meg, melyben a 3 tömeg rugókkal kapcsolódik egy fix talpazathoz.
4
Ez az alábbi ábrán is látható.
A mechanika törvényeiből tudjuk, hogy a tömeg és a gyorsulás szorzata egyenlő az alkalmazott erők összegével. Az itt használt erők a rugók általi visszaható erők (kötélerők). Ezek egyenlők a rugó merevségének és nyúlásának szorzatával. Az 𝑦1, 𝑦2, 𝑦3 irányú elmozdulások a t idő függvényei. Ezek első deriváltjai, 𝑦1, 𝑦2, 𝑦3, jelölik a sebességeket, és a második deriváltak, 𝑦1, 𝑦2, 𝑦3, pedig az 𝑚 1, 𝑚 2, 𝑚 3 súlyok gyorsulásait. Így az alábbi 3 egyenletet kapjuk: 𝑚1 𝑦1 + 𝑘1 𝑦1 + 𝑘2 𝑦1 − 𝑦2 = 0, 𝑚 2𝑦2 + 𝑘2 𝑦2 − 𝑦1 + 𝑘3 𝑦2 − 𝑦3 = 0, 𝑚 3 𝑦3 + 𝑘3 𝑦3 − 𝑦2 = 0,
(1)
amit mátrix formájában értelmezünk: 𝑀𝑦 + 𝐾𝑦 = 0,
(2)
ahol 𝑦1 𝑦 = 𝑦2 , 𝑦3
𝑚1 0 0 𝑀 = 0 𝑚2 0 , 0 0 𝑚3
𝑘1 + 𝑘2 −𝑘2 𝐾= 0
−𝑘2 𝑘2 + 𝑘3 −𝑘3
0 −𝑘3 . 𝑘3
Az 𝑀 mátrix a tömegmátrix, a 𝐾 mátrix pedig a merevség mátrix. Mindkét mátrix szimmetrikus. Az (1) egyenletrendszer azon pontos megoldásait keressük, amelyek időben periodikusak, hogy ábrázolni tudjuk majd a rendszer rezgéseit. Ezért a következő helyettesítést végezzük: 𝑦 𝑡 = 𝑦 0 𝑒 𝑖𝜔𝑡 ,
5
ahol i a képzetes egység, és ω a megoldás rezgési frekvenciája. Egyszerű számítás után megkapjuk, hogy ebben az esetben a gyorsulás 𝑦 𝑡 = −𝜔2 𝑦(𝑡), és ez leegyszerűsíthető a (2) felírás szerint a következőre: 𝐾𝑦 0 = 𝜔2 𝑀𝑦 0 .
(3)
Ha az összes súly 1-gyel egyenlő, akkor 𝑀 az egységmátrix, és (3) egy szokásos sajátérték probléma a 𝐾 mátrixra, azaz 𝑦 0a 𝐾 𝜔2 sajátértékhez tartozó sajátvektora. Ha a súlyok bármilyen értéket felvehetnek, akkor (3) egy általánosított sajátérték probléma, azaz 𝜔2 az 𝑀 −1/2𝐾𝑀 −1/2 mátrix sajátértéke. Természetesen, ha egy 𝑛-emeletes épületet tekintünk, akkor egy ugyanilyen mátrix problémát kapnánk, csak 𝑛 + 1-es méretben. Az 𝑀 mátrix mindig diagonális, a 𝐾 pedig mindig tridiagonális lesz.
Második példám egy képtömörítési módszert mutat be: „Képtömörítés SVD faktorizációval” Egy fekete-fehér képet tekinthetünk olyan derékszögű 𝐴 mátrixnak, melynek mérete azonos a pixelek számával. A mátrixba írt 𝑎 𝑖𝑗 adatokat a 0,1 intervallumból vesszük, ahol a 0 egy fehér pixelt, az 1 pedig egy fekete pixelt jelképez. A közbeeső értékek, 0 < 𝑎 𝑖𝑗 < 1, pedig a szürke különböző árnyalatainak felelnek meg. Tegyük fel, hogy a kép mérete olyan nagy, hogy nem tudjuk számítógépen tárolni vagy emailben elküldeni. Most megmutatjuk, hogy az SVD (singular value decomposition) hogyan minimalizálja a kép méretét úgy, hogy egy, a képhez látszólag közel hasonló képet felcserél az eredetivel. Egy 𝐴 ∈ ℂ𝑚×𝑛 𝑟 rangú mátrix SVD faktorizációja a következő: Σ 0 , 0 0 ahol 𝑈 ∈ ℂ𝑛 és 𝑉 ∈ ℂ𝑚 két unitér mátrix és Σ egy diagonális mátrix 𝜇 1,… , 𝜇 𝑟 elemekkel, 𝐴 = 𝑉Σ𝑈∗ é𝑠 Σ =
amelyekre 𝜇1 ≥ 𝜇 2 ≥ ⋯ ≥ 𝜇 𝑟 > 0. Ezek az értékek az 𝐴∗𝐴 (ahol 𝐴∗ az 𝐴 mátrix adjungáltja) sajátértékeinek pozitív négyzetgyökei, amelyeket 𝐴 „szinguláris értékeinek” (singular values) nevezünk. Tehát az SVD faktorizációt sajátérték-problémának is tekinthetjük. Ha 𝑈 és 𝑉 oszlopait 𝑢𝑖 és 𝑣𝑖 -ként jelöljük, akkor az 𝐴 mátrix SVD faktorizációja a következőképp is felírható: 𝑟
𝐴 = 𝑉Σ 𝑈∗ =
𝜇 𝑖 𝑣𝑖 𝑢∗𝑖 . 𝑖=1
6
Mivel a szinguláris értékek csökkenő sorrendben vannak, így 𝐴 könnyen közelíthető az első 𝑘 ≤ 𝑟 taggal. 𝑘
𝜇 𝑖 𝑣𝑖 𝑢∗𝑖 .
𝐴𝑘 = 𝑖=1
𝐴𝑘 lesz az 𝐴 legjobb közelítése a 𝑘 rangú mátrixok között. Természetesen 𝑘 sokkal kisebb, mint 𝑟, így 𝐴 közelítése 𝐴𝑘 -val sok memóriát spórol meg. Az 𝐴 ∈ ℝ𝑚×𝑛 mátrix tárolása 𝑚 × 𝑛 skalárt igényel. Viszont 𝐴𝑘 tárolásához elég 𝑘 darab 𝜇 𝑖 𝑣𝑖 ∈ ℂ𝑚 vektor és 𝑘 darab 𝑢𝑖 ∈ ℂ𝑛 vektor, azaz, 𝑘 𝑚 + 𝑛 skalár. Így tehát jobban megéri, ha 𝑘 kicsi, és megelégszünk 𝐴 egy közelítésével is. Az ábrán az eredeti kép 500×752 pixelből áll, így az ennek megfeleltetett 𝐴 mátrix mérete 500×752. Az eredeti kép mellett még további 3 kép látható, amelyek 3 különböző, 𝐴 -t közelítő 𝐴𝑘 mátrixnak felelnek meg. Az első közelítés 𝑘 = 10 értékkel nagyon elmosódott
képet ad, de 𝑘 = 20-ra a lényeg már látható. Viszont a 𝑘 = 60-ra kapott kép és az eredeti között nincs szinte semmi különbség, mégis a közelített értékkel a tárolási helyet az ötödére csökkentettük: 𝑘(𝑚 + 𝑛) 60(500 + 752) = ≈ 20%. 𝑚𝑛 500 × 752 Ennek a képtömörítési módszernek a fő számítási költsége az 𝐴 mátrix SVD faktorizációja. Mint korábban kifejtettük, az 𝐴 szinguláris értékei az 𝐴∗ 𝐴 mátrix sajátértékeinek a pozitív négyzetgyökei. Így az SVD faktorizáció egy lehetséges módszer egy mátrix sajátértékeinek és sajátvektorainak meghatározásához. Megjegyzem még, hogy természetesen képtömörítéshez léteznek más, sokkal hatékonyabb és olcsóbb algoritmusok is, mint az SVD faktorizáció.
7
A példák után ismertetem a szakdolgozatom tartalmát. Az első fejezetben felsorolom a definíciókat, tételeket, melyek ehhez a témakörhöz előismeretként szükségesek. A többi fejezetben pedig a különböző alkalmazásokat részletezem alaposabban. A fő fejezetek a direkt és iterációs módszerek és a legkisebb négyzetek módszere. Igyekszem példákkal is alátámasztani a módszerek hasznosságát az érthetőség érdekében.
8
1. Előismeretek 1.1. Mátrixok Tekintsünk egy 𝐴 ∈ 𝐹 𝑘 ×𝑛-es mátrixot az 𝐹 test fölött. Felsorolok néhány fontosabb mátrixtulajdonságot, illetve alapvető elnevezéseket, fogalmakat, melyek szükségesek lesznek a fejezetek megértéséhez. (Mindezek az 1 , 2 , 4 , és 6 forrásokból származnak.)
A mátrixok körében az összeadás asszociatív és kommutatív, a szorzás szintén asszociatív, viszont nem kommutatív. Továbbá a szorzás művelete disztributív az összeadásra nézve. Emiatt a következők igazak 𝜆 ∈ 𝐹, és 𝐴, 𝐵, 𝐶 mátrixokra: o
𝐴+𝐵 +𝐶 = 𝐴+ 𝐵+𝐶 ;
o 𝐴 + 𝐵 = 𝐵 + 𝐴; o 𝐴 𝐵𝐶 = 𝐴𝐵 𝐶; o 𝐴 𝐵 + 𝐶 = 𝐴𝐵 + 𝐴𝐶, 𝐴 + 𝐵 𝐶 = 𝐴𝐵 + 𝐴𝐶;
Ha két mátrix azonos méretű, akkor a köztük lévő relációkat ("=", "<", ">","≤", "≥") elemenként végezzük.
Azokat a mátrixokat, amelyeknek minden eleme nulla, nullmátrixnak nevezzük. Jelölése: 𝟎.
Az 𝐴 ∈ 𝐹 𝑘 ×𝑛 mátrix transzponáltján azt a 𝐵 ∈ 𝐹 𝑛×𝑘 mátrixot értjük, amelyre 𝑏𝑖𝑗 = 𝑎𝑗𝑖 . Jelölése: 𝐴𝑇 .
Az 𝐴 valós mátrixot szimmetrikusnak nevezzük, ha 𝐴𝑇 = 𝐴.
Azokat a mátrixokat, amelyeknek ugyanannyi sora van, ahány oszlopa, négyzetes mátrixoknak nevezzük.
𝐴 négyzetes mátrix esetén 𝐴 inverzének nevezzük azt a mátrixot, amellyel 𝐴-t balról vagy jobbról megszorozva, az egységmátrixot kapjuk eredményül. Jelölése: 𝐴−1 (𝐴𝐴−1 = 𝐸). Továbbá, ha egy mátrixnak van inverze, akkor a mátrixot reguláris mátrixnak nevezzük.
Egy 𝐴 valós négyzetes mátrixot ortogonálisnak nevezünk, ha van inverze és 𝐴−1 = 𝐴𝑇 . Jele: ⊥
9
Egy 𝐴 négyzetes mátrix diagonális, ha minden főátlón kívüli eleme nulla, azaz 𝑎 𝑖𝑗 = 0 ∀ 𝑖 ≠ 𝑗-re. Jelölése: 𝑑𝑖𝑎𝑔(𝑎, 𝑏, … , 𝑧), ahol 𝑎, 𝑏, … , 𝑧 a főátló elemei.
Egységmátrixnak hívjuk azt a speciális négyzetes diagonális mátrixot, melynek főátlójában csupa 1-esek szerepelnek (a többi eleme 0). Jelölése: 𝐸. Erre igaz: 𝐸𝐴 = 𝐴𝐸 = 𝐴.
Egy 𝐴 mátrixot felső háromszögmátrixnak nevezünk, ha a főátló alatti elemei mind nullák. Pl.: ∗ 0 𝐴= 0 0
∗ ∗ 0 0
∗ ∗ ∗ 0
∗ ∗ ∗ . ∗
A * jel helyén bármilyen 𝐹-beli szám állhat.
Egy 𝐴 mátrixot alsó háromszögmátrixnak nevezünk, ha a főátló feletti elemei mind nullák. Pl.: ∗ 𝐴 = ∗∗ ∗
0 0 0 ∗ 0 0 ∗ ∗ 0 ∗ ∗ ∗
Egy 𝐴 négyzetes mátrix determinánsát det(A)-val jelöljük. Egy négyzetes mátrixnak pontosan akkor van inverze, ha determinánsa nem nulla: det A ≠ 0. Fontos még a determinánsok szorzási szabálya: két 𝐴, 𝐵 négyzetes mátrixra: det(AB) = det(A) det(B).
Egy 𝑛 × 𝑘-s mátrix rangja a mátrix lineárisan független oszlopainak vagy sorainak maximális száma. Másképpen, a rang a mátrix oszlopvektorai vagy sorvektorai által kifeszített altér dimenziója az 𝐹 𝑘 vektortérben. Jelölése: 𝑟(𝐴).
Következnek a szükséges definíciók, tételek.
Definíció. Legyenek 𝑉1 és 𝑉2 ugyanazon 𝐹 kommutatív test feletti vektorterek. A 𝑉1 -ről 𝑉2 -be ható 𝒜 függvényt (homogén) lineáris leképezésnek nevezzük, ha művelettartó, azaz 1. minden 𝑢, 𝑣 ∈ 𝑉1 -re 𝒜 𝑢 + 𝑣 = 𝒜𝑢 + 𝒜𝑣 ; 2. minden 𝑢 ∈ 𝑉1 , 𝜆 ∈ 𝐹-re 𝒜 𝜆𝑢 = 𝜆(𝒜𝑢).
10
Definíció. Legyen 𝒜 lineáris leképezés 𝑉1 -ről 𝑉2 -be. Az 𝒜 leképezés képtere a képelemek halmaza, ezt 𝐼𝑚𝒜-val jelöljük. Tehát 𝐼𝑚𝒜 = 𝑦 ∈ 𝑉2 ∃𝑥 ∈ 𝑉1 𝒜𝑥 = 𝑦} = 𝒜𝑥 𝑥 ∈ 𝑉1 .
Definíció. Legyen 𝒜 lineáris leképezés 𝑉1 -ről 𝑉2 -be. Az 𝒜 leképezés magtere a 𝑉2 nullvektorára képződő elemek halmaza, ezt 𝐾𝑒𝑟𝒜-val jelöljük. Tehát 𝐾𝑒𝑟𝒜 = 𝑥 ∈ 𝑉1 𝒜𝑥 = 0 .
Megjegyzés. 𝐼𝑚𝒜 és 𝐾𝑒𝑟𝒜 alteret alkotnak. Legyen mostantól az alaptest ℂ (komplex számok teste).
Definíció. Legyen 𝐴 ∈ ℂ𝑛×𝑛 egy tetszőleges négyzetes mátrix. Ha egy nullvektortól különböző 0 ≠ 𝑣 ∈ ℂ𝑛 vektor és egy 𝜆 ∈ ℂ szám esetén teljesül az 𝐴𝑣 = 𝜆𝑣 egyenlőség, akkor a 𝑣 vektort a mátrix sajátvektorának, és a 𝜆 számot a sajátvektorhoz tartozó sajátértéknek nevezzük (és fordítva). Egy összetartozó sajátértéket és sajátvektort sajátpárnak nevezünk.
Tétel. 1. Minden sajátvektorhoz csak egy sajátérték tartozik. 2. Egy adott 𝜆 sajátértékhez tartozó összes sajátvektor és a 0 alteret alkotnak. Ezt az alteret a 𝜆-hoz tartozó sajátaltérnek nevezzük. Egy 𝑛 × 𝑛-es négyzetes mátrixnak pontosan 𝑛 darab sajátértéke van multiplicitással. Ezeket a det 𝐴 − 𝜆𝐸 = 0 egyenletből kapjuk. Az ezekhez a sajátértékekhez tartozó sajátvektorokat pedig az 𝐴 − 𝜆𝐸 𝑣 = 0 egyenletrendszer nem nulla (𝑣 ≠0) megoldásai adják.
11
Tétel. Egy adott 𝐴 ∈ ℂ𝑛×𝑛 mátrix esetén a 𝑝𝐴 𝜆 = det(𝐴 − 𝜆𝐸) polinomot a mátrix karakterisztikus polinomjának, a 𝑃𝐴 𝜆 = 0 egyenletet pedig karakterisztikus egyenletnek nevezzük.
Definíció. A (𝑉, ∙ ) párt normált térnek hívjuk, ha 𝑉 egy vektortér, és ∙ ∶ 𝑉 → ℂ egy adott függvény, ún. norma, az alábbi tulajdonságokkal: 1.
𝑥 = 0 ⇔ 𝑥 = 0,
2.
𝛼𝑥 = 𝛼 ∙ 𝑥 , ∀𝑥 ∈ 𝑉, ∀𝛼 ∈ 𝕂,
3.
𝑥 + 𝑦 ≤ 𝑥 + 𝑦 , ∀𝑥, 𝑦 ∈ 𝑉 (háromszög-egyenlőtlenség).
Például a ℂ𝑛 vektortér normált tér, ha az 𝑥 = 𝑥1 ,𝑥2 , … , 𝑥𝑛
𝑇
vektor normáját az alábbi
módon definiáljuk: 𝑥
𝑝
=
𝑝
𝑥1
𝑝
+ 𝑥2 𝑝 +.. . + 𝑥𝑛
𝑝
𝑝 = 1,2, … .
Az ilyen típusú normák közül a leggyakoribbak a következők:
1-es norma (oktaédernorma): 𝑥
2-es norma (euklideszi norma): 𝑥
maximumnorma: 𝑥
∞
1
= 𝑥1 + 𝑥2 +.. . + 𝑥𝑛 , 2
=
2
𝑥1 2 + 𝑥2 2 +.. . + 𝑥𝑛 2,
= max 𝑥1 , 𝑥2 ,… , 𝑥𝑛 .
Ezeket a normákat vektornormáknak nevezzük.
Tétel. Tegyük fel, hogy a ℂ𝑛 és a ℂ𝑚 normált terekben ugyanazt a vektornormát használjuk. Ekkor a korábban megismert vektornormák az alábbi mátrixnormákat indukálják:
Oktaédernorma 𝑝 = 1 : 𝐴
Maximumnorma 𝑝 = ∞ : 𝐴
∞
Euklideszi norma 𝑝 = 2 : 𝐴
2
1
= 𝑚𝑎𝑥𝑗 =1,…,𝑛
𝑚 𝑖 =1
= 𝑚𝑎𝑥𝑗 =1,…,𝑚 =
𝑎 𝑖𝑗 (oszlopösszegnorma),
𝑛 𝑗 =1
𝑎 𝑖𝑗 (sorösszegnorma),
𝜌(𝐴∗𝐴), ahol 𝜌 az 𝐴 mátrix spektrálsugara, és
𝐴∗ az 𝐴 mátrix transzponált konjugáltja.
12
Definíció. Legyen 𝐴 ∈ ℂ𝑛×𝑛 , 𝜆𝑖 az 𝐴 sajátértékei, 𝑖 = 1, … , 𝑛. Spektrálsugárnak hívjuk az abszolút értékben legnagyobb sajátértéket, 𝜌 𝐴 = max 𝜆𝑖 . 1≤𝑖≤𝑛
Definíció. 𝐴 ∈ ℂ𝑚 ×𝑛 mátrix adjungáltja 𝐴∗ = 𝐴𝑇 . 𝐴 önadjungált, ha 𝐴∗ = 𝐴.
Definíció. Egy 𝑚 × 𝑛-es mátrixnak pszeudoinverze egy olyan 𝐴† 𝑛 × 𝑚-es mátrix, amelyre 1. 𝐴𝐴† 𝐴 = 𝐴, 2. 𝐴† 𝐴𝐴† = 𝐴† , 3. 𝐴𝐴† önadjungált, 4. 𝐴† 𝐴 önadjungált.
Definíció. 𝐴 ∈ ℂ𝑛×𝑛 szimmetrikus, pozitív definit (SZPD) mátrix, ha 𝐴 = 𝐴𝑇 és 𝐴𝑥, 𝑥 > 0 ∀𝑥 ∈ ℂ𝑛 ,𝑥 ≠ 0.
Definíció. 𝐴 ∈ ℂ𝑛×𝑛 szigorúan diagonálisan domináns (SZDD), ha 𝑎 𝑖𝑖 >
𝑛 𝑗 =1 𝑗 ≠𝑖
𝑎 𝑖𝑗
𝑖 = 1, … , 𝑛.
Definíció. 𝐴 ∈ ℂ𝑛×𝑛 M-mátrix, ha
𝑎 𝑖𝑗 ≤ 0 ∀ 𝑖 ≠ 𝑗 és
∃𝑔 ∈ ℝ𝑛 ,𝑔 > 0: 𝐴𝑔 > 0.
Tétel. Legyen az 𝐴 ∈ ℝ𝑛×𝑛 mátrix olyan, hogy a főátlóján kívüli elemek nempozitívak. Ekkor 𝐴 pontosan akkor M-mátrix, ha van olyan 𝑔 > 0 vektor, mellyel 𝐴𝑔 > 0.
13
1.2. Lineáris egyenletrendszerek megoldhatósága Egy 𝑛 ismeretlenes, 𝑘 egyenletből álló lineáris egyenletrendszer általános alakját a következőképpen definiáljuk: 𝑎 11𝑥1 𝑎 21 𝑥1
+ +
⋯ + 𝑎 1𝑛 𝑥𝑛 ⋯ + 𝑎 2𝑛 𝑥𝑛
𝑎 𝑘1 𝑥𝑘
+
…
+
𝑎 𝑘𝑛 𝑥𝑛
= 𝑏1 = 𝑏2 ⋮ = 𝑏𝑘 ,
ahol 𝑎 𝑖𝑗 és 𝑏𝑗 adottak (𝑖 = 1, … , 𝑘 ; 𝑗 = 1, … , 𝑛), és keressük azokat az 𝑥𝑗 (𝑗 = 1, … , 𝑛) számokat, amelyek teljesítik a fenti egyenleteket. Az egyenletrendszer mátrixos alakban is felírható, ehhez bevezetjük az 𝑎 11 𝑎 12 ⋯ 𝑎 1𝑛 𝑎 21 𝑎 22 ⋯ 𝑎 2𝑛 𝐴= ⋮ 𝑎 𝑘1 𝑎 𝑘2 ⋯ 𝑎 𝑘𝑛 együtthatómátrixot. Így az egyenletrendszer a következőképp értelmezhető: 𝐴𝑥 = 𝑏, ahol 𝑏 az egyenletrendszer jobb oldalának értékeit tartalmazza egy 𝑘 elemű oszlopvektorban, 𝑥 pedig egy 𝑛 elemű oszlopvektor, amely a megoldás vektorát jelöli. Célunk: megkeresni azt az 𝑥-t, amellyel 𝐴𝑥 = 𝑏. Ennek az új alakban felírt egyenletrendszernek a megoldhatóságát írja le a következő tétel.
Tétel. 2 Egy 𝐴𝑥 = 𝑏 egyenletrendszer akkor és csak akkor megoldható, ha az 𝐴 együtthatómátrix és a 𝑏 vektorral kibővített együtthatómátrix rangja megegyezik: 𝑟 𝐴 = 𝑟(𝐴|𝑏). Ha az egyenletrendszer megoldható és 𝑟 𝐴 < 𝑛, akkor végtelen sok megoldás van, ha 𝑟 𝐴 = 𝑛, akkor egyértelmű a megoldás. A következő fejezetekben feltesszük, hogy az egyenletrendszerünk mátrixa négyzetes és valós 𝐴 ∈ ℝ𝑛×𝑛 , és csak egyértelmű megoldásunk van. Ennek szükséges és elégséges feltétele, hogy a mátrix determinánsa nullától különbözzön, det 𝐴 ≠ 0.
14
2. Direkt módszerek A lineáris egyenletrendszerek megoldásának és a sajátérték feladatoknak a két fő numerikus megoldási módja a direkt vagy az iterációs módszer. Elsőként a direkt módszerekkel foglalkozunk, mely alkalmazásánál a megoldást véges számú lépés után pontosan megkapjuk. Ezen belül pedig elsőként az egyik leggyakoribb eljárással, a Gauss-eliminációval ismerkedünk meg.
2.1. Gauss-elimináció Ezen eljárás során célunk, hogy az 𝐴𝑥 = 𝑏 egyenletrendszer együtthatómátrixát aritmetikai műveletekkel felső háromszögmátrixszá alakítsuk. Ez azért fontos, mert az ilyen alakú mátrixokból és a 𝑏 vektor segítségével könnyen megkapjuk az egyenletekbe visszafelé történő helyettesítésekkel a megoldást. Tekintsünk egy példát: 𝑥1 1 −1 2 5 𝐴 = 2 −2 1 , 𝑥 = 𝑥2 , 𝑏= 1 𝑥3 3 −2 7 20 𝑥 1 −1 2 5 1 2 −2 1 𝑥2 = 1 3 −2 7 𝑥3 20 A módszer során eltávolítjuk 𝑥1-et a második és harmadik egyenletből, majd 𝑥2-t a harmadik egyenletből. Ennek eredményeképpen, 𝑥3-t megkapjuk a harmadik egyenletből, ezután 𝑥2-t a másodikból, majd 𝑥1-t az elsőből. A következő műveleteket végezhetjük a mátrix sorai között:
egy sor nem nulla skalárral való szorzása;
egy sor konstansszorosának hozzáadása egy másik sorhoz;
két sor felcserélése.
Mindezt illusztrálva, elvégezzük az eljárást a példán. Az elimináció során a 𝑏-vel kiegészített együtthatómátrixból indulunk ki: 1 −1 2 5 2 −2 1 1 3 −2 7 20
15
Első lépésben mindig meghatározzuk a 𝑝 elemet, amely ekkor a mátrix főátlójának első eleme, azaz az 1,1 indexelésű elem lesz, jelen esetben 𝑝 = 1. Ez azért szükséges, mert e szerint az elem szerint fogjuk elvégezni a sorok közti műveleteket. Célunk, hogy a 𝑝 elem alatti elemeket kinullázzuk, hisz egy felső háromszögmátrixot szeretnénk kapni. Ehhez a következő sorműveleteket hajtjuk végre az 1. lépésben: 2
a 2. sorból kivonjuk az 1. sor
–szeresét
2
a 3. sorból kivonjuk az 1. sor –szeresét
3
3
𝑝
𝑝
𝑝
𝑝
= =
𝑎 21 𝑝 𝑎 31 𝑝
2
, azaz: 2. sor – 1. sor ∙ , 1 3
, azaz: 3. sor – 1. sor ∙ . 1
Így az 1. lépés után ezt az eredményt kapjuk: 1 −1 2 5 0 0 −3 −9 0 1 1 5 A felső háromszögmátrixos alak elérése érdekében felcseréljük a 2. és 3. sort. 1 −1 2 5 0 1 1 5 0 0 −3 −9 Ezek után a 2. lépésben ismét meghatározzuk a 𝑝 elemet, ami most a főátló második eleme lesz, tehát a 2,2 indexelésű elem, azaz 𝑝 = 1. A 𝑝 alatt 0 található, így nem kell már sorok közti műveleteket végeznünk a kinullázáshoz. Áttérhetünk a 3. lépésre, viszont ott már nincs teendőnk, mivel a felső háromszög alak előállt az együtthatómátrixban. Utolsó lépésben csak vissza kell helyettesítenünk az egyenletekbe, lentről fölfelé haladva. Ezzel megkapjuk a megoldást: 𝑥1 = 1, 𝑥2 = 2, 𝑥3 = 3. A következő tételek segítségével meghatározzuk, mikor alkalmazható a Gauss-elimináció.
2.1.1. Tétel. 2 A Gauss-módszer pontosan akkor hajtható végre, ha az 𝐴 mátrix egyik főminorja sem zérus, azaz 𝑑𝑒𝑡 𝐴 1:𝑘, 1: 𝑘
≠ 0 (𝑘 = 1, … , 𝑛).
A tétel bizonyításában 𝐴(1) a kiinduló együtthatómátrixot, 𝐴(2) pedig az első lépés után kapott együtthatómátrixot jelöli (a többi hasonlóan 𝐴 együtthatómátrixok elemei.
16
𝑛
(1)
(2)
(𝑛)
-ig). Az 𝑎 𝑖𝑗 , 𝑎 𝑖𝑗 ,… , 𝑎 𝑖𝑗
pedig ezen
Bizonyítás. 2 A Gauss-elimináció során az egyes sorokból kivonjuk más sorok számszorosait. Ez az eljárás nem változtatja meg a determinánst. Tehát det 𝐴 1:1,1: 1 det 𝐴 1: 2,1: 2
1
= det 𝐴 2
= det 𝐴
1:1,1: 1
1:2,1: 2
(1)
= 𝑎 11 ≠ 0, (1) (2)
= 𝑎 11 𝑎 22 ≠ 0,
⋮ det 𝐴 1: 𝑛, 1: 𝑛
= det 𝐴
𝑛
1:𝑛, 1: 𝑛
(1)
(𝑛)
= 𝑎 11 … 𝑎 𝑛𝑛 ≠ 0. (𝑛)
Az utolsó feltétel a visszahelyettesítés miatt kell, hiszen ez úgy kezdődik, hogy 𝑎 𝑛𝑛 -nel osztanunk kell az 𝑥𝑛 ismeretlen kiszámításához. Ebből következik az állítás.
Az alábbi tétel szintén arra ad feltételt, hogy mikor alkalmazható a Gauss-módszer.
2.1.2. Tétel. 5 Ha az 𝐴 mátrix 1. szimmetrikus, pozitív definit mátrix (SZPD) vagy 2. szigorúan diagonálisan domináns (SZDD) vagy 3. M-mátrix, akkor a Gauss-módszer végrehajtható.
2.2. LU-felbontás Amint azt az előző fejezetben láthattuk, egy mátrix Gauss-eliminációval felső háromszögmátrix alakra hozható. Ezt kibővítve kapjuk az alábbi tételt.
2.2.1. Tétel. 3 Legyen 𝐴 = 𝑎 𝑖𝑗 , 𝑖, 𝑗 = 1, … , 𝑛 egy olyan 𝑛 × 𝑛-es mátrix, melynek minden 𝑘-ad rendű főminorja reguláris. Ekkor létezik olyan 𝐿 normált (főátlóiban 1-esek állnak) alsó háromszögmátrix és egy 𝑈 felső háromszögmátrix, melyre 𝐴 = 𝐿𝑈.
17
Megjegyzés. Ha egy reguláris mátrixnak létezik LU-felbontása, akkor az LU-felbontása egyértelmű. A felbontásban szereplő 𝑈 felső háromszögmátrix megegyezik a Gauss-elimináció során kapott felső háromszögmátrixszal, az 𝐿 alsó háromszögmátrixot pedig az elimináció során használt skalárszorzók segítségével tudjuk meghatározni. Szemléltetésül az alábbi példában megadom az 𝐴 mátrix LU-felbontását. (A példa 6 -ból származik.) 𝐴=
2 1 1 4 1 0 −2 2 1
Gauss-elimináció: 2 1 1 2 1 1 2 1 1 4 1 0 ⟶ 0 −1 −2 ⟶ 0 −1 −2 . −2 2 1 0 3 2 0 0 −4 2 1 1 Tehát 𝑈 = 0 −1 −2 , a kapott felső háromszögmátrix. Az 𝐿-et a lépésenkénti 0 0 −4 skalárszorzókból kaphatjuk meg. Ismét kijelöljük a 𝑝 elemet, akárcsak a Gauss-eliminációnál. Az első lépésben 𝑝 = 2 (= 𝑎 11) és a következő műveleteket végeztük el:
(magát az 1. sort változatlanul hagytuk)
a 2. sorból kivontuk az 1. sor -szeresét, azaz: 2. sor-1. sor∙ = 2. sor −2∙1.sor
a 3. sorból kivontuk az 1. sor
4
4
𝑝
2
−2
−2
𝑝
2
-szeresét, azaz: 3. sor-1. sor∙
= 3.sor +1∙1.sor.
0 Ezek alapján határozzuk meg az első lépés utáni 𝑣1 vektort: 𝑣1 = −2 . A vektorban lévő 1 értékek azokat a műveleteket és skalárszorosokat jelölik, amelyeket az adott sor Gausseliminálásához alkalmaztunk. A 0 jelöli az első sor változatlanul hagyását, a −2-t úgy kaptuk, hogy a 2. sorhoz viszonyítva nézzük az első sorral végzett műveletet, (mivel a 2. sorból kivontuk az 1. sor kétszeresét), az 1-t pedig hasonlóan, a harmadik sorhoz viszonyítva az első sorral végzett műveletet jelöli (3. sorhoz hozzáadtuk az első sor egyszeresét). A második lépésben 𝑝 = −1(= 𝑎 22 ), és a következő műveleteket végeztük el:
1. sort változatlanul hagytuk
2. sort változatlanul hagytuk
3. sorból kivonjuk a 2. sor -szeresét, azaz: 3. sor-2. sor∙
3
3
𝑝
−1
18
=3. sor +3∙2. sor.
0 Így ezek alapján: 𝑣2 = 0 , mivel az első két sort változatlanul hagytuk, a harmadik sorhoz 3 pedig hozzáadtuk a második sor háromszorosát. E két vektorból kapjuk a következő mátrixokat: 1 0 0 1 0 0 𝐿1 𝑣1 = −2 1 0 , 𝐿2 𝑣2 = 0 1 0 . 1 0 1 0 3 1 Tehát a 𝑣1 , 𝑣2 vektorokat beleírtuk egy-egy normált alsó háromszögmátrixba. Ezekből 𝐿-et a következő képlettel kaphatjuk meg: 𝐿 = 𝐿1 𝑣1 Azaz, kiszámoljuk az 𝐿 1 𝑣1
−1,
𝐿2 𝑣2
−1
−1 ∙
𝐿 2 𝑣2
−1 .
inverzeket, majd szorzatukból megkapjuk az 𝐿 alsó
háromszögmátrixot. 1 0 0 1 0 0 2 1 0 , 𝐿 2 𝑣2 −1 = 0 1 0 −1 0 1 0 −3 1 1 0 0 1 0 0 1 0 0 𝐿= 2 1 0 0 1 0 = 2 1 0 −1 0 1 0 −3 1 −1 −3 1 1 0 0 Tehát 𝐿 = 2 1 0 . Most már ellenőrizhetjük az 𝐴 = 𝐿𝑈 felbontás helyességét. −1 −3 1 1 0 0 2 1 1 2 1 1 𝐿𝑈 = 2 1 0 0 −1 −2 = 4 1 0 −1 −3 1 0 0 −4 −2 2 1 𝐿1 𝑣1
−1
=
Ennek a mátrixszorzásnak az eredménye tényleg megegyezik 𝐴-val. Az LU-felbontás több szempontból is hasznos, például megkönnyíti a mátrix determinánsának kiszámolását. Ugyanis az 𝐴 mátrix determinánsa megegyezik az 𝑈 determinánsával (det 𝐴 = det 𝑈), melyet könnyedén megkapunk 𝑈 főátlóbeli elemeinek szorzataként. Másrészt inverz számításnál is sokkal kevesebb műveletet igényel, ha először meghatározzuk az LU-felbontást, majd megoldjuk a lineáris egyenletrendszereket.
19
2.3. Cholesky-felbontás A Cholesky-módszer az LU-faktorizációhoz hasonló felbontás, viszont ebben az esetben az 𝐴 mátrix szimmetrikus, pozitív definit. Ha a mátrix SZPD, úgy kevesebb műveletet kell elvégeznünk a faktorizáció során. Ehhez kapcsolódóan két tételt ismertetek, melyek közül az elsőnél csak a szimmetrikusság feltétele is elegendő.
2.3.1. Tétel. 2 Szimmetrikus 𝐴 ∈ ℝ𝑛×𝑛 mátrix esetén egyértelműen létezik egy 𝐿 normált alsó háromszögmátrix és egy 𝐷 diagonális mátrix, melyekkel 𝐴 = 𝐿𝐷𝐿𝑇 .
2.3.2. Tétel. (Cholesky-felbontás) 2 Tegyük fel, hogy 𝐴 egy szimmetrikus, pozitív definit mátrix. Ekkor létezik pontosan egy olyan pozitív diagonálisú 𝐺 alsó háromszögmátrix, mellyel 𝐴 = 𝐺𝐺 𝑇. Tekintsünk egy példát a felbontásra 6 . Kétféle módszerrel is megoldható a Choleskyfelbontás: LU-felbontásból kiindulva, vagy gyökvonásos módszerrel. Most a gyökvonásos algoritmusra mutatok példát. Legyen 𝐴 az alábbi SZPD mátrix, ami egyenlő egy alsó háromszögmátrix és egy felső háromszögmátrix szorzatával, mivel ezek jelölik 𝐺𝐺 𝑇 -t. (Egy alsó háromszögmátrix transzponáltja egy felső háromszögmátrix.) 𝑔11 0 5 0 −1 𝐴 = 0 6 1 = 𝑔21 𝑔22 𝑔31 𝑔32 −1 1 6
0 0 𝑔33
𝑔11 𝑔12 0 𝑔22 0 0
𝑔13 𝑔23 𝑔33
A fenti egyenlet alapján kifejezzük 𝐴 oszlopainak elemeit a 𝐺 és 𝐺 𝑇 elemeivel. A szimmetria miatt elég csak egy alsó vagy felső háromszög alakban kiszámolni 𝐴 elemeit. (Most az alsó háromszög alak szerint számolunk.) Fontos még megjegyezni, hogy 𝐴 szimmetriája miatt 𝑔12 = 𝑔21 , 𝑔13 = 𝑔31 és 𝑔23 = 𝑔32 . 1. oszlop: 2 ⇒𝑔 = 5, 5 = 𝑔11 11
0 = 𝑔21 𝑔11 = 𝑔21 5 ⇒ 𝑔21 = 0 ,
20
−1 = 𝑔31 𝑔11 = 𝑔31 5 ⇒ 𝑔31 =
−1 5
=−
5 5
,
2. oszlop: 2 2 2 2 6 = 𝑔21 𝑔12 + 𝑔22 = 𝑔21 + 𝑔22 = 0 + 𝑔22 ⇒ 𝑔22 = 6 ,
1 = 𝑔31 𝑔12 + 𝑔32 𝑔22 = 𝑔31 𝑔21 + 𝑔32 𝑔22 = 𝑔32 6 ⇒ 𝑔32 =
1 6
=
6 6
,
3. oszlop: 1
1
5
6
2 = 𝑔2 + 𝑔2 + 𝑔2 = + + 𝑔2 6 = 𝑔31 𝑔13 + 𝑔32 𝑔23 + 𝑔33 31 32 33 33 ⇒ 𝑔33 =
Ezek alapján 𝐺 =
5 0 −
5 5
0 6
0 0
6
13 30
6
5 , és 𝐺 𝑇 =
0 0
30
0 6 0
−
13 30
=
13 30 30
.
5 5 6
6
, szorzatuk pedig
13 30 30
tényleg megegyezik az 𝐴 mátrixszal. Ha a megoldáshoz a másik módszert választjuk, akkor az LU-felbontásból indulunk ki, tehát ismerjük az 𝑈 és 𝐿 mátrixokat. A szimmetrikus mátrixokra vonatkozó tétel képletével felírjuk a felbontást: 𝐴 = 𝐿𝐷𝐿𝑇 , ahol 𝐷 = 𝑑𝑖𝑎𝑔(𝑈). Ezután ugyanúgy az 𝐴 = 𝐺𝐺 𝑇 felbontást használjuk, amit a 𝐺 = 𝐿𝐷1/2 segítségével fejezünk ki.
21
3. Legkisebb négyzetek módszere
A legkisebb négyzetek módszerének létrejötte azon alapszik, hogy szükség volt egy általános megoldási módszerre az 𝐴𝑥 = 𝑏 egyenletrendszer megoldásához abban az esetben, amikor nincs az egyenletrendszernek klasszikus értelemben vett megoldása, vagyis, amikor 𝑏 nincs benne 𝐴 értékkészletében. (A fejezetet 3 alapján dolgoztam ki.) Tehát, a cél egy olyan 𝑥 vektort találni, amely az 𝐴𝑥-hez a legközelebb eső vektor. Többféle norma áll rendelkezésünkre az 𝐴𝑥 és 𝑏 közötti távolság megadásához, de a legegyszerűbb választás (ami miatt a módszer neve, hogy „legkisebb négyzetek”) az euklideszi vektornorma. Más szóval, a legkisebb négyzetek módszer lényege, hogy megtalálja a megoldást, 𝑥 ∈ ℝ𝑘 -t, a következő minimalizálási problémára: 𝑏 − 𝐴𝑥
𝑛
= min𝑦 ∈ ℝ𝑘 𝑏 − 𝐴𝑦
𝑛
,
ahol 𝐴 ∈ ℝ𝑛×𝑘 egy 𝑛 sorból és 𝑘 oszlopból álló mátrix, 𝑏 egy ℝ𝑛 -beli vektor és ∙ az euklideszi normát ℝ𝑛 -ben ( 𝑎 =
(3.1) 𝑛
jelöli
𝑎, 𝑎 ).
Mint tudjuk, négyzetes mátrix esetén 𝑛 = 𝑝, és emellett, ha 𝐴 reguláris, akkor létezik egy egyértelmű minimalizáló 𝑥 = 𝐴−1 𝑏 vektor, és a minimum egyenlő nullával. Ilyen esetben a legkisebb négyzetek módszere gyakorlatilag egyet jelent a lineáris egyenletrendszer megoldásával. Ha 𝐴 nem reguláris vagy 𝑛 ≠ 𝑝, akkor ez esetben a legkisebb négyzetek módszer általános módszert ad a lineáris egyenletrendszerek megoldására, nem négyzetes vagy szinguláris mátrixokra nézve. Ha az 𝐴𝑥 = 𝑏 egyenletrendszernek létezik megoldása, akkor ez a megoldása a legkisebb négyzetek módszernek is. Ennek (az állításnak) a megfordítása nem igaz, ezt láthatjuk is az alábbi geometriai megközelítésből. A legkisebb négyzetek módszere geometriai szempontból is értelmezhető: meg akarjuk találni a 𝑏 vektor ortogonális vetületét az 𝐴 értelmezési tartományán. Ez a meghatározás helyes, mivel 𝐴𝑥 a legközelebbi vektor 𝑏-hez 𝐼𝑚(𝐴)-ban. Az ortogonális vetítés egy jól ismert tulajdonsága, hogy 𝑏 − 𝐴𝑥 ortogonális 𝐼𝑚(𝐴)-ra. A 3.1. ábrán látható a 𝑏 vektor és 𝐴𝑥 ortogonális vetítése a vektor alterére, 𝐼𝑚(𝐴)-ra. Ennek következtében világos, hogy (3.1) mindig ad legalább egy megoldást, viszont az 𝐴𝑦 = 𝑏 egyenletrendszernek nincs megoldása, ha 𝑏 nincs benne 𝐼𝑚(𝐴)-ban.
22
3.1.ábra: Legkisebb négyzetek módszere: 𝑏 vetítése 𝐼𝑚(𝐴)-ra. 3 Fontos még megemlíteni, hogy a másik fő alkalmazása a legkisebb négyzetek módszerének az adatillesztés.
3.1. Főbb eredmények Továbbra is tekintsük a (3.1)-ben található legkisebb négyzetek problémát: keressük azt az 𝑥 ∈ ℝ𝑘 megoldást, ami minimalizálja a 𝑏 − 𝐴𝑦 oszlopból álló mátrix, és 𝑏 ∈ ℝ𝑛 , ∙
𝑛
𝑛
normát, ahol 𝐴 ∈ ℝ𝑛×𝑘 egy 𝑛 soros 𝑘
pedig az euklideszi normát jelöli.
3.1.1. Lemma. 𝟑 Egy 𝑥 ∈ ℝ𝑘 vektor megoldása a legkisebb négyzetek módszerének akkor és csak akkor, ha kielégíti az úgynevezett normálegyenletet: 𝐴∗ 𝐴𝑥 = 𝐴∗𝑏,
(3.2)
ahol 𝐴∗ az 𝐴 transzponált konjugáltja, és 𝐴∗ 𝐴 egy 𝑘 méretű négyzetes mátrix.
Bizonyítás. 𝟑 Legyen 𝑥 ∈ ℝ𝑘 a (3.1) megoldása. Például: 𝑏 − 𝐴𝑥
2 𝑛
≤ 𝑏 − 𝐴𝑦
2 𝑛
,
∀𝑦 ∈ ℝ𝑘 .
Bármely 𝑧 ∈ ℝ𝑘 -ra és bármely 𝑡 ∈ ℝ-re legyen 𝑦 = 𝑥 + 𝑡𝑧 . Ekkor 𝑏 − 𝐴𝑥
2 𝑛
≤ 𝑏 − 𝐴𝑥
2 𝑛
+ 2𝑡 𝐴𝑥 − 𝑏 , 𝐴𝑧 + 𝑡 2 𝐴𝑧
Ebből az alábbira következtetünk:
23
2. 𝑛
0 ≤ 2 sgn(𝑡) 𝐴𝑥 − 𝑏 , 𝐴𝑧 + 𝑡 𝐴𝑧
2 𝑛
,
ami arra utal, hogy ha t tart a 0-hoz, akkor 𝐴𝑥 − 𝑏, 𝐴𝑧 = 0 ∀𝑧 ∈ ℝ𝑘 . Ezután arra jutunk, hogy 𝐴∗ 𝐴𝑥 − 𝐴∗𝑏 = 0. Ezt megfordítva, ha 𝑥 a (3.2)-beli normálegyenlet megoldása, akkor 𝐴𝑥 − 𝑏, 𝐴𝑧 = 0,
∀𝑧 ∈ ℝ𝑘 .
Tehát 𝑏 − 𝐴𝑥
2 𝑛
≤ 𝑏 − 𝐴𝑦
2 𝑛
∀𝑦 = 𝑥 + 𝑡𝑧 ∈ ℝ𝑘 .
Vagyis, 𝑥 megoldása (3.1)-nek is.
3.1.2. Tétel. 𝟑 Minden 𝐴 ∈ ℝ𝑛×𝑘-s mátrixra a (3.2)-beli normálegyenletnek mindig létezik legalább egy megoldása. Továbbá, ez a megoldás egyértelmű akkor és csak akkor, ha 𝐾𝑒𝑟(𝐴) = 0 . A normálegyenlet partikuláris megoldása kifejezhető 𝐴 egy pszeudoinverzének (𝐴† ) elemeivel.
3.1.3. Tétel. 𝟑 Az 𝑥𝑏 = 𝐴† 𝑏 vektor megoldása a (3.1)-ben definiált legkisebb négyzetek módszernek. Ha (3.1)-nek több megoldása is van, akkor közülük 𝑥𝑏 az az egyértelmű megoldás, ami a legkisebb normával rendelkezik. Például, ∀𝑥 ≠ 𝑥𝑏 -re, melyre 𝐴𝑥𝑏 − 𝑏
2
= 𝐴𝑥 − 𝑏 2, a
következőt kapjuk 𝑥𝑏
2
< 𝑥 2.
Bizonyítás. 𝟑 Minden 𝑥 ∈ ℝ𝑘 -re 𝐴𝑥 − 𝑏 felbontható a következőképpen: 𝐴𝑥 − 𝑏 = 𝐴 𝑥 − 𝑥𝑏 − 𝐸 − 𝐴𝐴† 𝑏. Ez a felbontás ortogonális, mivel 𝐴 𝑥 − 𝑥𝑏 ∈ 𝐼𝑚(𝐴) és 𝐸 − 𝐴𝐴† 𝑏 ∈ 𝐼𝑚(𝐴)⊥ , mert 𝐴𝐴† az ortogonális vetítés mátrixa ℂ𝑚 -nek 𝐼𝑚(𝐴)-ra. Ebből a felbontásból arra következtethetünk, hogy 𝐴𝑥 − 𝑏
2 2
= 𝐴𝑥 − 𝐴𝑥𝑏
2 2
+ 𝐴𝑥𝑏 − 𝑏
24
2 2
2
≥ 𝐴𝑥𝑏 − 𝑏 2 ,
(3.3)
ami alapján láthatjuk, hogy 𝑥𝑏 valóban megoldása (3.1)-nek. Továbbá, ha
𝐴𝑥𝑏 − 𝑏
2
=
𝐴𝑥 − 𝑏 2, akkor (3.3)-mal látható, hogy 𝐴𝑥 = 𝐴𝑥𝑏 és 𝑧 = 𝑥 − 𝑥𝑏 ∈ 𝐾𝑒𝑟(𝐴). Ezzel a következő felbontást kapjuk 𝑥-re: 𝑥 = 𝑧 + 𝑥𝑏 . Ez a felbontás ortogonális, mivel 𝑧 ∈ 𝐾𝑒𝑟(𝐴) és 𝑥𝑏 = 𝐴† 𝑏 ∈ (𝐾𝑒𝑟(𝐴))† . Emiatt, ha 𝑥 ≠ 𝑥𝑏 , 𝑥
2 2
= 𝑧
2 2
+ 𝑥𝑏
2 2
> 𝑥𝑏
2. 2
Megjegyzés. Az 𝑥𝑏 = 𝐴† 𝑏 vektornak egyszerű geometriai értelmezése van: 𝑥𝑏 a (𝐾𝑒𝑟(𝐴))† -nak az az egyértelmű megoldás vektora, aminek az 𝐴 mátrix szerinti képe megegyezik 𝑏-nek 𝐼𝑚(𝐴)-ra való vetítésével. Végül tekintsünk egy példát a módszerre. 10 Van 3 hegy, melyeket 𝑥1 ,𝑥2 és 𝑥3-mal jelölünk. Egy helyről mérve őket rendre a következő magasságokat kaptuk: 2474 láb, 3882 láb és 4834 láb (1 láb=30,48 cm). Viszont 𝑥1-ről tekintve a másik 2 hegyre, 𝑥2 1422 lábbal, 𝑥3 pedig 2354 lábbal magasabbnak tűnik. Ha pedig 𝑥2-ről nézzük a harmadik hegyet, akkor onnan 𝑥3 950 lábbal tűnik magasabbnak. Ezek alapján felírhatjuk a kiinduló egyenletrendszert:
𝐴𝑥 = 𝑏,
Minimalizálni akarjuk az
𝐴=
1 0 0 1 0 0 −1 1 −1 0 0 −1
0 0 1 0 1 1
𝑥1 𝑥2 𝑥3
2474 3882 4834 = 1422 2354 950
=𝑏
𝐴𝑥 − 𝑏 2-t. A normálegyenlet alakjának elérése érdekében
mindkét oldalt megszorozzuk 𝐴∗-gal. 1 0 0 −1 −1 0 𝐴 = 0 1 0 1 0 −1 0 0 1 0 1 1 1 0 0 0 1 0 −1 −1 0 3 −1 −1 0 0 1 = −1 3 −1 1 0 −1 −1 1 0 0 1 1 −1 −1 3 −1 0 1 0 −1 1 2474 3882 0 0 −1 −1 0 −1302 4834 = 4354 1 0 1 0 −1 1422 0 1 0 1 1 8138 2354 950 ∗
𝐴∗ 𝐴 =
1 0 0 0 1 0 0 0 1
1 𝐴∗ 𝑏 = 0 0
25
𝐴∗𝐴𝑥 = 𝐴∗ 𝑏 3 −1 −1 𝑥1 −1302 𝑥 = −1 3 −1 2 4354 𝑥3 −1 −1 3 8138 Ebből az egyenletrendszerből kiszámítva a megoldást: 𝑥1 = 2472, 𝑥2 = 3886, 𝑥3 = 4832 értékeket kapjuk.
26
4. Iterációs módszerek Egy iterációs módszer alkalmazásánál végtelen sok lépésre van szükség ahhoz, hogy a pontos megoldáshoz jussunk. Ez úgy lehetséges, hogy az iteráció során a megoldást olyan 𝑥𝑘 konvergens sorozatokkal közelítjük, amelyek határértéke az 𝐴𝑥 = 𝑏 egyenletrendszer egyértelmű megoldása (𝑥 ∗ ), és mivel végtelen lépés szükséges, így az iterációk száma (𝑘) a +∞-be tart. Az 𝐴 mátrixot továbbra is négyzetesnek és regulárisnak tekintjük. A lineáris iterációs módszerek általános alakja a következőképp írható fel: 𝑥 (𝑘+1) = 𝐵𝑥 (𝑘) + 𝑓 , ahol a 𝐵 mátrixot iterációs mátrixnak nevezzük és 𝑥 (𝑖) az 𝑥 vektor 𝑖-edik iterációját jelöli. Az alábbi definíció a konzisztencia fogalmát ismerteti.
Definíció. Az 𝑥 (𝑘+1) = 𝐵𝑥 (𝑘) + 𝑓 iterációt az 𝐴𝑥 = 𝑏 egyenletrendszerrel konzisztensnek hívjuk, ha 𝑥 ∗ = 𝐵𝑥 ∗ + 𝑓 . (Az 𝑥 ∗ az egyenletrendszer megoldása.) Egy iterációs módszer a következő felbontáson alapul: 𝐴 = 𝑀 − 𝑁, ahol 𝑀 reguláris. Ekkor az iterációs eljárás az alábbi módon definiálható: ( 3 alapján) 𝑎𝑑𝑜𝑡𝑡 𝑥(0) ∈ ℝ𝑛 𝑘𝑒𝑧𝑑𝑒𝑡𝑖 𝑣𝑒𝑘𝑡𝑜𝑟 𝑀𝑥 (𝑘+1) = 𝑁𝑥 (𝑘) + 𝑏 ∀𝑘 ≥ 1.
(4.1)
Az iteráció konvergenciájára a következő tétel ad feltételt.
4.1.1. Tétel. 3 A (4.1)-ben definiált iterációs módszer akkor és csak akkor konvergens, ha az 𝑀 −1 𝑁 mátrix spektrálsugara kielégíti a következő feltételt: 𝜌(𝑀 −1 𝑁) < 1. A tétel bizonyításához szükséges az alábbi definíció:
Definíció. Az 𝑟 (𝑘) = 𝑏 − 𝐴𝑥 (𝑘) vektort (illetőleg az 𝑒 (𝑘) = 𝑥 (𝑘 ) − 𝑥 -et) a 𝑘-adik iteráció hibájának nevezzük. Nyilvánvaló, hogy az iteráció akkor konvergens, ha 𝑒 (𝑘) a 0-hoz tart.
27
A 4.1.1. tétel bizonyítása:
Bizonyítás. 3 Az 𝑒 (𝑘) hiba megadható indukciós eljárással a következőképpen: 𝑒 (𝑘) = 𝑥 (𝑘) − 𝑥 = 𝑀−1 𝑁𝑥 (𝑘−1) + 𝑀 −1 𝑏 − 𝑀−1 𝑁𝑥 + 𝑀−1 𝑏 = 𝑀 −1 𝑁 𝑥(𝑘−1) − 𝑥 = = 𝑀 −1 𝑁𝑒 (𝑘−1) . Mivel 𝑒 (𝑘) = (𝑀−1 𝑁)𝑘 𝑒 (0), és célunk, hogy a hiba 0-hoz tartson, ezért lim𝑘⟶+∞ 𝑒 (𝑘) = 0 minden 𝑒 (0)-ra, ami akkor és csak akkor teljesül, ha 𝜌(𝑀−1 𝑁) < 1.
Tehát ha a megoldást közelítő sorozatok konvergálnak egy határhoz, pontosabban 𝑥 ∗-hoz, és az iterációk száma tart a végtelenbe, akkor a következőre jutunk: 𝑀 − 𝑁 𝑥 = 𝐴𝑥 = 𝑏. Továbbá a (4.1)-ben definiált sorozattal ekvivalens a következő felírás: 𝑥 (𝑘 +1) = 𝑀 −1 𝑁𝑥 (𝑘) + 𝑀 −1 𝑏. Ez esetben az 𝑀−1 𝑁 mátrixot nevezzük az iterációs módszer iterációs mátrixának. Könnyen láthatjuk, hogy ez a felírás megegyezik az általános alakban felírt egyenlettel, ha 𝐵 = 𝑀−1 𝑁 és 𝑓 = 𝑀 −1 𝑏. Az alfejezetekben a klasszikus iterációs módszereket mutatom be példákkal és összehasonlítást is végzek a módszerek gyorsasága között.
4.1. Jacobi-iteráció Az előző 𝑥 (𝑘+1) = 𝑀−1 𝑁𝑥 (𝑘) + 𝑀 −1 𝑏 alakból az 𝑀 = 𝐷 és 𝑁 = 𝐷 − 𝐴 helyettesítéssel kaphatjuk meg a Jacobi-iteráció alakját: 𝑥 (𝑘+1) = 𝐷−1 𝐷 − 𝐴 𝑥 (𝑘) + 𝐷−1 𝑏
⇒
𝑥 (𝑘+1) = (𝐸 − 𝐷−1 𝐴)𝑥(𝑘) + 𝐷−1 𝑏.
Itt 𝐸 − 𝐷−1 𝐴 jelöli a Jacobi-iteráció iterációs mátrixát, melyet 𝐵𝐽 -vel fogunk jelölni. A felbontásban szereplő 𝐷 az 𝐴 együtthatómátrix diagonális mátrixát jelöli. Itt 𝐴 az 𝐴 = 𝐿 + 𝐷 + 𝑈 alakban írható fel, ahol 𝐿 és 𝑈 nem az LU-felbontás mátrixai, hanem 𝐿 a diagonális alatti elemek -1-szereseinek, 𝑈 pedig a diagonális fölötti elemek -1-szereseinek mátrixa. A
28
számolásokat a Jacobi-iteráció komponensenkénti alakjával érdemesebb elvégezni, ami a következőképpen írható fel: 𝑥𝑖
𝑘+1
1 =− 𝑎 𝑖𝑖
𝑛
𝑎 𝑖𝑗 𝑥𝑗 𝑘 − 𝑏 𝑖 , 𝑖 = 1, … , 𝑛. 𝑗 =1,𝑗≠𝑖
Példaként tekintsük az alábbi egyenletrendszert, és végezzünk el egy lépést a Jacobiiterációval. (A kiinduló adatok a 11 -ből származnak.) 5 1 0 𝐴= 1 2 0 , 0 1 3 Először kiszámítjuk a 𝐵𝐽 iterációs mátrixot.
7 1 𝑏 = 5 , 𝑥 (0) = 1 2 1
1 0 0 5 −1 𝐵𝐽 = 𝐸 − 𝐷 𝐴 = 0 1 0 − 0 0 0 1 0 0 −1/5 −1/2 0 𝐵𝐽 = 0 −1/3
0 0 2 0 0 3 0 0 0
−1
5 1 0 1 2 0 0 1 3
Ezután ellenőrizzük a konvergencia feltételét, tehát a 𝜌(𝐵𝐽 ) < 1 feltételt. Ehhez meg kell határoznunk 𝐵𝐽 sajátértékeit, ugyanis ezek közül abszolút értékben a legnagyobbra lesz szükségünk, mivel 𝜌 𝐵𝐽 = max λ (ahol 𝜆 a 𝐵𝐽 mátrix sajátértékét jelöli). −𝜆 −1/5 0 1 𝜆 0 = −𝜆3 + 0 + 0 − 0 − 0 + 𝜆 = −𝜆3 + det(𝐵𝐽 − 𝜆𝐸) = det −1/2 −𝜆 10 10 0 −1/3 −𝜆 𝜆 𝜆 1 10 10 =0 ⇒ = 𝜆3 → 𝜆1 = 0 ⇒ = 𝜆2 → 𝜆2 = , 𝜆3 = − 10 10 10 10 10
⇒ −𝜆3 +
Tehát 𝜌 𝐵𝐽 =
10 10
≈ 0,316 < 1, azaz az iteráció konvergens. Most már elkezdhetjük az 𝑥 (1)
vektor kiszámítását. (1)
𝑥1
𝑥 (1) = 𝑥2(1)
és a szükséges képlet 𝑥𝑖 𝑘+1 = −
(1)
𝑥3 Így: 1 5 1 (1) 𝑥2 = − 2 1 (1) 𝑥3 = − 3 𝑥11 = −
1∙ 1 +0 ∙ 1 −7 =
6 5
1∙ 1 +0 ∙ 1 −5 = 2 0∙ 1 +1 ∙ 1 −2 =
1 3
29
1 𝑎 𝑖𝑖
𝑘 𝑛 𝑗 =1,𝑗 ≠𝑖 𝑎 𝑖𝑗 𝑥𝑗
− 𝑏𝑖 .
Tehát egy lépés után az
𝑥 (1)
6/5 = 2 vektort kapjuk. 1/3
A következő fejezetben nézünk egy példát a Gauss-Seidel-iterációval is. Következtetéseket csak ezután vonunk le.
4.2. Gauss-Seidel-iteráció Az 𝑀 = 𝐷 − 𝐿 és 𝑁 = 𝑈 helyettesítéssel kapjuk meg a Gauss-Seidel-iteráció alakját: 𝑥 (𝑘 +1) = (𝐷 − 𝐿)−1 𝑈𝑥 (𝑘) + (𝐷 − 𝐿)−1 𝑏. Ezesetben a (𝐷 − 𝐿)−1 𝑈 mátrixot nevezzük iterációs mátrixnak és 𝐵𝐺𝑆 -sel jelöljük. Az iteráció komponensenkénti alakja a következő: (𝑘+1) 𝑥𝑖
1 =− 𝑎 𝑖𝑖
𝑖−1
𝑛 (𝑘+1) 𝑎 𝑖𝑗 𝑥𝑗
𝑗 =1
(𝑘)
+
𝑎 𝑖𝑗 𝑥𝑗
− 𝑏𝑖 .
𝑗 =𝑖+1 (𝑘+1)
A fő különbség a Jacobi és a Gauss-Seidel iterációk között, hogy az utóbbi 𝑥𝑖
megadásához felhasználja a már addig kiszámolt 1, … , 𝑖 − 1 indexű 𝑥 (𝑘+1) elemeket, míg a Jacobi-iteráció csak az 𝑥 (𝑘) vektort használja fel 𝑥 (𝑘+1) elemeinek megadásához. Tekintsük most az alábbi példát ( 9 ) és végezzünk el egy lépést Gauss-Seidel-iterációval. 2 1 0 𝐴= 0 4 2 , 0 0 1 𝐵𝐺𝑆 = (𝐷 − 𝐿)−1 𝑈 =
1 1 𝑏 = 1 , 𝑥 (0) = 2 1 0
2 0 0 0 0 0 − 0 4 0 0 0 0 0 0 1 0 0 0
−1
0 1 0 0 0 2 0 0 0
0 1/2 0 𝐵𝐺𝑆 = 0 0 1/2 0 0 0 Most is ellenőriznünk kell a konvergencia feltételét. −𝜆 1/2 0 det(𝐵𝐺𝑆 − 𝜆𝐸) = det 0 −𝜆 1/2 = −𝜆 𝜆2 − 0 = −𝜆3 0 0 −𝜆 ⇒ −𝜆3 = 0 ⟶ 𝜆 = 0
30
Tehát 𝜌 𝐵𝐺𝑆 = 0 < 1, azaz az iteráció konvergens minden kezdővektorra. Továbbá ha a spektrálsugár 0-val egyenlő, akkor az iterációról tudjuk, hogy véges. Így elkezdhetjük az első lépés kiszámítását. (1)
𝑥1
(𝑘 +1)
𝑥 (1) = 𝑥2(1) , szükséges képlet: 𝑥𝑖
=−
(1)
1 𝑎 𝑖𝑖
(𝑘 +1) 𝑖−1 + 𝑗 =1 𝑎 𝑖𝑗 𝑥𝑗
(𝑘 ) 𝑛 𝑗 =𝑖+1 𝑎 𝑖𝑗 𝑥𝑗
− 𝑏𝑖
𝑥3 1 2 1 (1) 𝑥2 = − 4 𝑥11 = −
(1)
𝑥3 = −
1 1
1∙ 2 +0 ∙ 0 −1 = −
1 2
0∙ −
1 1 +2∙0 −1 = 2 4
0∙ −
1 1 +0∙ −1 =1 2 4
−1/2 1/4 . 1 Összehasonlítva az előző fejezetből a Jacobi-iteráció, és a most számolt Gauss-Seidel-
Tehát Gauss-Seidel-iterációval egy lépés után 𝑥 (1) =
iteráció spektrálsugarait, láthatjuk, hogy a Gauss-Seidel-iterációhoz tartozó a kisebb. Minél kisebb a spektrálsugár, annál gyorsabban tart az iteráció a megoldáshoz. Tehát ez esetben a 2. módszer, azaz a Gauss-Seidel-iteráció konvergál gyorsabban a megoldáshoz. Ez a megállapítás azonban nem általánosítható, tudunk olyan példát is mutatni, amelyben a Jacobiiteráció konvergál gyorsabban. A konvergencia vizsgálatához szükséges és elégséges feltételt adott az iterációs mátrix spektrálsugarára szabott feltétel ellenőrzése. Ugyanakkor
elégséges feltételt ad a
konvergenciára, ha az iterációs mátrix valamely normában kisebb, mint 1. Viszont a konvergencia feltételeinek ellenőrzését nem mindig kell elvégeznünk, ugyanis bizonyos típusú együtthatómátrixokra a konvergencia adott. Így csak azt kell megállapítani, hogy 𝐴 tényleg az éppen adott tulajdonságú együtthatómátrixnak felel-e meg. E szempontból a következő esetek lehetségesek.
4.2.1. Tétel. 2 Ha az egyenletrendszer együtthatómátrixa M-mátrix, akkor a Jacobi és a Gauss-Seideliterációk konvergálnak az egyenletrendszer megoldásához tetszőleges kezdeti vektor esetén.
31
4.2.2. Tétel. 2 Szigorúan diagonálisan domináns (SZDD) együtthatómátrixok esetén a Jacobi és GaussSeidel-iterációk minden kezdővektor esetén az egyenletrendszer megoldásához konvergálnak.
4.2.3. Tétel. 2 Ha az 𝐴 együtthatómátrix szimmetrikus és pozitív definit (SZPD), akkor a Gauss-Seideliteráció konvergál az egyenletrendszer megoldásához tetszőleges kezdeti vektor esetén.
Fontos még meghatározni, hogy meddig tart egy iterációs eljárás, hány lépés után kell leállnunk. Erre ad feltételt a Banach-féle fixponttétel: 𝑥
(𝑘)
−𝑥
𝐿𝑘 ≤ 𝑥 (1) − 𝑥 (0) . 1−𝐿
∗
A mi esetünkben 𝐿 a 𝐵 iterációs mátrix normája, és szükséges, hogy 𝐵 < 1, azaz az iterációs mátrix valamely normában kisebb legyen, mint 1. Ennek segítségével és az első iterációs lépés eredményével megadható, hogy hány lépésre van szükségünk adott pontosság eléréséhez. Tekintsünk egy példát erre: hány iterációs lépést kell megtennünk ahhoz, hogy 10−2 nagyságú hibával határozzuk meg az egyenlet megoldását? 6 Tegyük fel, hogy az alábbiakat már kiszámoltuk: 𝑥 (0) 𝐵𝐽
−1/5 −1/5 0 (1) (1) 9/4 43/20 = 0 , Jacobi-iterációval: 𝑥 = , Gauss-Seidel-iterációval: 𝑥 = 0 37/9 35/12 3
∞
= 0,75 = < 1 és 4
𝐵𝐺𝑆
2
∞
= 0,4 = < 1 5
Oldjuk meg először Jacobi-iterációval a feladatot, ekkor 𝐿 = 𝐵𝐽 3 𝑘 𝑥 (𝑘) − 𝑥 ∗ ∞ ≤ 4 3 1− 4 3 4 3 4
−1/5 0 9/4 − 0 37/9 0
𝑘
∙4∙ 𝑘
<
∞
.
< 10 −2 ∞
37 < 10−2 9
9 100 ∙ 37 ∙ 4
9 14800 ≈ 25,74 ⇒ 𝑘 = 26 𝑘> 3 log 4 Tehát 26 lépést kell megtennünk Jacobi-iterációval az adott pontosságú közelítés eléréséhez. log
32
Most oldjuk meg a feladatot Gauss-Seidel-iterációval, ekkor 𝐿 = 𝐵𝐺𝑆 2 𝑘 𝑥 (𝑘) − 𝑥 ∗ ∞ ≤ 5 2 1− 5 2 5 2 5
𝑘
−1/5 0 43/20 − 0 35/12 0
∞.
< 10 −2 ∞
5 35 ∙ ∙ < 10−2 3 12
𝑘
<
36 100 ∙ 5 ∙ 35
36 17500 ≈ 6,75 ⇒ 𝑘 = 7 𝑘> 2 log 5 pontosság mellett, a Gauss-Seidel-iterációval csak 7 lépést kell elvégeznünk a log
Így 10 −2
megoldás közelítéséhez.
4.3. Richardson-iteráció Ennél az iterációs módszernél feltesszük, hogy az együtthatómátrix szimmetrikus és pozitív definit (SZPD), valamint definiálunk egy 𝑝 paramétert, amire 𝑝 ∈ ℝ, 𝑝 ≠ 0. A Richardson iteráció alakját a lineáris egyenletrendszer általános alakjából a következő átalakításokkal kapjuk meg: 𝐴𝑥 = 𝑏 ⇒ 𝑝𝐴𝑥 = 𝑝𝑏 ⇒ 0 = −𝑝𝐴𝑥 + 𝑝𝑏 ⇒ 𝑥 = 𝑥 − 𝑝𝐴𝑥 + 𝑝𝑏 ⇒ 𝑥 = (𝐸 − 𝑝𝐴)𝑥 + 𝑝𝑏 Tehát megszoroztuk az egyenletrendszert 𝑝-vel, 0-ra rendeztünk, majd mindkét oldalhoz hozzáadtunk 𝑥-et. Vagyis a Richardson-iteráció alakja: 𝑥 (𝑘+1) = 𝐸 − 𝑝𝐴 𝑥
𝑘
+ 𝑝𝑏.
Ez esetben az iterációs mátrix 𝐸 − 𝑝𝐴 , melyet 𝐵𝑅 -rel jelölünk. A főtengelytétel alapján tudjuk, hogy ha 𝐴 szimmetrikus, akkor a sajátértékei valósak. Mivel 𝐴 nemcsak szimmetrikus, hanem pozitív definit is, ezért tudjuk, hogy sajátértékei mind pozitívak. Az iterációs mátrix sajátértékéről most teszünk egy állítást, amit szeretnénk, hogy teljesüljön, majd megpróbáljuk bizonyítani sejtésünket.
33
4.3.1. Állítás. 8 Az 𝐸 − 𝑝𝐴 mátrix sajátértéke 1 − 𝑝𝜆𝑖 .
Bizonyítás. 8 A sajátérték, sajátvektor definíciójából a következőnek kell teljesülni: 𝐴𝑣𝑖 = 𝜆𝑖 𝑣𝑖
⇒
𝐸 − 𝑝𝐴 𝑣𝑖 = 1 − 𝑝𝜆𝑖 𝑣𝑖 .
Csak a bal oldali tagot alakítjuk át: 𝐸 − 𝑝𝐴 𝑣𝑖 = 𝑣𝑖 − 𝑝𝐴𝑣𝑖 = 𝑣𝑖 − 𝑝𝜆𝑖 𝑣𝑖 = 1 − 𝑝𝜆𝑖 𝑣𝑖 . Tehát beláttuk, hogy a bal oldali kifejezés egyenlő a jobb oldalival, azaz 1 − 𝑝𝜆𝑖 tényleg
𝐸 − 𝑝𝐴 sajátértéke.
Ezzel az állítással már megkaptuk, hogy 𝜌 𝐵𝑅 = max |1 − 𝑝𝜆𝑖 | < 1. Tegyük fel, hogy 𝜆1 ≥ 𝜆2 ≥ ⋯ ≥ 𝜆𝑛 > 0, azaz ∀ 𝜆𝑖 > 0. Ekkor 𝜆1 a legnagyobb, ezért erre nézzük meg, hogy mikor teljesülne a spektrálsugárra vonatkozó becslés. −1 < 1 − 𝑝𝜆1 < 1 Rendezzük a bal oldali egyenlőtlenséget 𝑝-re. −1 < 1 − 𝑝𝜆1 ⟶ 𝑝 =
2 𝜆1
Az alábbi tétel ugyanezt a képletet adja 𝑝 -re, emellett az iteráció konvergenciájához fontos feltételt ad.
4.3.2. Tétel. 7 Legyen 𝐴 szimmetrikus, pozitív definit mátrix. Ekkor az 𝑥 (𝑘+1) = 𝐸 − 𝑝𝐴 𝑥 𝑝 ∈ 0; 𝑝𝑜𝑝𝑡 =
2 𝜌 𝐴 2
𝜆 1 +𝜆 𝑛
𝑘
+ 𝑝𝑏 iteráció
esetén tetszőleges 𝑥 (0) mellett konvergens. Ekkor az optimális paraméter és 𝜌 𝐵𝑅 (𝑝𝑜𝑝𝑡 ) =
𝜆 1 −𝜆 𝑛 𝜆 1 +𝜆 𝑛
, ahol 𝜌 𝐴 = 𝜆1 ≥ 𝜆2 ≥ ⋯ ≥ 𝜆𝑛 > 0 az 𝐴 mátrix
sajátértékei. Nézzünk meg egy példát, amiben az optimális 𝑝 paramétert keressük. 9 𝐴=
2 −1 −1 2
Látható, hogy 𝐴 szimmetrikus, mivel 𝐴 = 𝐴𝑇 , továbbá a mátrixnak egyik főminorja sem zérus, ezért 𝐴 pozitív definit. Tehát 𝐴 SZPD mátrix, így alkalmazható a Richardson iteráció. Először meg kell határoznunk az együtthatómátrix spektrálsugarát, hogy tudjuk, 𝑝-t mely 34
intervallumból kell választanunk. A spektrálsugár kiszámításához az 𝐴 karakterisztikus polinomja szükséges, hogy megkapjuk a gyököket. 𝑝𝐴 𝜆 = (2 − 𝜆)2 − 1 = 4 − 4𝜆 + 𝜆2 − 1 = 𝜆2 − 4𝜆 + 3 𝜆2 − 4𝜆 + 3 = 0 ⟶ 𝜆1 = 3, Tehát 𝜌 𝐴 = 3, így 𝑝 ∈ 0;
2 3
. Az optimális 𝑝: 𝑝𝑜𝑝𝑡 =
35
𝜆2 = 1 2
3+1
1
= . 2
Összefoglalás Dolgozatomban ismertettem a lineáris egyenletrendszerek két fontosabb numerikus megoldási módját, a direkt és iterációs módszereket, valamint a legkisebb négyzetek módszerét. Direkt módszereken belül megvizsgáltuk a Gauss-elimináció, LU-felbontás és a Cholesky-felbontás alkalmazását, az iterációs módszerek közül pedig a Jacobi-, Gauss-Seidelés Richardson-iterációval ismerkedtünk meg. Mindezen módszerek alkalmazásánál fontos szerepe volt a mátrixoknak, hiszen a példákban ezek segítségével kaptuk meg az adott problémára a megoldást. Emellett többször is előfordult, hogy az éppen használt módszer alkalmazhatóságához vagy a feladat megoldhatóságához egy adott mátrixtulajdonságra volt szükség. Ezek alapján is látható, hogy az algebrából ismert mátrixok az analízis területein is rendkívül jól alkalmazhatóak. Fontos még megemlíteni, hogy a szakdolgozatban említetteken kívül a lineáris egyenletrendszerek megoldásához természetesen még számtalan más módszer is létezik.
36
Köszönetnyilvánítás Ezúton szeretném megköszönni témavezetőmnek, Fialowski Alice tanárnőnek, hogy segítőkészségével,
türelmével,
hasznos
tanácsaival
hozzájárult
a
szakdolgozatom
elkészítéséhez. Továbbá köszönöm családom bíztatását, szaktársaim segítőkészségét és mindazoknak, akik támogattak a szakdolgozat megírása során.
37
Irodalomjegyzék 1 Freud Róbert: Lineáris algebra, ELTE Eötvös Kiadó, 2006 2 Faragó István - Horváth Róbert: Numerikus módszerek, TYPOTEX, 2011 3 Grégoire Allaire, Sidi Mahmoud Kaber: Numerical Linear Algebra, Springer 4 Stoyan Gisbert, Takó Galina: Numerikus módszerek I., TYPOTEX 5 Faragó István: Alkalmazott analízis II. előadás anyaga, 6 Fekete Imre: Alkalmazott analízis II. gyakorlat anyaga 7 László Lajos: Numerikus módszerek I. előadás anyaga (IK), http://numanal.inf.elte.hu/~laszlo/ea6.pdf 8 Lócsi Levente: Numerikus módszerek I. gyakorlat anyaga (IK) 9 Bozsik József, Krebsz Anna: Numerikus módszerek példatár, ELTE IK, 2010 10 http://cs.wpunj.edu/~kaufmanl/cs480/leastsquarescorr.ppt 11 http://www.stud.u-szeged.hu/Betyar.Gabor/jacobi.pdf
38
Nyilatkozat Név: Borostyán Dóra ELTE Természettudományi Kar, szak: Matematika BSc ETR azonosító: BODSABT.ELTE NEPTUN azonosító: AZWB6V Szakdolgozat címe: Lineáris algebra és mátrixok alkalmazása a numerikus analízisben
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.
Budapest, 2013. 05.
________________________ a hallgató aláírása
39
