Masarykova univerzita Brno Fakulta pˇr´ırodovˇedeck´a Katedra aplikovan´e matematiky
Line´ arn´ı syst´ emy se speci´ aln´ımi maticemi Diplomov´a pr´ace
kvˇeten 2006
Jaroslava Ben´aˇckov´a
Podˇ ekov´ an´ı Vu ´vodu bych r´ada podˇekovala vedouc´ı diplomov´e pr´ace prof. RNDr. Ivanˇe Horov´e, CSc. z Katedry aplikovan´e matematiky PˇrF MU v Brnˇe za peˇcliv´e pˇreˇcten´ı textu, cenn´e rady, pˇripom´ınky k pr´aci a za trpˇelivost. D´ale bych chtˇela podˇekovat sv´ ym rodiˇc˚ um za veˇskerou podporu, kter´e se mi v pr˚ ubˇehu studia dostalo.
Prohl´ aˇ sen´ı ˇ Cestnˇ e prohlaˇsuji, ˇze jsem svou diplomovou pr´aci vypracovala samostatnˇe a pouˇzila jsem pouze uvedenou literaturu.
V Brnˇe dne 25. kvˇetna 2006
Obsah ´ Uvod
5
1 Z´ akladn´ı pojmy 1.1 Matice a operace s nimi . . . . . . . . 1.2 Grafy a matice . . . . . . . . . . . . . 1.3 Pˇr´ım´e metody ˇreˇsen´ı soustav line´arn´ıch 1.3.1 Gaussova eliminaˇcn´ı metoda . . 1.3.2 Metoda LU-rozkladu . . . . . . 1.3.3 Metoda nejvˇetˇs´ıho sp´adu . . . . 2 Speci´ aln´ı matice 2.1 Nez´aporn´e matice . . . . . . . . . . 2.2 Stochastick´e matice . . . . . . . . . 2.2.1 Dvojitˇe stochastick´e matice 2.3 Symetrick´e a hermitovsk´e matice . 2.4 M-matice . . . . . . . . . . . . . . 2.5 Stabiln´ı matice . . . . . . . . . . . 2.6 P´asov´e matice . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . . . rovnic . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . .
. . . . . . .
. . . . . .
. . . . . . .
3 Iteraˇ cn´ı metody ˇ reˇ sen´ı soustav line´ arn´ıch rovnic 3.1 Iteraˇcn´ı metody . . . . . . . . . . . . . . . . . . . 3.2 Jacobiho metoda (metoda prost´e iterace) . . . . . 3.3 Gauss-Seidelova metoda . . . . . . . . . . . . . . 3.4 Superrelaxaˇcn´ı metoda . . . . . . . . . . . . . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
. . . . . . .
. . . .
. . . . . .
6 6 7 9 9 12 13
. . . . . . .
17 17 24 25 27 29 31 32
. . . .
36 36 38 39 43
4 Singul´ arn´ı line´ arn´ı syst´ emy 48 4.1 Semikonvergentn´ı matice . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Semiiteraˇcn´ı metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 Z´ avˇ er
53
Literatura
54
Pˇ r´ıloha
55
4
´ Uvod V posledn´ı dobˇe st´ale roste z´ajem o ˇreˇsen´ı velk´ ych soustav line´arn´ıch algebraick´ ych soustav jak s ˇr´ıdk´ ymi, tak s pln´ ymi maticemi nebo o ˇreˇsen´ı u ´loh na vlastn´ı ˇc´ısla matic. Podstatn´ y vliv na pˇrehodnocen´ı star´ ych numerick´ ych metod line´arn´ı algebry mˇely modern´ı poˇc´ıtaˇce, kter´e podn´ıtily z´ajem o nov´e algoritmy, jeˇz se hod´ı k automatizovan´emu prov´adˇen´ı v´ ypoˇct˚ u. Z´akladn´ı skupinu metod line´arn´ı algebry tvoˇr´ı metody pˇr´ım´e. Pˇr´ımou metodou se obvykle rozum´ı metoda, kter´a umoˇzn ˇuje z´ıskat ˇreˇsen´ı u ´lohy pomoc´ı koneˇcn´eho poˇctu operac´ı. Tyto metody hraj´ı v numerick´e line´arn´ı algebˇre d˚ uleˇzitou roli. Klasick´ ymi pˇr´ıklady jsou Gaussova eliminace nebo metoda LU-rozkladu. Druhou skupinu tvoˇr´ı gradientn´ı metody jako napˇr´ıklad metoda sdruˇzen´ ych gradient˚ u nebo metoda stˇr´ıdav´ ych smˇer˚ u. Pˇr´ım´e metody jsou vˇsak pˇri ˇreˇsen´ı nˇekter´ ych u ´loh vˇetˇsinou m´alo efektivn´ı. Velmi d˚ uleˇzit´ ym prostˇredkem ˇreˇsen´ı soustav line´arn´ıch rovnic jsou tzv. iteraˇcn´ı metody. Ve sv´e diplomov´e pr´aci se tedy hlavnˇe zab´ yv´am iteraˇcn´ımi metodami pro ˇreˇsen´ı soustav s r˚ uzn´ ymi speci´aln´ımi typy matic. V prvn´ı kapitole jsou definov´any z´akladn´ı pojmy, kter´e budou potˇreba v dalˇs´ım v´ ykladu. A to definice matice a r˚ uzn´e operace s nimi, d´ale souvislost graf˚ u a matic a nakonec n´aznak pˇr´ım´ ych metod pro hled´an´ı ˇreˇsen´ı soustav line´arn´ıch rovnic. Ve druh´e kapitole jsou pops´any r˚ uzn´e speci´aln´ı typy matic, jejich vlastnosti, hled´an´ı vlastn´ıch ˇc´ısel a vektor˚ u. Nakonec se zab´ yv´am konkr´etn´ımi iteraˇcn´ımi metodami pro r˚ uzn´e matice soustav. V posledn´ı kapitole je to konkr´etnˇe pˇr´ıpad semiiteraˇcn´ı metody pro singul´arn´ı syst´emy. V´ yklad je doplnˇen pˇr´ıklady, na kter´ ych jsou jasnˇe uk´az´any teoretick´e poznatky. V pˇr´ıloze pak najdeme programy vytvoˇren´e pro v´ ypoˇcetn´ı syst´em Matlab a k nˇekter´ ym pˇr´ıklad˚ um obr´azky z´ıskan´e pomoc´ı tˇechto program˚ u.
5
Kapitola 1 Z´ akladn´ı pojmy 1.1
Matice a operace s nimi
V prvn´ı kapitole jsou uvedeny z´akladn´ı pojmy, kter´e se t´ ykaj´ı matic a nejd˚ uleˇzitˇejˇs´ı operace s maticemi. Definice 1.1.1. Matice typu (m, n) je soustava mn ˇc´ısel uspoˇr´adan´ych do m ˇr´adk˚ u a n sloupc˚ u a11 . . . a1n . .. ... A = .. . . am1 . . . amn
ˇ ıslo aik je prvek matice na m´ıstˇe (i, k), tj. v i -t´em Oznaˇcujeme ji A = (aik ). C´ ˇr´adku a k -t´em sloupci. Jsou-li prvky re´aln´e, ˇr´ık´ame, ˇze matice je re´aln´a , jsou-li komplexn´ı, matice je komplexn´ı. n-rozmˇern´y vektor je potom matice typu (n, 1). Je-li m = n, pak mluv´ıme o ˇctvercov´e matici n-t´eho ˇr´adu. Jestliˇze a = (a1 , . . . , an )T a b = (b1 , . . . , bn )T , pak jejich skal´ arn´ı souˇcin je T
ha, bi = a b =
n X
aj b j .
j=1
D´ale uvedeme r˚ uzn´e speci´aln´ı typy a vlastnosti matic. ˇ Rekneme, ˇze ˇctvercov´a matice A je regul´arn´ı, jestliˇze k n´ı existuje matice B tak, ˇze AB = BA = I. Jestliˇze takov´a matice neexistuje, potom A se naz´ yv´a singul´arn´ı. −1 Matice B je pak inverzn´ı matice k A a oznaˇcujeme ji A . ˇ Ctvercov´ a matice A se naz´ yv´a ryze ˇr´adkovˇe diagon´ alnˇe dominantn´ı (resp. matice s pˇrevl´ adaj´ıc´ı diagon´alou), jestliˇze |aii | >
n X j=1 j6=i
|aij |,
i = 1, . . . , n.
ˇ Rekneme, ˇze ˇctvercov´a matice A je podobn´a matici B, existuje-li regul´arn´ı matice 6
ˇ T tak, ˇze A = T BT −1 . Rekneme, ˇze matice A je kongruentn´ı s matic´ı B, jestliˇze existuje permutaˇcn´ı matice P takov´a, ˇze plat´ı B = P AP T . Pˇritom permutaˇcn´ı matice je takov´a, kter´a m´a v kaˇzd´em ˇr´adku a kaˇzd´em sloupci jedin´ y nenulov´ y prvek, rovn´ y jedn´e. Definice 1.1.2. Necht’ A je ˇctvercov´a matice. Nenulov´y vektor x se naz´yv´a vlastn´ı vektor matice A, plat´ı-li Ax = λx pro nˇejak´e ˇc´ıslo λ. Toto λ se naz´yv´a vlastn´ı ˇc´ıslo matice A odpov´ıdaj´ıc´ı vlastn´ımu vektoru x. Mnoˇzina vˇsech vlastn´ıch ˇc´ısel matice se naz´yv´a spektrum matice, ozn. σ(A). Spektr´ aln´ım polomˇerem pak naz´yv´ame ˇc´ıslo ρ(A) = max{|λ|, λ ∈ σ(A)}. Oznaˇcen´ı: Necht’ A = (aik ), matice typu (n, n). Vynech´ame-li v t´eto matici i -t´ y ˇr´adek a k -t´ y sloupec, dostaneme matici typu (n-1, n-1), kterou oznaˇc´ıme Aik . Definice 1.1.3. Determinant ˇctvercov´e matice A = (aik ) n-t´eho ˇr´adu je ˇc´ıslo X detA = sgn(P )a1k1 . . . ankn , P =(k1 ,...,kn )
kde P jsou vˇsechny permutace index˚ u 1, . . . , n a sgn(P ) je znam´enko permutace P. Pˇripomeˇ nme nyn´ı z´akladn´ı poznatek z line´arn´ı algebry. Soustavou n line´arn´ı rovnic o n nezn´am´ ych rozum´ıme soustavu Ax = b, kde A, matice soustavy, je . . . a1n .. ... . . . . amn
a11 .. A= .
am1
a b = (b1 , . . . , bn )T .
Vˇ eta 1.1.1 (Frobenius). Syst´em Ax = b je ˇreˇsiteln´y pr´avˇe tehdy, kdyˇz hodnost matice A je rovna hodnosti rozˇs´ıˇren´e matice syst´emu (A|b).
1.2
Grafy a matice
S maticemi pomˇernˇe u ´zce souvis´ı pojem graf. Nˇekter´e vlastnosti matice m˚ uˇzeme dokonce z grafu matice odvozovat. Proto zde uvedeme nˇekolik z´akladn´ı definic a vˇet t´ ykaj´ıc´ıch se graf˚ u a matic. ~ Definice 1.2.1. Orientovan´y graf G(V, H) je uspoˇr´adan´a dvojice koneˇcn´ych mnoˇzin V, H, kde H je tvoˇrena nˇekter´ymi uspoˇr´adan´ymi dvojicemi prvk˚ u z V: H ⊂ V × V . 7
Prvky mnoˇziny V se naz´ yvaj´ı uzly (vrcholy, body), prvky mnoˇziny H hrany. Jestliˇze hrana zaˇc´ın´a i konˇc´ı ve stejn´em uzlu, mluv´ıme o smyˇcce. Posloupnost uzl˚ u (u1 , . . . , uk ) takov´a, ˇze (u1 , u2 ), (u2 , u3 ), . . ., (uk−1 , uk ) jsou hrany, se naz´ yv´a spojen´ı; neopakuj´ı-li se uzly, oznaˇc´ıme ji jako dr´aha. Pˇritom poˇcet hran ve spojen´ı nazveme d´elkou spojen´ı. ~1 = (V1 , H1 ) je podgrafem grafu G ~2 = (V2 , H2 ), jestliˇze je V1 ⊂ V2 a H1 ⊂ H2 . G ~ je spojen´ı uzl˚ Definice 1.2.2. Cyklus v grafu G u (u1 , . . . , um , u1 ). D´elkou cyklu je pak ˇc´ıslo m. Pˇritom pˇredpokl´ad´ame, ˇze vˇsechny uzly jsou navz´ajem r˚ uzn´e. Orientovan´ y graf nazveme acyklick´y, jestliˇze neobsahuje ˇz´adn´ y cyklus. ~ = (V, H) je orientovan´y Vˇ eta 1.2.1 (Z´akladn´ı vˇeta o acyklick´ ych grafech). Necht’ G graf, |V | = n. Potom n´ asleduj´ıc´ı vlastnosti jsou ekvivalentn´ı: ~ je acyklick´y. 1. Graf G ~ ′ grafu G ~ m´a vlastnost, ˇze v nˇem existuje uzel, 2. Kaˇzd´y nepr´ azdn´y podgraf G ~ ′ nejde ˇza´dn´a hrana. do nˇehoˇz v G ~ 3. Existuje oˇc´ıslov´an´ı mnoˇziny uzl˚ u V ˇc´ısly 1, 2, . . . , n takov´e, ˇze kaˇzd´ a hrana v G jde z uzlu s menˇs´ım ˇc´ıslem do uzlu s ˇc´ıslem vˇetˇs´ım. ~ nem´a ˇza´dnou smyˇcku a pro kaˇzdou dvojici r˚ 4. G uzn´ych uzl˚ u u, v bud’ neexistuje ~ dr´aha z u do v nebo neexistuje v G ~ dr´aha z v do u. vG D˚ u k a z. viz [3] ~ se naz´ ~ dr´aha do kaˇzGraf G yv´a silnˇe souvisl´y, jestliˇze existuje z kaˇzd´eho uzlu v G ~ d´eho jin´eho uzlu grafu G. A = (aik ) je ˇctvercov´a matice ˇr´adu n. Oznaˇcme N mnoˇzinu index˚ u 1, 2, . . . , n. ~ Matici A pˇriˇrad´ıme orientovan´ y graf G(A) o n uzlech takto: ~ G(A) = (N, H), kde H je mnoˇzina dvojic (i, k), i ∈ N , k ∈ N , pro kter´e aik 6= 0. Hranovˇe ohodnocen´y graf je graf, v nˇemˇz je kaˇzd´e hranˇe pˇriˇrazena jej´ı hodnota. ~ = Orientovan´emu grafu tedy m˚ uˇzeme naopak pˇriˇradit tzv. uzlovou matici U (G) (uik ): uik = 1, je-li v grafu hrana z uzlu i do uzlu k a uik = 0, nen´ı-li v grafu takov´a ˇ ad uzlov´e matice je roven poˇctu uzl˚ hrana. R´ u grafu. ˇ Definice 1.2.3. Ctvercov´ a matice A se naz´yv´a reducibiln´ı (rozloˇziteln´a), je-li tvaru ! A1 B , 0 A2 kde A1 , A2 jsou ˇctvercov´e matice ˇr´adu alespoˇ n 1, anebo lze-li A pˇrev´est na tento tvar ˇ permutac´ı ˇr´adk˚ u a sloupc˚ u. Ctvercov´ a matice se naz´yv´a ireducibiln´ı (nerozloˇziteln´a), nen´ı-li reducibiln´ı. 8
ˇ Vˇ eta 1.2.2. Ctvercov´ a matice je ireducibiln´ı, pr´avˇe kdyˇz jej´ı orientovan´y graf je silnˇe souvisl´y. Kaˇzdou ˇctvercovou matici lze permutac´ı ˇr´adk˚ u a sloupc˚ u (odpov´ıdaj´ıc´ı permutaˇcn´ı matici P) pˇrev´est na horn´ı blokovˇe troj´ uheln´ıkov´y tvar, jehoˇz diagon´ aln´ı bloky jsou uˇz ireducibiln´ı: A11 A12 . . . A1r 0 A22 . . . A2r T P AP = .. . . .. .. . . . . . 0 0 . . . Arr D˚ u k a z. viz. [3]
Tento tvar matice se naz´ yv´a norm´ aln´ı tvar. Druh´a ˇc´ast z t´eto vˇety m´a rozs´ahl´e aplikace v numerick´e matematice, zejm´ena pˇri ˇreˇsen´ı rovnic a v´ ypoˇctu vlastn´ıch ˇc´ısel matic. Definice 1.2.4. Koneˇcn´y neorientovan´y graf (kr´atce graf ) G = (V, H) je uspoˇr´adan´a dvojice koneˇcn´ych mnoˇzin (V, H). V je mnoˇzina uzl˚ u, H mnoˇzina nˇekter´ ych neuspoˇr´adan´ ych dvojic prvk˚ u z V , tzv. neorientovan´ych hran. Sled v grafu G je posloupnost uzl˚ u u1 , . . . , us , kde kaˇzd´e dva po sobˇe jdouc´ı uzly uk , uk+1 , pro k = 1, . . . , s − 1 jsou spojeny hranou v G. Cesta v grafu G je sled, v nˇemˇz se ˇz´adn´e dva uzly neopakuj´ı. Souvisl´y graf je graf, v nˇemˇz mezi kaˇzd´ ymi dvˇema r˚ uzn´ ymi uzly existuje cesta.
1.3
Pˇ r´ım´ e metody ˇ reˇ sen´ı soustav line´ arn´ıch rovnic
V t´eto diplomov´e pr´aci se zab´ yv´am pˇrev´aˇznˇe iteraˇcn´ımi metodami pro ˇreˇsen´ı soustav line´arn´ıch rovnic. Proto bych zde r´ada kr´atce pˇripomnˇela i z´akladn´ı pˇr´ım´e metody, kter´e jsou tak´e d˚ uleˇzit´e pro ˇreˇsen´ı tˇechto soustav.
1.3.1
Gaussova eliminaˇ cn´ı metoda
Pˇredpokl´adejme, ˇze syst´em line´arn´ıch rovnic je zaps´an v maticov´em tvaru: Ax = b. Pˇredpokl´ad´ame, ˇze matice koeficient˚ u na lev´e stranˇe nen´ı singul´arn´ı a ˇze koeficient a11 nen´ı nulov´ y. Prvn´ı ˇr´adek vztahu dˇel´ıme koeficientem a11 . Pak postupnˇe pro vˇsechny hodnoty indexu i = 2, . . . , n od i-t´eho ˇr´adku odeˇc´ıt´ame nov´ y prvn´ı
9
ˇr´adek, vyn´asoben´ y koeficientem ai1 . Dost´av´ame x1 1 a112 . . . a11n 1 1 0 a22 . . . a2n x2 . . . . . . . .. .. . . . . 0 a1n2 . . . a1nn
xn
a
=
b11 b12 .. . b1n
,
ij kde a1ij = a11 , pro j = 1, . . . , n. V pˇr´ıpadˇe, ˇze by koeficient a11 byl nulov´ y, je vˇzdy moˇzn´e vymˇenit poˇrad´ı rovnic tak, aby nov´ y koeficient nulov´ y nebyl. Stejn´ ym zp˚ usobem pokraˇcujeme d´ale. Nakonec (po proveden´em n-t´em kroku) m´a v´ ysledn´ y syst´em line´arn´ıch rovnic tvar: 1 b1 x1 1 a112 a113 . . . a11n 2 2 2 0 1 a23 . . . a2n x2 b2 3 3 = b3 . 0 0 x 1 . . . a 3 3n .. . . .. .. .. .. .. . . . . . . . bnn xn 0 0 0 ... 1
Vˇsechny prvky matice pod hlavn´ı diagon´alou jsou nulov´e, diagon´aln´ı prvky jsou jedniˇcky. Horn´ı index vˇzdy ukazuje pˇr´ısluˇsn´ y krok, ve kter´em byl koeficient z´ısk´an. ˇ Reˇsen´ı nyn´ı vypoˇc´ıt´ame pomoc´ı vztah˚ u: xi =
bii
−
n X
aiij xj ,
j=i+1
i = n, n − 1, . . . , 1.
ˇ ste pomoc´ı Gaussovy eliminaˇcn´ı metody soustavu line´ Pˇ r´ıklad 1.3.1. Reˇ arn´ıch rovnic. 7, 9x1 + 5, 6x2 + 5, 7x3 − 7, 2x4 = 6, 68
8, 5x1 − 4, 8x2 + 0, 8x3 + 3, 5x4 = 9, 95
4, 3x1 + 4, 2x2 − 3, 2x3 + 9, 3x4 = 8, 60
3, 2x1 − 1, 4x2 − 8, 9x3 + 3, 3x4 = 1, 00. ˇ sen´ı: Soustavu rovnic si pˇrep´ıˇseme do maticov´eho tvaru a postupnˇe upravuReˇ jeme, jak bylo uvedeno v pˇredchoz´ı teoretick´e ˇc´asti.
7, 9 5, 6 5, 7 −7, 2 8, 5 −4, 8 0, 8 3, 5 4, 3 4, 2 −3, 2 9, 3 3, 2 −1, 4 −8, 9 3, 3
x1 x2 x3 x4
=
x1 1 0, 70886 0, 72152 −0, 91139 0 −10, 82531 −5, 33292 11, 24682 x2 0 1, 15190 −6, 30254 13, 21898 x3 x4 0 −3, 66835 −11, 20886 6, 21645 10
6, 68 9, 95 8, 60 1, 00
=
. 0, 84557 2, 76266 4, 96405 −1, 70582
x1 0, 84557 1 0, 70886 0, 72152 −0, 91139 0 1 0, 49263 −1, 03894 x2 −0, 25520 = 0 0 −6, 8700 14, 41573 x3 5, 25801 −2, 64198 0 0 −9, 40172 2, 40525 x4 0, 84557 x1 1 0, 70886 0, 72152 −0, 91139 0 1 0, 49263 −1, 03894 x2 −0, 25520 = 0 0 1 −2, 09836 x3 −0, 76536 −9, 83768 x4 0 0 0 −17, 32294 1 0, 70886 0, 72152 −0, 91139 x1 0, 84557 0 1 0, 49263 −1, 03894 x2 −0, 25520 = 0 0 1 −2, 09836 x3 −0, 76536
0
0
0
1
x4
0, 56790
Z t´eto soustavy jiˇz velice jednoduˇse z´ısk´ame ˇreˇsen´ı p˚ uvodn´ı soustavy line´arn´ıch rovnic: x1 = 0, 9671, x2 = 0, 1248, x3 = 0, 4263, x4 = 0, 5679. V´ ysledek m˚ uˇzeme zkontrolovat dosazen´ım do soustavy rovnic. ˇ sen´ı nedourˇcen´ Reˇ ych soustav, neboli soustav se singul´arn´ı matic´ı koeficient˚ u, pomoc´ı Gaussovy eliminaˇcn´ı metody je analogick´e jako pˇri soustavˇe line´arn´ıch rovnic s regul´arn´ı matic´ı koeficient˚ u. Pˇ r´ıklad 1.3.2. Najdˇete ˇreˇsen´ı soustavy rovnic: x1 + x2 − 2x3 + x4 = 1
x1 − 3x2 + x3 + x4 = 0
4x1 − x2 − x3 − x4 = 1
4x1 + 3x2 − 4x3 − x4 = 2. ˇ sen´ı: Soustavu rovnic si opˇet pˇrep´ıˇseme v maticov´em tvaru a postupnˇe upReˇ ravujeme. 1 x1 1 1 −2 1 1 −3 1 1 x2 0 = 4 −1 −1 −1 x3 1 2 x4 4 3 −4 −1 1 1 −2 1 x1 1 0 −4 3 0 x2 −1 = 0 −5 7 −5 x3 −3 0 −1 4 −5 x4 −2 1 1 −2 1 x1 1 0 1 −3 0 x 1 2 4 4 = 7 0 0 13 −4 x −5 3 4 −5 x4 − 74 0 0 13 4 11
1 0 0 0
1 −2 1 1 − 43 0 0 1 − 20 13 0 0 0
x1 x2 x3 x4
1
1 4 = 7 − 13 0
Hodnost matice koeficient˚ u je 3. Proto poloˇz´ıme x4 = u, kde u je parametr. Obecn´e 1 2 ˇreˇsen´ı soustavy je tedy jednoparametrick´e a m´a tvar: x1 = 13 + 12 u, x2 = − 13 + 15 u, 13 13 7 20 x3 = − 13 + 13 u, x4 = u.
1.3.2
Metoda LU - rozkladu
Z´akladem t´eto metody ˇreˇsen´ı syst´em˚ u line´arn´ıch rovnic je rozklad ˇctvercov´e matice na souˇcin doln´ı a horn´ı troj´ uheln´ıkov´e matice. Libovolnou ˇctvercovou matici m˚ uˇzeme takto jednoznaˇcnˇe rozloˇzit, jestliˇze jsou d´any nenulov´e diagon´aln´ı prvky jedn´e z hledan´ ych troj´ uheln´ıkov´ ych matic a nejsou-li z´aroveˇ n hlavn´ı subdeterminanty rozkl´adan´e matice nulov´e. Napˇr´ıklad pro n = 4 tedy plat´ı A = LU , u11 u12 u13 u14 l11 0 0 0 a11 a12 a13 a14 a 21 a22 a23 a24 l21 l22 0 0 0 u22 u23 u24 . = a31 a32 a33 a34 l31 l32 l33 0 0 0 u33 u34 a41 a42 a43 a44
l41 l42 l43 l44
0
0
0
u44
Naˇs´ım u ´kolem je tedy urˇcit prvky lij a uij , pro dan´a i a j, a ty urˇc´ıme ze vztah˚ u, kter´e dostaneme rozn´asoben´ım tˇechto matic. T´ım dostaneme 16 rovnic pro 20 nezn´am´ ych a ˇctyˇri nezn´am´e veliˇciny si tedy m˚ uˇzeme libovolnˇe vybrat. Nejˇcastˇeji se pouˇz´ıv´a v´ ybˇer l11 = l22 = l33 = l44 = 1 nebo u11 = u22 = u33 = u44 = 1. Pˇri ˇreˇsen´ı soustavy line´arn´ıch rovnic Ax = b metodou LU-rozkladu vyj´adˇr´ıme matici A jako souˇcin horn´ı troj´ uheln´ıkov´e matice U a doln´ı troj´ uheln´ıkov´e matice L. Pak tedy m´ame LU x = b ˇ sen´ı soustav a oznaˇc´ıme-li U x = y, m˚ uˇzeme ˇreˇsit dvˇe soustavy U x = y a Ly = b. Reˇ ˇ sen´ı soustavy U x = y najdeme snadno, protoˇze L i U jsou troj´ uheln´ıkov´e matice. Reˇ je jiˇz ˇreˇsen´ı p˚ uvodn´ı soustavy line´arn´ıch rovnic. Pˇ r´ıklad 1.3.3. Metodou LU-rozkladu ˇreˇste n´ asleduj´ıc´ı soustavu line´ arn´ıch rovnic. 0, 500x1 + 1, 000x3 + 1, 023x4 = 4, 725 1, 500x1 + 1, 000x2 + 3, 702x4 = 3, 402 1, 273x1 − 2, 752x2 + 3, 208x3 − 1, 305x4 = 2, 709
2, 000x1 + 1, 000x2 + 1, 000x3 + 4, 007x4 = 1, 231. 12
ˇ sen´ı: Nejdˇr´ıve si soustavu pˇrep´ıˇseme do Reˇ 0, 5 0 1 1, 023 1, 5 1 0 3, 702 1, 273 −2, 752 3, 208 −1, 305 2 1 1 4, 007
maticov´eho tvaru x1 4, 725 x 3, 402 2 = x3 2, 709
1, 231
x4
.
Najdeme matice L a U (pˇriˇcemˇz vol´ıme lii = 1, i = 1, . . . , n) a pak ˇreˇs´ıme soustavy: 1 0 0 0 y1 4, 725 3 1 0 0 y2 3, 402 = 2, 546 −2, 752 1 0 y3 2, 709 4 1 0 1 y4 1, 231
a
0, 5 0 0 0
0 1 1, 023 1 −3 0, 633 0 −7, 594 −2, 168 0 0 −0, 718
x1 x2 x3 x4
=
4, 725 −10, 733 −38, 968 −6, 896
.
ˇ sen´ı prvn´ı soustavy je y = (4, 725, −10, 773, −38, 968, −6, 896)T a po dosazen´ı Reˇ tohoto vektoru do druh´e soustavy dostaneme jiˇz ˇreˇsen´ı p˚ uvodn´ıho syst´emu rovnic x1 = −14, 980, x2 = −9, 682, x3 = 2, 390 a x4 = 9, 604. Speci´aln´ım pˇr´ıpadem metody LU-rozkladu je Cholesk´eho metoda (metoda odmocnin). Tato metoda se pouˇz´ıv´a pro syst´emy line´arn´ıch rovnic tvaru Ax = b, kde matice A je symetrick´a. Symetrickou matici lze totiˇz vyj´adˇrit jako souˇcin dvou navz´ajem transponovan´ ych troj´ uhlen´ıkov´ ych matic a rovnice Ax = b pak pˇrejde na tvar TTTx = b a ˇreˇsen´ı pˇrejde na ˇreˇsen´ı dvou soustav T T y = b a T x = y.
1.3.3
Metoda nejvˇ etˇ s´ıho sp´ adu
Vˇ eta 1.3.1. Jestliˇze je A symetrick´a a pozitivnˇe definitn´ı matice typu (n, n) a x a b jsou vektory z Rn , potom kvadratick´ a forma F (x) = hAx, xi − 2hb, xi dosahuje minim´aln´ı hodnoty pro v´ybˇer x = x∗ pr´avˇe tehdy, kdyˇz je x∗ ˇreˇsen´ım soustavy Ax = b. D˚ u k a z. Mˇejme ˇreˇsen´ı x∗ soustavy Ax∗ = b a ∆x 6= o. Pak F (x∗ + ∆x) − F (x∗ ) =
= hA(x∗ + ∆x), x∗ + ∆xi − 2hb, x∗ + ∆xi − hAx∗ , x∗ i − 2hb, x∗ i = 13
= hAx∗ , x∗ i + hAx∗ , ∆xi + hA∆x, x∗ i +hA∆x, ∆xi− | {z } stejn´ e v´ yrazy
∗
−2hb, x i − 2hb, ∆xi − hAx∗ , x∗ i + 2hb, x∗ i = ∗ = 2hAx |{z}, ∆xi − 2hb, ∆xi + hA∆x, ∆xi = =b
= 2hb, ∆xi − 2hb, ∆xi + hA∆x, ∆xi = = hA∆x, ∆xi > 0,
nebot’ A je pozitivnˇe definitn´ı. Odtud jiˇz plyne F (x∗ + ∆x) > F (x∗ ). Naopak mˇejme F (ˆ x) = min F (x) a chceme dok´azat, ˇze Aˆ x = b. Necht’ t ∈ R, v ∈ Rn . Jestliˇze definujeme funkci F (ˆ x + tv), pak minimum funkce F najdeme pro t = 0. Tedy F (ˆ x + tv) = hA(ˆ x + tv), xˆ + tvi − 2hb, xˆ + tvi = = hAˆ x, xˆi + 2thAˆ x, vi + t2 hAv, vi − 2hb, xˆi − 2thb, vi. Nyn´ı zderivujeme podle promˇenn´e t a poloˇz´ıme rovno 0. dF (ˆ x + tv) =0 dt 2hAˆ x, vi + 2thAv, vi − 2hb, vi = 0. Protoˇze v´ıme, ˇze v bodˇe t = 0 nast´av´a minimum F , dostaneme hAˆ x, vi = hb, vi,
∀v ∈ Rn
a tedy Aˆ x = b.
ˇ sen´ı u Reˇ ´lohy Ax = b je tedy pro symetrickou a pozitivnˇe definitn´ı matici A moˇzn´e nahradit minimalizac´ı F (x). Zvol´ıme poˇc´ateˇcn´ı aproximaci minima x0 a vybereme vhodn´ y smˇer v 0 . Vektor x1 nyn´ı urˇc´ıme podle vzorce xk+1 = xk + δk v k , kde δk je konstanta. v k a δk vol´ıme tak, aby platilo F (xk+1 ) < F (xk ). Jestliˇze posloupnost x0 , x1 , x2 , . . . konverguje, mus´ı b´ yt jej´ı limita ˇreˇsen´ım syst´emu Ax = b. Metoda nejvˇetˇs´ıho sp´adu, neboli gradientn´ı metoda, je zaloˇzena na v´ ybˇeru vektoru v v kaˇzd´em kroku tak, aby leˇzel ve smˇeru nejvˇetˇs´ı zmˇeny funkce F , tj. ve smˇeru gradientu t´eto funkce: grad(F ) = (∂Fx1 , . . . , ∂Fxn )T . Kvadratickou formu F a jej´ı parci´aln´ı derivaci m˚ uˇzeme zapsat ve tvaru: F =
n X
i,j=1
aij xi xj − 2 14
n X i=1
bi xi ,
∂Fxi = 2
n X i=1
Oznaˇc´ıme-li reziduum ri =
− 12 ∂Fxi ,
aij xj − 2bi .
tj. ri = bi −
Pn
i=1
aij xj , potom plat´ı
grad(F ) = −2r,
kde r = (r1 , . . . , rn )T . V t´eto metodˇe zvol´ıme v k jako reziduum rk , kde rik
= bi −
n X
aij xkj , i = 1, . . . , n.
i=1
Koeficient λk vybereme nyn´ı tak, aby hodnota F (xk+1 ) = F (xk +λk rk ) byla minim´aln´ı. F (xk + λk rk ) = hA(xk + λk rk ), xk + λk rk i − 2hb, xk + λk rk i = = hAxk , xk i + 2λk hAxk , rk i + λ2k hArk , rk i − 2hb, xk i − 2λk hb, rk i = |rk = b − Axk a tedy Axk = b − rk | = hb, xk i − hrk , xk i + 2λk hb, rk i − 2λk hrk , rk i+ +λ2k hArk , rk i − 2hb, xk i − 2λk hb, rk i = = −hb, xk i − hrk , xk i − 2λk hrk , rk i + λ2k hArk , rk i. Nyn´ı zderivujeme podle λi a poloˇz´ıme rovno 0. T´ım dostaneme minimum. −2hri , ri i + 2λi hAri , ri i = 0 λi =
hri , ri i . hAri , ri i
Dost´av´ame tedy iteraˇcn´ı vzorec metody nejvˇetˇs´ıho sp´adu: xk+1 = xk +
hrk , rk i k r , hArk , rk i
kde k = 0, 1, . . .. Pˇ r´ıklad 1.3.4. Metodou nejvˇetˇs´ıho sp´ adu najdˇete ˇreˇsen´ı syst´emu rovnic 4x1 − x2 = 2
−x1 + 4x2 − x3 = 6
−x2 + 4x3 = 2.
ˇ sen´ı: Nejdˇr´ıve si zvol´ıme poˇc´ateˇcn´ı aproximaci x0 = (0, 0, 0)T . Rezidua budou Reˇ m´ıt tvar r1 = 2 − 4x1 + x2 ,
r2 = 6 + x1 − 4x2 + x3 ,
r3 = 2 + x2 − 4x3 . 15
Potom r0 = (2, 6, 2)T a urˇc´ıme koeficient λ0 hr0 , r0 i . λ0 = = 0, 3438. hAr0 , r0 i Prvn´ı iterace tedy bude x1 = x0 + λ0 r0 = (0, 6876, 2, 0628, 0, 6876)T . . Stejn´ ym postupem dostaneme r1 = (1, 3125, −0, 8750, 1, 3125)T , λ1 = 0, 1964 a tedy x2 = x1 + λ1 r1 = (0, 9453, 1, 8906, 0, 9453)T , . r2 = (0, 1094, 0, 3281, 0, 1094)T , λ2 = 03438 a tedy x3 = x2 + λ2 r2 = (0, 9829, 2, 0034, 0, 9829)T . Pro srovn´an´ı, pˇresn´e ˇreˇsen´ı soustavy je x∗ = (1, 2, 1)T .
16
Kapitola 2 Speci´ aln´ı matice 2.1
Nez´ aporn´ e matice
Tato kapitola se zab´ yv´a ˇctvercov´ ymi nez´aporn´ ymi maticemi, tedy takov´ ymi, kter´e maj´ı vˇsechny prvky nez´aporn´e. Nejdˇr´ıve zavedeme tato oznaˇcen´ı: A ≥ B, A > B,
jestliˇze aij ≥ bij pro vˇsechna i, j, jestliˇze aij > bij pro vˇsechna i, j.
Nyn´ı m˚ uˇzeme uv´est n´asleduj´ıc´ı definice a nˇekter´e vˇety o nez´aporn´ ych matic´ıch. ˇ Definice 2.1.1. Rekneme, ˇze matice A je nez´ aporn´a, kdyˇz A ≥ 0, tedy vˇsechny jej´ı prvky jsou nez´ aporn´e. Matice A je kladn´a, kdyˇz A > 0, tedy vˇsechny jej´ı prvky jsou kladn´e. Nez´aporn´e matice maj´ı nˇekter´e zˇrejm´e vlastnosti: Vˇ eta 2.1.1. Mˇejme dvˇe nez´ aporn´e matice A, B. Pak plat´ı: 1. Jsou-li A a B t´ehoˇz typu, je A + B opˇet nez´ aporn´a matice. 2. Lze-li A a B n´ asobit, je jejich souˇcin AB nez´ aporn´a matice. 3. Souˇcin AB kladn´e matice A a nez´ aporn´e nenulov´e matice B je nez´aporn´a nenulov´ a matice. D˚ u k a z. Zˇrejm´ y.
17
Z toho vypl´ yv´a, ˇze podobn´a tvrzen´ı plat´ı i pro n´asoben´ı matice a vektoru (vektor je matice typu (n, 1)). D´ale je uvedena jeˇstˇe jin´a definici reducibiln´ı a ireducibiln´ı matice, neˇz jak´a je uvedena v kapitole Grafy a matice. Definice 2.1.2. Matice A ≥ 0 je reducibiln´ı, jestliˇze ∃i, j
∀k ∈ N : (Ak )ij = 0.
Matice A ≥ 0 je ireducibiln´ı, jestliˇze ∀i, j
∃k ∈ N : (Ak )ij > 0.
Pozn´ amka 2.1.1. (Ak )ij oznaˇcuje prvek na pozici (i, j) v k-t´e mocninˇe matice A, tj. Ak . Zaved’me nyn´ı pojem struktury nenulov´ ych prvk˚ u matice, kter´a se nejˇcastˇeji vyˇsetˇruje pomoc´ı tzv. booleovsk´ ych matic. Definice 2.1.3. Dvˇe matice A = (aik ) a B = (bik ) maj´ı stejnou strukturu nenulov´ych prvk˚ u, jestliˇze pro vˇsechna i, k plat´ı: aik 6= 0,
⇔
bik 6= 0.
Struktura nenulov´ ych prvk˚ u mocniny ˇctvercov´e nez´aporn´e matice A souvis´ı tak´e ~ se sledy v orientovan´em grafu G(A) t´eto matice. Tato souvislost je pops´ana v n´asleduj´ıc´ı vˇetˇe. Vˇ eta 2.1.2. Je-li A ˇctvercov´a nez´ aporn´a matice a k ∈ N, pak v mocninˇe Ak je na ~ m´ıstˇe (i, j) nenulov´y prvek, pr´avˇe kdyˇz existuje v orientovan´em grafu G(A) sled d´elky k z uzlu i do uzlu j. D˚ u k a z. Vˇetu dok´aˇzeme matematickou indukc´ı. Pro k = 1 vˇeta zˇrejmˇe plat´ı. Pˇredpokl´adejme tedy nyn´ı, ˇze vˇeta plat´ı pro k a dok´aˇzeme, ˇze plat´ı i pro k + 1. Oznaˇcme A = (apq ), Ak = B = (bpq ), Ak+1 = C = (cpq ). Protoˇze tedy C = AB, pak cij =
X
bip apj .
p
~ Necht’ v orientovan´em grafu G(A) existuje sled d´elky k + 1 z i do j: (i, p1 , . . . , pk , j). Podle indukˇcn´ıho pˇredpokladu je bipk > 0. Protoˇze tak´e apk j > 0 a vˇsechny sˇc´ıtance v pˇredchoz´ı sumˇe jsou nez´aporn´e, je i cij > 0, tj. v Ak+1 na m´ıstˇe (i, j) nenulov´ y (tedy kladn´ y) prvek. Pak je v sumˇe nˇekter´ y z ˇclen˚ u kladn´ y, napˇr. bit atj > 0. To znamen´a, ˇze bit > 0 a atj > 0 a podle indukˇcn´ıho pˇredpokladu tedy existuje v grafu ~ ~ G(A) sled d´elky k z i do t. A protoˇze (i, t) je hrana v G(A), existuje i sled d´elky k + 1 z i do j. 18
Vˇ eta 2.1.3. Je-li A nez´aporn´a nerozloˇziteln´a matice n-t´eho ˇr´adu a k0 , k1 , . . . , kn−1 kladn´a ˇc´ısla, pak matice k0 I + k1 A + k2 A2 + . . . + kn−1 An−1 je kladn´a. Speci´alnˇe je (I + A)n−1 > 0. P i D˚ u k a z. Protoˇze matice n−1 cet nez´aporn´ ych matic, staˇc´ı dok´azat, i=0 ki A je souˇ p ˇze pro libovoln´a pevn´a i a j je u nˇekter´e mocniny A , 0 ≤ p ≤ n − 1, na m´ıstˇe (i, j) kladn´ y prvek. To plat´ı pro i = j, nebot’ I m´a na tomto m´ıstˇe kladn´ y prvek. Je-li ~ i 6= j, pak ze siln´e souvislosti grafu G(A) plyne, ˇze existuje dr´aha z i do j. D´elka ~ d t´eto dr´ahy nepˇrev´ yˇs´ı n − 1, nebot’ v G(A) je n uzl˚ u a dr´aha obsahuje jen r˚ uzn´e d uzly. Podle pˇredchoz´ı vˇety je vˇsak prvek na m´ıstˇe (i, j) matice A , 1 ≤ d ≤ n − 1, kladn´ y. T´ım je prvn´ı tvrzen´ı dok´az´ano. Druh´e tvrzen´ı je d˚ usledkem prvn´ıho, protoˇze s pouˇzit´ım binomick´e vˇety dostaneme ! ! n − 1 n − 1 (I + A)n−1 = I + A + ... + Ak + . . . + An−1 . 1 k
Nyn´ı m˚ uˇzeme uv´est hlavn´ı vˇetu o nez´aporn´ ych matic´ıch. Vˇ eta 2.1.4 (Perronova-Frobeniova vˇeta o nez´aporn´ ych matic´ıch). Necht’ A je ˇctvercov´a nez´ aporn´a ireducibiln´ı matice n-t´eho ˇr´adu, n > 1. Pak spektr´ aln´ı polomˇer ρ(A) je kladn´e jednoduch´e vlastn´ı ˇc´ıslo matice A a tomuto vlastn´ımu ˇc´ıslu odpov´ıd´a kladn´y ˇadn´emu jin´emu vlastn´ımu ˇc´ıslu matice A uˇz neodpov´ıd´a nez´aporn´y vlastn´ı vektor. Z´ vlastn´ı vektor. K d˚ ukazu t´eto vˇety potˇrebujeme jeˇstˇe jedno lemma, kter´e je zde uvedeno bez d˚ ukazu. Lemma 2.1.1 (Perron). Je-li A > 0, pak ρ(A) je kladn´e vlastn´ı ˇc´ıslo matice A. Tomuto vlastn´ımu ˇc´ıslu odpov´ıd´a jedin´y line´ arnˇe nez´ avisl´y vlastn´ı vektor, kter´y pˇritom lze volit kladn´y. D˚ u k a z. (Perronovy-Frobeniovy vˇety): Nejdˇr´ıve oznaˇcme m(A) matici (|aik |) a nazvˇeme ji modul matice A a tedy m(x) = (|x1 |, . . . , |xn |)T je modul vektoru x. Nyn´ı mˇejme ireducibiln´ı matici A. Matice (I + A)n−1 je kladn´a a tedy i matice T (I + AT )n−1 = (I + A)n−1 je kladn´a. Podle Perronova lemmatu existuje vektor y > 0 tak, ˇze T (I + A)n−1 y = ρ ((I + A)n−1 )T y neboli
y T (I + A)n−1 = ρ (I + A)n−1 y T .
Necht’ d´ale λ je vlastn´ı ˇc´ıslo matice A takov´e, ˇze |λ| = ρ(A). 19
Necht’ x je nˇekter´ y vlastn´ı vektor odpov´ıdaj´ıc´ı vlastn´ımu ˇc´ıslu λ (tedy x 6= 0), tj. Ax = λx. Plat´ı pak |λ|m(x) ≤ Am(x) neboli ρ(A)m(x) ≤ Am(x). Obdobnˇe ρ2 (A)m(x) ≤ ρ(A)Am(x) = Aρ(A)m(x) ≤ A2 m(x), a obecnˇe ρk (A)m(x) ≤ Ak m(x), N´asob´ıme-li k-tou z nerovnost´ı ˇc´ıslem
n−1 k
k = 1, 2, . . . . !
a seˇcteme pro k = 1, . . . , n − 1,
spolu s rovnost´ı m(x) = Im(x) dostaneme (1 + ρ(A))n−1 m(x) ≤ (I + A)n−1 m(x). N´asobme zleva kladn´ ym vektorem y T . M´ame pak (1 + ρ(A))n−1 (y T m(x)) ≤ y T (I + A)n−1 m(x). Prav´a strana je rovna ρ((I + A)n−1 )(y T m(x)). Protoˇze y T m(x) je kladn´e ˇc´ıslo, dost´av´ame (1 + ρ(A))n−1 ≤ ρ((I + A)n−1 ). Matice (I + A)n−1 m´a vlastn´ı ˇc´ısla tvaru (1 + α)n−1 , kde α prob´ıh´a vlastn´ı ˇc´ısla matice A. To znamen´a, ˇze existuje vlastn´ı ˇc´ıslo µ matice A tak, ˇze |(1 + µ)n−1 | = ρ((I + A)n−1 ). Pˇritom vˇsak |µ| ≤ ρ(A). Dosazen´ım do pˇredchoz´ı nerovnosti plyne (1 + ρ(A))n−1 ≤ |(1 + µ)n−1 | neboli 1 + ρ(A) ≤ |1 + µ| ≤ 1 + |µ| ≤ 1 + ρ(A). Protoˇze lev´a strana spl´ yv´a s pravou, plat´ı vˇsude rovnost. Speci´alnˇe odtud plyne, ˇze µ ≥ 0, 20
a tedy µ = ρ(A). Rovnost tak´e plat´ı ve vˇsech nerovnostech, kter´e jsme dˇr´ıve sˇc´ıtali. Speci´alnˇe pro k = 1 je Am(x) = ρ(A)m(x) neboli Am(x) = µm(x) a tak´e (I + A)n−1 m(x) = (1 + µ)n−1 m(x) = ρ((I + A)n−1 )m(x). Podle Perronova lemmatu je m(x) > 0. K vlastn´ımu ˇc´ıslu µ existuje jedin´ y line´arnˇe nez´avisl´ y vlastn´ı vektor. Nav´ıc ρ(A) > 0, protoˇze A je nenulov´a matice (n > 1!). Zb´ yv´a dok´azat, ˇze ρ(A) je jednoduch´e vlastn´ı ˇc´ıslo matice A. To vypl´ yv´a z tzv. Schurovy vˇety [A ˇctvercov´a matice, λ jej´ı vlastn´ı ˇc´ıslo, pak λ je jednoduch´e, pr´avˇe kdyˇz jsou splnˇeny dvˇe podm´ınky: a) existuje jedin´ y line´arnˇe nez´avisl´ y vlastn´ı vektor matice A a tedy tak´e jedin´ y line´arnˇe nez´avisl´ y vlastn´ı vektor matice AT ; b) o vektorech u a v z a) plat´ı v T u 6= 0]. Vlastn´ımu ˇc´ıslu ρ(A) totiˇz odpov´ıd´a jedin´ y line´arnˇe nez´avisl´ y vlastn´ı vektor matice A, napˇr. u, kter´ y je pˇritom kladn´ y. Pˇr´ısluˇsn´ y vlastn´ı T vektor v matice A , kter´a je tak´e nerozloˇziteln´a, lze volit kladn´ y v > 0. Proto je T jejich souˇcin hv, ui = v u kladn´ y a z Schurovy vˇety plyne jednoduchost vlastn´ıho ˇc´ısla ρ(A). Nakonec uk´aˇzeme, ˇze ˇz´adn´emu jin´emu vlastn´ımu ˇc´ıslu uˇz neodpov´ıd´a nez´aporn´ y vlastn´ı vektor. Necht’ tedy Az = ξz, z ≥ 0 a ξ 6= ρ(A). Podle dok´azan´eho existuje k matici AT kladn´ y vlastn´ı vektor w > 0: AT w = ρ(A)w. Potom vˇsak wT Az = wT ξz = ξ(wT z), ale tak´e wT Az = ρ(A)(wT z). Celkem tedy (ρ(A) − ξ)(wT z) = 0, coˇz vzhledem k ρ(A) − ξ 6= 0 a wT z > 0 je spor. Z Perronovy-Frobeniovy vˇety tak´e plyne n´asleduj´ıc´ı. Vezmˇeme nejmenˇs´ı kruh se stˇredem v poˇc´atku, kter´ y obsahuje vˇsechna vlastn´ı ˇc´ısla nez´aporn´e ireducibiln´ı matice - tento kruh m´a polomˇer rovn´ y spektr´aln´ımu polomˇeru matice. Pak vlastn´ı ˇc´ısla t´eto matice jsou v komplexn´ı rovinˇe um´ıstˇena tak, ˇze na hranici uvaˇzovan´eho 21
kruhu je vˇzdy na kladn´e poloose vlastn´ı ˇc´ıslo. To vˇsak neznamen´a, ˇze na hranici tohoto kruhu uˇz nemohou b´ yt jin´a vlastn´ı ˇc´ısla. Plat´ı tak´e zaj´ımav´a vˇeta: Vˇ eta 2.1.5. Necht’ A je nez´ aporn´a ireducibiln´ı matice n-t´eho ˇr´adu. Necht’ h je pˇrirozen´e ˇc´ıslo. Pak tyto vlastnosti matice A a ˇc´ısla h jsou ekvivalentn´ı: 1. Existuje pr´avˇe h r˚ uzn´ych vlastn´ıch ˇc´ısel matice A, jejichˇz absolutn´ı hodnota je rovna ρ(A). 2. Existuje permutaˇcn´ı matice P tak, ˇze P AP T 0 A12 0 0 A23 0 . .. .. T . P AP = . . . 0 0 0 Ah1 0 0
m´a tvar 0 0 .. .
... ... ...
. . . Ah−1,h ... 0
se ˇctvercov´ymi diagon´aln´ımi bloky, a ˇza´dnou permutaˇcn´ı matic´ı nelze A pˇrev´est na obdobn´y tvar s v´ıce neˇz h blokov´ymi ˇr´adky.
~ 3. Nejvˇetˇs´ı spoleˇcn´y dˇelitel d´elek vˇsech cykl˚ u orientovan´eh grafu G(A) matice A je h. 4. Je-li (−1)n λn + kn1 λn1 + kn2 λn2 + . . . + kns λns charakteristick´y polynom matice A, kn1 6= 0, kn2 6= 0, . . . , kns 6= 0, n > n1 > . . . > ns ≥ 0, pak nejvˇetˇs´ı spoleˇcn´y dˇelitel ˇc´ısel n − n1 , n1 − n2 , . . . , ns−1 − ns je h. 5. O spektru σ(Z) matice Z plat´ı h = max{k ∈ N, σ(e
2πi k
A) = σ(A)}.
D˚ u k a z. Oznaˇc´ıme-li ht , t = 1, . . . , 5, ˇc´ıslo h v t-t´em tvrzen´ı, staˇc´ı dok´azat, ˇze h1 ≤ h2 , h2 ≤ h3 , h3 ≤ h4 , h4 ≤ h5 , h5 ≤ h1 . Cel´ y d˚ ukaz najdeme v [3]. D˚ usledek 2.1.1. Je-li A nez´ aporn´a ireducibiln´ı matice a m´a-li A pr´avˇe h r˚ uzn´ych vlastn´ıch ˇc´ısel o absolutn´ı hodnotˇe ρ(A), pak tato vlastn´ı ˇc´ısla tvoˇr´ı v rovinˇe komplexn´ıch ˇc´ısel vrcholy pravideln´eho h-´ uheln´ıku o stˇredu v poˇc´atku, jehoˇz jeden vrchol je ρ(A). Vˇsechna tato vlastn´ı ˇc´ısla jsou jednoduch´a. Pˇ r´ıklad 2.1.1. h = 1 → v komplexn´ı rovinˇe je jen bod ρ(A), h = 2 → v komplexn´ı rovinˇe jsou body ρ(A) a −ρ(A). Pozn´ amka 2.1.2. Ireducibiln´ı nez´ aporn´a matice, pro kterou je h = 1, se naz´yv´a primitivn´ı. Je-li h > 1, naz´yv´a se imprimitivn´ı a h je index imprimitivity. Existuje jeˇstˇe jin´a definice primitivn´ı matice. Matice A ≥ 0 se naz´ yv´a primitivn´ı, jestliˇze ∃k ∈ N : Ak > 0. Matice A ≥ 0 se naz´ yv´a imprimitivn´ı, jestliˇze ∀k ∈ N ∃i, j : (Ak )ij = 0. 22
Vˇ eta 2.1.6. Je-li A nez´aporn´a ˇctvercov´a matice, pak ρ(A) je vlastn´ım ˇc´ıslem matice A a existuje nez´aporn´y vlastn´ı vektor matice A, odpov´ıdaj´ıc´ı tomuto vlastn´ımu ˇc´ıslu. D˚ u k a z. viz [3] Vlastn´ı ˇc´ıslo ρ(A) nez´aporn´e matice A se zpravidla naz´ yv´a Perronovo vlastn´ı ˇc´ıslo a odpov´ıdaj´ıc´ı nez´aporn´ y vlastn´ı vektor Perron˚ uv vlastn´ı vektor. Vˇ eta 2.1.7. Necht’ A je nez´ aporn´a matice n-t´eho ˇr´adu, n ≥ 2. Potom n´ asleduj´ıc´ı podm´ınky jsou ekvivalent´ı: 1. An−1 = 0, 2. existuje k ∈ N tak, ˇze Ak = 0, ~ 3. orientovan´y graf G(A) matice A je acyklick´y, 4. existuje permutaˇcn´ı matice P tak, ˇze P AP T je horn´ı troj´ uheln´ıkov´a matice s nulami na hlavn´ı diagon´ ale, 5. ρ(A) = 0. D˚ u k a z. Je tˇreba dok´azat implikace: 1. ⇒ 2., 2. ⇒ 3., 3. ⇒ 4., 4. ⇒ 5. a 5. ⇒ 1. 1. ⇒ 2. Zˇrejm´e. ~ 2. ⇒ 3. Obsahuje-li graf G(A) matice A cyklus d´elky s a je-li uzel j v tomto cyklu, je j-t´ y diagon´aln´ı prvek kaˇzd´e z matic As , A2s , A3s , . . . kladn´ y. Nem˚ uˇze tedy pro k ˇz´adn´e k platit A = 0. 3. ⇒ 4. Z teorie graf˚ u v´ıme, ˇze kaˇzd´a siln´a komponenta acyklick´eho grafu je jednouzlov´a; existuje tedy permutaˇcn´ı matice P takov´a, ˇze P AP T je horn´ı troj´ uhlen´ıkov´a ~ matice. Diagon´aln´ı prvky jsou t´emˇeˇr vˇsechny rovny nule, nebot’ by v G(A) existovala smyˇcka, coˇz je cyklus d´elky 1. 4. ⇒ 5. Protoˇze P AP T je horn´ı troj´ uheln´ıkov´a matice s nulov´ ymi diagon´aln´ımi T prvky, jsou vˇsechna vlastn´ı ˇc´ısla matice P AP nulov´a a tedy i vˇsechna vlastn´ı ˇc´ısla matice A jsou nulov´a. Proto je ρ(A) = 0. 5. ⇒ 1. Z 5. vypl´ yv´a, ˇze A m´a vˇsechna vlastn´ı ˇc´ısla nulov´a. Proto v Jordanovˇe tvaru A jsou vˇsechny Jordanovy bloky tvaru 0 1 0 ... 0 0 0 1 ... 0 . . . . .. .. .. . . ... . 0 0 0 ... 1 0 0 0 ... 0 Pˇritom je vˇzdy (n−1)-n´ı mocnina tohoto Jordanova bloku nulov´a. Proto An−1 = 0.
23
Vˇ eta 2.1.8. Necht’ A je nez´ aporn´a ˇctvercov´a matice, kter´ a m´a kladn´y vlastn´ı vektor. Potom je tento vektor Perron˚ uv vlastn´ı vektor, tj. vlastn´ı ˇc´ıslo odpov´ıdaj´ıc´ı tomuto vektoru je ρ(A). D˚ u k a z. Oznaˇcme u kladn´ y vlastn´ı vektor matice A, takˇze u > 0 a Au = cu. Pak existuje nez´aporn´ y vlastn´ı vektor v matice AT , odpov´ıdaj´ıc´ı ρ(AT ) = ρ(A). Plat´ı tedy v T Au = cv T u a tak´e v T Au = ρ(A)v T u. Protoˇze v T u > 0, plyne z porovn´an´ı obou rovnost´ı, ˇze c = ρ(A).
2.2
Stochastick´ e matice
Tzv. stochastick´e matice jsou d˚ uleˇzitou podtˇr´ıdou nez´aporn´ ych matic. Nejdˇr´ıve uvaˇzujme n stav˚ u s1 , s2 , . . . , sn nˇejak´eho procesu. Pˇredpokl´adejme, ˇze pravdˇepodobnost pohybu procesu ze stavu si do stavu sj je ˇcasovˇe nez´avisl´a, oznaˇcme ji aij . Takov´ y proces se naz´ yv´a koneˇcn´ y homogenn´ı markovsk´ y ˇretˇezec. Nyn´ı matice A = (aij ) je stochastick´a. ˇ Definice 2.2.1. Ctvercov´ a matice A ˇr´adu n se naz´yv´a (ˇr´adkovˇe) stochastick´a, jestliˇze splˇ nuje n X aij ≥ 0, aij = 1, i, j = 1, . . . , n. j=1
Nyn´ı oznaˇcme e vektor sam´ ych jedniˇcek, tj. e = (1, 1, . . . , 1)T . Vˇ eta 2.2.1. Nez´aporn´a matice A je stochastick´a pr´avˇe tehdy, kdyˇz e je vlastn´ı vektor matice A pˇr´ısluˇsn´y vlastn´ı hodnotˇe 1. P D˚ u k a z. Nejprvˇe mˇejme nez´apornou, stochastickou matici A. Z definice plat´ı j aij = 1 a tuto sumu m˚ uˇzeme zapsat Ae = e. Pak e je vlastn´ı vektor matice A pˇr´ısluˇsn´ y vlastn´ımu ˇc´ıslu 1. Naopak mˇejme nez´apornou matici s vlastn´ım vektorem e pˇr´ısluˇsn´ ym vlastn´ımu ˇc´ıslu P 1, tj. Ae = e. Z toho plyne aij ≥ 0 a j aij = 1, coˇz jsou podm´ınky z definice stochastick´e matice. Druhou podm´ınku z definice stochastick´e matice lze tedy podle d˚ ukazu Vˇety 2.2.1 zapsat tak´e takto: Ae = e a podle Vˇety 2.1.8 je e Perron˚ uv vektor nez´aporn´e matice A.
24
Vˇ eta 2.2.2. Je-li P nez´aporn´a matice, jej´ıˇz Perron˚ uv vektor u je kladn´y a ρ(P ) > 0, pak existuje diagon´aln´ı matice D s kladn´ymi diagon´ aln´ımi prvky a ˇc´ıslo k > 0 tak, ˇze matice A = kDP D−1 je stochastick´a. Pˇritom k = (ρ(P ))−1 . D˚ u k a z. Necht’ P u = pu, u > 0. Existuje tedy diagon´aln´ı matice D s kladn´ ymi T −1 diagon´aln´ımi prvky tak, ˇze Du = e = (1, . . . , 1) . Necht’ k = (ρ(P )) . Pak Ae = kDP D−1 e = kDP D−1 Du = kDP u = kDρ(P )u = Du = e, a pˇritom A ≥ 0. Volba k je jednoznaˇcn´a, nebot’ 1 = ρ(A) = kρ(DP D−1 ) = kρ(P ). Vˇ eta 2.2.3. Souˇcin stochastick´ych matic t´ehoˇz ˇr´adu je opˇet stochastick´a matice. D˚ u k a z. Necht’ A a B jsou stochastick´e matice t´ehoˇz ˇr´adu. Protoˇze aij ≥ 0,
i, j = 1, . . . , n,
bkl ≥ 0,
k, l = 1, . . . , n,
pak bude platit i aij bkl ≥ 0. Druhou podm´ınku z definice stochastick´e matice jsme zapsali Ae = e a tedy Be = e. Pak pro souˇcin matic plat´ı ABe = Ae = e a t´ım je dok´az´ano, ˇze AB je stochastick´a matice.
2.2.1
Dvojitˇ e stochastick´ e matice
Dalˇs´ı d˚ uleˇzitou podtˇr´ıdou nez´aporn´ ych matice jsou tzv. dvojitˇe stochastick´e matice. ˇ Definice 2.2.2. Ctvercov´ a matice A se naz´yv´a dvojitˇe stochastick´a, jsou-li stochasT tick´e matice A i A , tj. jestliˇze plat´ı: X X aik ≥ 0, aik = 1, aik = 1, i
k
neboli A ≥ 0,
Ae = e,
AT e = e.
O dvojitˇe stochastick´ ych matic´ıch plat´ı nˇekter´e vˇety. Vˇ eta 2.2.4. Perronovo vlastn´ı ˇc´ıslo dvojitˇe stochastick´e matice je 1 a odpov´ıd´a mu jak u A, tak u AT vlastn´ı vektor e. D˚ u k a z. Je zˇrejm´ y ze z´apisu Ae = e, AT e = e. 25
Vˇ eta 2.2.5. Souˇcin dvojitˇe stochastick´ych matic je opˇet dvojitˇe stochastick´a matice. D˚ u k a z. Mˇejme A, B dvojitˇe stochastick´e matice. Podobnˇe jako v analogick´e vˇetˇe o stochastick´ ych matic´ıch je aij bkl ≥ 0 a ABe = e. Zb´ yv´a dok´azat, ˇze (AB)T e = e. (AB)T e = B T AT e = B T e = e.
Vˇ eta 2.2.6. Je-li A dvojitˇe stochastick´a matice n-t´eho ˇr´adu, existuje permutace (i1 , . . ., in ) index˚ u 1, . . . , n takov´ a, ˇze a1i1 a2i2 . . . anin > 0. Jin´ymi slovy, existuje nenulov´y ˇclen v rozvoji det A z Definice determinantu 1.1.3. D˚ u k a z. Oznaˇcme P1 , . . ., PN pro N = n! permutaˇcn´ı matice ˇr´adu n. Matice Pi , i = 1, . . . , N jsou zˇrejmˇe dvojitˇe stochastick´e a kaˇzdou dvojitˇe stochastickou matici A tedy lze vyj´adˇrit ve tvaru A=
N X
λi ≥ 0,
λi Pi ,
i=1
N X
λi = 1.
i=1
V tomto vyj´adˇren´ı je aspoˇ n jeden koeficient λi kladn´ y. M´a-li pˇr´ısluˇsn´a matice Pi jedniˇcky na m´ıstech (1, i1 ), (2, i2 ), . . ., (n, in ), je ak,ik > 0 pro k = 1, . . . , n a plat´ı tedy dokazovan´e tvrzen´ı. Nyn´ı kaˇzd´e re´aln´e ˇctvercov´e matici A = (aik ) n-t´eho ˇr´adu pˇriˇrad´ıme bod euklie o souˇradnic´ıch (a11 , a12 , . . ., a1n , a21 , a22 , . . ., dovsk´eho n2 -rozmˇern´eho prostoru E a2n , . . ., an1 , . . . , ann ). Vˇ eta 2.2.7 (Birkhoffova vˇeta). Mnoˇzinˇe dvojitˇe stochastick´ych matice n-t´eho ˇr´adu e mnohostˇen, jehoˇz vrcholy jsou pr´avˇe vˇsechny body odpov´ıdaj´ıc´ı odpov´ıd´a v prostoru E permutaˇcn´ım matic´ım. Mnohostˇen je mnoˇzina bod˚ u tvaru N X
λi Ai ,
i=1
kde λi ≥ 0,
PN
i=1
λi = 1 a Ai jsou matice.
D˚ u k a z. viz [3]
26
2.3
Symetrick´ e a hermitovsk´ e matice
Definice 2.3.1. Matice A se naz´yv´a symterick´a, plat´ı-li AT = A neboli plat´ı-li pro A = (aik ), ˇze aik = aki pro vˇsechna i, k. Matice A je hermitovsk´ a, je-li A∗ = A neboli plat´ı pro A = (aik ), ˇze aik = a ¯ki pro vˇsechna i, k. Kaˇzd´e ˇctvercov´e matici A m˚ uˇzeme pˇriˇradit hermitovskou matici 1 (A + A∗ ); 2 t´eto matici budeme ˇr´ıkat symetrick´a ˇc´ast matice A a oznaˇcovat ji ReA. Pˇ r´ıklad 2.3.1. Symetrick´a matice tˇret´ıho ˇr´adu: 1 2 3+i −i 5 . 2 3+i 5 −6 Hermitovsk´a matice tˇret´ıho ˇr´adu:
2 − 31 2 − 3i 5 i − 31 . 2 + 3i −i −4
Definice 2.3.2. Matice A se naz´yv´a ortogon´aln´ı, plat´ı-li AAT = I. Matice A se naz´yv´a unit´arn´ı, plat´ı-li AA∗ = I. Nyn´ı uvedeme vˇetu, kterou budeme potˇrebovat v d˚ ukazu hlavn´ıch vˇet o symetrick´ ych a hermitovsk´ ych matic´ıch. Vˇ eta 2.3.1. Vlastn´ı ˇc´ısla hermitovsk´e (tedy i symetrick´e) matice jsou re´aln´a. D˚ u k a z. Necht’ λ je vlastn´ı ˇc´ıslo matice A a u odpov´ıdaj´ıc´ı vlastn´ı vektor. Pak z vlastnosti skal´arn´ıho souˇcinu hAu, ui = hu, A∗ ui = hu, Aui = hAu, ui je hAu, ui re´aln´e. Protoˇze Au = λu, je hAu, ui = hλu, ui = λhu, ui. Nebot’ hu, ui je re´aln´e, je i λ re´aln´e. Vˇ eta 2.3.2 (hlavn´ı vˇeta o symetrick´ ych matic´ıch). Kaˇzdou re´alnou symetrickou matici A lze vyj´adˇrit ve tvaru A = CDC T , kde C je ortogon´aln´ı a D re´aln´a diagon´ aln´ı matice. Pˇritom diagon´ aln´ı prvky matice D jsou vlastn´ı ˇc´ısla matice A, sloupce matice C vlastn´ı vektory matice A (k-t´y sloupec matice C odpov´ıd´a k-t´emu diagon´ aln´ımu prvku matice D). 27
D˚ u k a z. Matice A m´a podle Vˇety 2.3.1 re´aln´a vlastn´ı ˇc´ısla. Proto m˚ uˇzeme A T vyj´adˇrit ve tvaru A = CT C , kde C je ortogon´aln´ı a T re´aln´a horn´ı troj´ uheln´ıkov´a matice. Protoˇze A je symetrick´a, je T tak´e symetrick´a. Proto je matice T diagon´aln´ı, oznaˇcme ji D = diag{λ1 , . . . , λn }. Piˇsme matici C v blokov´em z´apisu C = (C1 , . . . , Cn ), kde n je ˇr´ad matice C a Ck je k-t´ y sloupec. T Ze vztahu A = CDC plyne AC = CD neboli λ1 0 . . . 0 0 λ2 . . . 0 A(C1 , . . . , Cn ) = (C1 , . . . , Cn ) . .. . . .. . . . . . . 0
0
. . . λn
Podle vˇety o n´asoben´ı blokov´ ych matic je tedy
(AC1 , . . . , ACn ) = (λ1 C1 , . . . , λn Cn ) neboli AC1 = λ1 C1 , . . . , ACn = λn Cn . k-t´ y sloupec Ck matice C je tedy skuteˇcnˇe vlastn´ım vektorem matice A odpov´ıdaj´ıc´ım k-t´emu diagon´aln´ımu prvku matice D. Vˇ eta 2.3.3 (hlavn´ı vˇeta o hermitovsk´ ych matic´ıch). Kaˇzdou hermitovskou matici A lze vyj´adˇrit ve tvaru A = U DU ∗ , kde U je unit´arn´ı a D re´aln´a diagon´ aln´ı matice. Pˇritom diagon´ aln´ı prvky matice D jsou vlastn´ı ˇc´ısla matice A, sloupce U vlastn´ı vektory matice A (k-t´y sloupec matice U odpov´ıd´a k-t´emu diagon´aln´ımu prvku matice D). D˚ u k a z. Existuje unit´arn´ı matice U a horn´ı troj´ uheln´ıkov´a matice D tak, ˇze ∗ A = U DU . Protoˇze matice A je hermitovsk´a, je tak´e D hermitovsk´a. Proto je D diagon´aln´ı, a protoˇze je kaˇzd´ y diagon´aln´ı prvek matice D roven sv´emu komplexnˇe sdruˇzen´emu ˇc´ıslu (z definice hermitovsk´e matice), je D re´aln´a diagon´aln´ı matice. Oznaˇc´ıme-li Uk k-t´ y sloupec matice U a n ˇr´ad matice A, je U = (U1 , . . . , Un ). Necht’ D = diag{λ1 , . . . , λn }. Ze vztahu A = U DU ∗ plyne AU = U D neboli v blokov´em tvaru λ1 0 . . . 0 0 λ2 . . . 0 A(U1 , . . . , Un ) = (U1 , . . . , Un ) .. .. . . . .. . . . . 0 0 . . . λn Proto
(AU1 , . . . , AUn ) = (λ1 U1 , . . . , λn Un ), tj. AU1 = λ1 U1 , . . . , AUn = λn Un . Uk je tedy vlastn´ı vektor matice A odpov´ıdaj´ıc´ı vlastn´ımu ˇc´ıslu λk . 28
2.4
M-matice
Na zaˇc´atek uvedeme definici pozitivnˇe definitn´ı a semidefinitn´ı matice. Budeme ji v dalˇs´ım potˇrebovat. Definice 2.4.1. Hermitovsk´a (popˇr. symetrick´a) matice A = (aik ) se naz´yv´a pozitivnˇe definitn´ı, plat´ı-li pro kaˇzd´y nenulov´y komplexn´ı (re´aln´y) vektor x, ˇze hAx, xi > 0 neboli
X
aik x¯i xk > 0.
i,k
Hermitovsk´a (popˇr. symetrick´a) matice A = (aik ) se naz´yv´a pozitivnˇe semidefinitn´ı, plat´ı-li pro kaˇzd´y komplexn´ı (re´aln´y) vektor x, ˇze hAx, xi ≥ 0 neboli
X i,k
aik x¯i xk ≥ 0.
Nejdˇr´ıve zavedeme toto oznaˇcen´ı: Tˇr´ıdou Zn (n ≥ 1) budeme rozumˇet mnoˇzinu vˇsech ˇctvercov´ ych re´aln´ ych matic n-t´eho ˇr´adu, jejichˇz vˇsechny nediagon´aln´ı prvky jsou nekladn´e: Zn = {A = (aik ), typu(n, n); aik ≤ 0, i 6= k}. Sjednocen´ı vˇsech mnoˇzin Zn oznaˇc´ıme Z: [ Z= Zn . n=1,2,...
Definice 2.4.2. Matice tˇr´ıdy K jsou matice A v Z, kter´e (pro pˇr´ısluˇsn´e n) splˇ nuj´ı jednu z ekvivalentn´ıch podm´ınek: 1. existuje vektor x ≥ 0 tak, ˇze Ax > 0; 2. existuje vektor x > 0 tak, ˇze Ax > 0; 3. existuje diagon´aln´ı matice D s kladn´ymi diagon´ aln´ımi prvky tak, ˇze matice P AD = (wik ) splˇ nuje wii > k6=i |wik | pro vˇsechna i; 4. kdykoliv je B ∈ Zn a B ≥ A, je B regul´arn´ı matice;
5. kaˇzd´e re´aln´e vlastn´ı ˇc´ıslo kaˇzd´e hlavn´ı podmatice matice A je kladn´e; 6. vˇsechny hlavn´ı minory matice A jsou kladn´e; 7. souˇcty vˇsech hlavn´ıch minor˚ u k-t´eho ˇr´adu matice A jsou kladn´e pro k = 1, . . . , n; 29
8. kaˇzd´e re´aln´e vlastn´ı ˇc´ıslo matice A je kladn´e; 9. existuje matice C ≥ 0 a ˇc´ıslo k > ρ(C) tak, ˇze A = kI − C; 10. existuje rozklad A = P −Q matice A takov´y, ˇze P −1 ≥ 0, Q ≥ 0 a ρ(P −1 Q) < 1; 11. A je regul´arn´ı matice a A−1 ≥ 0; 12. detA(Nk ) > 0, k = 1, . . . , n, kde Nk = {1, 2, . . . , k} a detA(Nk ) jsou hlavn´ı minory matice A; 13. existuje doln´ı troj´ uheln´ıkov´a matice R a horn´ı troj´ uheln´ıkov´a matice S, obˇe s kladn´ymi diagon´aln´ımi prvky, tak, ˇze A = RS; 14. existuje doln´ı troj´ uheln´ıkov´a matice R ∈ Zn a horn´ı troj´ uheln´ıkov´a matice S ∈ Zn , obˇe s kladn´ymi diagon´ aln´ımi prvky, tak, ˇze A = RS; 15. pro kaˇzd´y rozklad A = P − Q matice A slpˇ nuj´ıc´ı P −1 ≥ 0 a Q ≥ 0 plat´ı ρ(P −1 Q) < 1; 16. existuje diagon´aln´ı matice D s kladn´ymi diagon´ aln´ımi prvky tak, ˇze matice 1 −1 T B = DAD m´a symetrickou ˇc´ast 2 (B + B ) pozitivnˇe definitn´ı; 17. existuje diagon´aln´ı matice H s kladn´ymi diagon´ aln´ımi prvky tak, ˇze matice 1 C = AH m´a symetrickou ˇc´ast 2 (C + C T ) pozitivnˇe definitn´ı; 18. re´aln´a ˇc´ast kaˇzd´eho vlastn´ıho ˇc´ısla matice A je kladn´a. Matice tˇr´ıdy K oznaˇcujeme tak´e jako M-matice. Definice 2.4.3. K0 nazveme tˇr´ıdu vˇsech matic v Z, kter´e splˇ nuj´ı (pro odpov´ıdaj´ıc´ı n) jednu z n´ asleduj´ıc´ıch ekvivalentn´ıch podm´ınek: 1. A + εI ∈ K pro vˇsechna ε > 0; 2. kaˇzd´e re´aln´e vlastn´ı ˇc´ıslo kaˇzd´e hlavn´ı podmatice matice A je nez´ aporn´e; 3. vˇsechny hlavn´ı minory matice A jsou nez´ aporn´e; 4. souˇcty vˇsech hlavn´ıch minor˚ u k-t´eho ˇr´adu matice A jsou nez´ aporn´e pro k = 1, . . . , n; 5. kaˇzd´e re´aln´e vlastn´ı ˇc´ıslo matice A je nez´ aporn´e; 6. existuje matice C ≥ 0 a ˇc´ıslo k ≥ ρ(C) tak, ˇze A = kI − C; 7. existuje diagon´aln´ı matice D s kladn´ymi diagon´ aln´ımi prvky tak, ˇze symetrick´a 1 T −1 ˇc´ast 2 (B + B ) matice B = DAD je pozitivnˇe semidefinitn´ı;
30
8. existuje diagon´aln´ı matice H s kladn´ymi diagon´ aln´ımi prvky takov´ a, ˇze syme1 T trick´a ˇc´ast 2 (C + C ) matice C = AH je pozitivnˇe semidefinitn´ı; 9. kaˇzd´e vlastn´ı ˇc´ıslo matice A m´a nez´ apornou re´alnou ˇc´ast. Vˇ eta 2.4.1. Plat´ı K ⊂ K0 . Matice z K0 je v K, pr´avˇe kdyˇz je regul´arn´ı. D˚ u k a z. Necht’ A ∈ K. Pak A ∈ Z a A splˇ nuje vlastnost 6. z Definice 2.4.2. Pak vˇsak splˇ nuje vlastnost 3. z Definice 2.4.3, tj. A ∈ K0 . Necht’ nyn´ı A ∈ K0 . Je-li A singul´arn´ı, nen´ı A ∈ K podle vlastnosti 11. Definice 2.4.2. Je-li A regul´arn´ı, pak nem´a vlastn´ı ˇc´ıslo 0. Podle 5. z Definice 2.4.3 je kaˇzd´e re´aln´e vlastn´ı ˇc´ıslo matice A nez´aporn´e, tedy kladn´e. Proto A splˇ nuje 8. z Definice 2.4.2 a A ∈ K.
2.5
Stabiln´ı matice
V t´eto ˇc´asti se budeme zab´ yvat tzv. stabiln´ımi maticemi. Ty se pouˇz´ıvaj´ı v aplikac´ıch napˇr. pˇri ˇreˇsen´ı soustav diferenci´aln´ıch rovnic nebo pˇri ˇreˇsen´ı r˚ uzn´ ych ekonomick´ ych probl´em˚ u. ˇ Definice 2.5.1. Ctvercov´ a matice A se naz´yv´a stabiln´ı, maj´ı-li vˇsechna vlastn´ı ˇc´ısla matice A z´ apornou re´alnou ˇc´ast, a pozitivnˇe stabiln´ı, maj´ı-li vˇsechna vlastn´ı ˇc´ısla matice A kladnou re´alnou ˇc´ast. Pozn´ amka 2.5.1. Zˇrejmˇe tedy A je stabiln´ı, pr´avˇe kdyˇz −A je pozitivnˇe stabiln´ı. Charakterizace stabiln´ıch matice je uvedena v n´asleduj´ıc´ı vˇetˇe: Vˇ eta 2.5.1 (hlavn´ı vˇeta o pozitivnˇe stabiln´ıch matic´ıch). Necht’ A je (obecnˇe komplexn´ı) ˇctvercov´a matice. Oznaˇcme B = (A+I)−1 (A−I). Potom jsou tyto vlastnosti matice A navz´ajem ekvivalentn´ı: 1. matice A je pozitivnˇe stabiln´ı; 2. matice B existuje a plat´ı ρ(B) < 1; 3. matice B existuje a pro kaˇzdou pozitivnˇe definitn´ı (hermitovskou) matici C existuje (hermitovsk´ a) pozitivnˇe definitn´ı matice P tak, ˇze P − BP B ∗ = C (B ∗ je konjugovan´a, tj. transponovan´a a komplexnˇe sdruˇzen´ a, matice k matici B); 4. matice B existuje a pro kaˇzdou pozitivnˇe definitn´ı matici C ˇrada C + BCB ∗ + B 2 C(B ∗ )2 + B 3 C(B ∗ )3 + . . . konverguje;
31
5. matice B existuje a ˇrada I + BB ∗ + B 2 (B ∗ )2 + B 3 (B ∗ )3 + . . . konverguje; 6. matice B existuje a existuje pozitivnˇe definitn´ı matice Q tak, ˇze Q − BQB ∗ = I; 7. matice B existuje a existuje pozitivnˇe definitn´ı matice R tak, ˇze R − BRB ∗ je pozitivnˇe definitn´ı; 8. existuje pozitivnˇe definitn´ı matice S tak, ˇze AS + SA∗ je pozitivnˇe definitn´ı; 9. pro kaˇzdou pozitivnˇe definitn´ı matici G existuje pozitivnˇe definitn´ı matice V tak, ˇze AV + V A∗ = G; 10. existuje pozitivnˇe definitn´ı matice V tak, ˇze AV + V A∗ = I; 11. existuje pozitivnˇe definitn´ı matice W tak, ˇze matice W AW −1 + W −1 A∗ W je pozitivnˇe definitn´ı; 12. existuje regul´arn´ı matice T tak, ˇze symetrick´a ˇc´ast matice T AT −1 , tj. 1 Re(T AT −1 ) = (T AT −1 + (T ∗ )−1 A∗ T ∗ ), 2 je pozitivnˇe definitn´ı. D˚ u k a z. viz [3]
2.6
P´ asov´ e matice
V numerick´e matematice se pˇri ˇreˇsen´ı obyˇcejn´ ych a parci´aln´ıch diferenci´aln´ıch rovnic diskretizaˇcn´ımi metodami a pˇri numerick´em ˇreˇsen´ı probl´emu vlastn´ıch ˇc´ısel matic ˇcasto vyskytuj´ı tzv. p´asov´e matice. Je-li A = (aik ) ˇctvercov´a matice n-t´eho ˇr´adu, oznaˇcme p = max(p0 , 0), q = max(q0 , 0),
p0 = max{k − i; i, k, aik 6= 0}, q0 = − min{k − i; i, k, aik 6= 0}.
ˇ ısla p a q jsou ˇs´ıˇrky minim´aln´ıho naddiagon´aln´ıho a poddiagon´aln´ıho p´asu, v nˇemˇz C´ jsou obsaˇzeny vˇsechny nenulov´e nediagon´aln´ı prvky matice. 32
ˇ Definice 2.6.1. Rekneme, ˇze matice A je (p, q)-p´asov´a, jsou-li p a q ˇc´ısla definovan´ a v´yˇse. D˚ uleˇzit´e speci´aln´ı pˇr´ıpady jsou doln´ı troj´ uheln´ıkov´a matice, pro kterou je p = 0, horn´ı troj´ uheln´ıkov´a matice, pro kterou je q = 0 a diagon´aln´ı matice, pro niˇz plat´ı p = 0, q = 0. Pˇ r´ıklad 2.6.1. Matice (1, 1)-p´asov´a je tˇr´ıdiagon´aln´ı matice, matice (2, 2)-p´asov´a je pˇetidiagon´aln´ı matice, matice pro kterou je q = 1 je matice v tzv. Hessenbergovˇe tvaru. V n´asleduj´ıc´ıch vˇet´ach jsou uvedeny z´akladn´ı vlastnosti a hlavn´ı v´ yznam p´asov´ ych matic. Vˇ eta 2.6.1. Necht’ A1 , A2 jsou ˇctvercov´e matice n-t´eho ˇr´adu, necht’ A1 je (p1 , q1 )p´asov´a, A2 je (p2 , q2 )-p´asov´a matice. Pak AT1 je (q1 , p1 )-p´asov´a, A1 + A2 je (p3 , q3 )p´asov´a, kde p3 ≤ max(p1 , p2 ), q3 ≤ max(q1 , q2 ), a A1 A2 je (p4 , q4 )-p´asov´a, kde p4 ≤ min(n − 1, p1 + p2 ),
q4 ≤ min(n − 1, q1 + q2 ).
D˚ u k a z. Prvn´ı dvˇe tvrzen´ı jsou zˇrejm´a. Tˇret´ı plyne z toho, ˇze o ˇc´ıslech p40 , p10 , p20 z definice p, q v zaˇc´atku kapitoly plat´ı p40 = max{k − i; i, k,
n X j=1
(2)
(1) (2)
aij ajk 6= 0} ≤ (1)
≤ max{k − j; j, k, ajk 6= 0} + max{j − i; i, j, aij 6= 0} = p10 + p20 , a ovˇsem tak´e p40 ≤ n − 1, a obdobnˇe pro q4 . Vˇ eta 2.6.2. Necht’ A je ˇctvercov´a silnˇe regul´arn´ı (p, q)-p´asov´a matice. Pak existuje rozklad A = BC matice A na souˇcin (0, q)-p´asov´e matice B a (p, 0)-p´asov´e matice C. D˚ u k a z. Existuje rozklad silnˇe regul´arn´ı matice A na souˇcin BC doln´ı troj´ uheln´ıkov´e matice B a horn´ı troj´ uheln´ıkov´e matice C. Necht’ tato matice B je (0, q1 )p´asov´a, matice C (p1 , 0)-p´asov´a. Podle Vˇety 2.6.1 je p ≤ p1 , q ≤ q1 . Staˇc´ı dok´azat, ˇze p = p1 , q = q1 . Oznaˇc´ıme-li A = (aik ), B = (bik ), C = (cik ), plat´ı bii 6= 0, cii 6= 0, i = 1, . . . , n, nebot’ obˇe matice B a C jsou regul´arn´ı. Proto je p1 = max{k − i; i, k, cik 6= 0}. Existuje tedy prvek cik 6= 0 tak, ˇze i ≤ k, a pˇritom p1 = k − i, a tak´e cjk = 0 pro vˇsechna j, 1 ≤ j < i. Plat´ı pak pro tuto dvojici i, k aik =
i X j=1
bij cjk = bii cik 6= 0. 33
Proto p ≥ k − i = p1 , takˇze skuteˇcnˇe p = p1 . Obdobnˇe je q1 = − min{k − i; i, k, bik 6= 0}. Existuje proto prvek bik 6= 0, i ≥ k, takov´ y, ˇze q1 = i − k, a pˇritom bij = 0 pro vˇsechna j, 1 ≤ j < k. Pak aik =
k X j=1
bij cjk = bik ckk 6= 0,
odkud q ≥ i − k = q1 , tj. q = q1 . Vˇ eta 2.6.3. Necht’
A=
α1 γ1 0 . . . 0 0 γ1 α2 γ2 . . . 0 0 0 γ2 α3 . . . 0 0 .. .. .. . . . . . .. . . . . . 0 0 0 . . . αn−1 γn−1 0 0 0 . . . γn−1 αn
a ξ ∈ R. Utvoˇrme rekurentnˇe posloupnost δ0 , δ1 , . . ., δn takto: δ0 = 1,
δ1 = α1 − ξ,
2 δk = (αk − ξ)δk−1 − γk−1 δk−2 ,
k = 2, . . . , n.
Nen´ı-li ˇza´dn´e z ˇc´ısel δ0 , . . ., δn rovno nule, m´a matice A tolik vlastn´ıch ˇc´ısel vˇetˇs´ıch neˇz ξ, kolik je znam´enkov´ych shod v posloupnosti δ0 , δ1 , . . ., δn . Pozn´ amka 2.6.1. Znam´enkov´a shoda nastane, kdyˇz plat´ı sgnδk = sgnδk+1 . D˚ u k a z. (vˇety) Matice A m´a zˇrejmˇe tolik vlastn´ıch ˇc´ısel vˇetˇs´ıch neˇz ξ, kolik m´a matice A − ξI kladn´ ych vlastn´ıch ˇc´ısel. Tento poˇcet je u symetrick´e matice P roven poˇctu znam´enkov´ ych shod v posloupnosti 1, detP (N1 ), . . . , detP (Nn ), kde detP (Nk ), Nk = {1, 2, . . . , k}, k = 1, . . . , n jsou hlavn´ı minory matice P . T´eto vˇety lze vyuˇz´ıt k nalezen´ı vlastn´ıch ˇc´ısel tˇr´ıdiagon´aln´ı symetrick´e matice metodou p˚ ulen´ı interval˚ u. Nejdˇr´ıve najdeme interval [α, β], v nˇemˇz leˇz´ı vˇsechna vlastn´ı ˇc´ısla matice A. Pak zvol´ıme ξ = 21 (α + β) a zjist´ıme, kolik m´a matice A vlastn´ıch ˇc´ısel v intervalu [ξ, β] (a tedy tak´e kolik v intervalu [α, ξ]). Pokud je v nˇekter´em z interval˚ u alespoˇ n jedno vlastn´ı ˇc´ıslo, interval se opˇet rozp˚ ul´ı. Pokraˇcujeme tak dlouho, dokud nedos´ahneme poˇzadovan´e pˇresnosti aproximace vlastn´ıho ˇc´ısla. Vˇ eta 2.6.4. Re´aln´a tˇr´ıdiagon´aln´ı matice A = (aik ) n-t´eho ˇr´adu splˇ nuj´ıc´ı ak,k+1 ak+1,k > 0 k = 1, . . . , n − 1, m´a n re´aln´ych jednoduch´ych vlastn´ıch ˇc´ısel.
34
D˚ u k a z. Provedeme indukc´ı vzhledem k n. Pro n = 1 vˇeta plat´ı. Necht’ n > 1 a pˇredpokl´adejme, ˇze vˇeta plat´ı pro matice ˇr´adu n − 1. Matici A lze pˇrev´est na matici A˜ tvaru α1 −1 . . . 0 0 0 0 −β1 α2 . . . . .. ... ... ... .. . , 0 0 . . . αn−1 −1 0 0 . . . −βn−1 αn
pro kterou βi > 0, i = 1, . . . , n−1. Protoˇze A˜ ∈ Z, existuje ˇc´ıslo k tak, ˇze kI +A˜ ∈ K. ˜ jemuˇz odpov´ıd´a kladn´ Pak existuje re´aln´e jednoduch´e vlastn´ı ˇc´ıslo matice kI + A, y vlastn´ı vektor. Pro matici A to znamen´a, ˇze existuje jej´ı jednoduch´e re´aln´e vlastn´ı ˜ ˇc´ıslo ω, jemuˇz odpov´ıd´a kladn´ y vlastn´ı vektor y. Nyn´ı dostaneme, ˇze matice A, kter´a uˇz nem´a vlastn´ı ˇc´ıslo ω, opˇet splˇ nuje pˇredpoklady vˇety a m´a ˇr´ad n − 1. Podle indukˇcn´ıho pˇredpokladu m´a matice A˜ n−1 re´aln´ ych jednoduch´ ych vlastn´ıch ˇc´ısel. To jsou vˇsak i vlastn´ı ˇc´ısla matice A, kter´a proto m´a spolu s ω n re´aln´ ych jednoduch´ ych vlastn´ıch ˇc´ısel. Vˇ eta 2.6.5. Necht’ A = (aik ) je re´aln´a tˇr´ıdiagon´aln´ı matice n-t´eho ˇr´adu, n ≥ 2, splˇ nuj´ıc´ı ai,i+1 ai+1,i > 0. Jsou-li λ1 > λ2 > . . . > λn re´aln´a vlastn´ı ˇc´ısla matice A, pak pro kaˇzd´y re´aln´y vlastn´ı vektor z = (z1 , . . . , zn )T matice A plat´ı 1. z1 6= 0, zn 6= 0; 2. je-li zk = 0, pak ak−1,k ak,k+1 zk−1 zk+1 < 0; 3. vynech´ame-li v posloupnosti z1 , a12 z2 , a12 a23 z3 , . . ., a12 a23 . . . an−1,n zn nuly, pak poˇcet znam´enkov´ych zmˇen je pr´avˇe r − 1, jestliˇze vektor z odpov´ıd´a vlastn´ımu ˇc´ıslu λr . D˚ u k a z. viz [3]
35
Kapitola 3 Iteraˇ cn´ı metody ˇ reˇ sen´ı soustav line´ arn´ıch rovnic 3.1
Iteraˇ cn´ı metody
V numerick´e matematice vˇetˇsinou nevystaˇc´ıme s algoritmy, kter´e d´avaj´ı v´ ysledky ve tvaru vzorce. To se st´av´a napˇr´ıklad v pˇr´ıpadˇe, ˇze v soustavˇe rovnic je poˇcet rovnic a poˇcet nezn´am´ ych velik´ y a uˇzit´ı pˇr´ım´ ych metod by bylo pracn´e a velmi sloˇzit´e. Nˇekdy je tedy tˇreba uˇz´ıt takov´eho numerick´eho v´ ypoˇctu, kter´ y umoˇzn ˇuje naj´ıt pˇribliˇznou hodnotu v´ ysledku pˇreruˇsen´ım v´ ypoˇctu ve vhodn´em okamˇziku. Iteraˇcn´ımi metodami ˇreˇs´ıme rovnice nebo soustavu rovnic tvaru x = f (x), kde x znamen´a nezn´amou, popˇr. vektor a f funkci, popˇr. vektorovou funkci. Zvol´ıme poˇc´ateˇcn´ı aproximaci x0 a poˇc´ıt´ame postupnˇe x1 , x2 , x3 , . . . ze vzorce xk+1 = f (xk ),
k = 0, 1, . . . .
Konverguje-li posloupnost ˇc´ısel, popˇr. vektor˚ u xk , je limita (je-li funkce f v okol´ı limity spojit´a) ˇreˇsen´ım dan´eho probl´emu. Probl´em obvykle je naj´ıt takovou u ´pravu rovnice g(x) = 0 na tvar x = f (x), aby uveden´a iteraˇcn´ı metoda d´avala konvergentn´ı posloupnost a ta aby nav´ıc konvergovala dostateˇcnˇe rychle. Nejdˇr´ıve si definujme r˚ uzn´e typy norem matice, kter´e budeme potˇrebovat k urˇcen´ı, zda konkr´etn´ı iteraˇcn´ı metoda konverguje. Definice 3.1.1. Necht’ je d´ana matice a11 . . . a1n . . A = .. . . . .. . an1 . . . ann
Potom definujeme ˇr´adkovou normu matice: kAk∞ = max i
36
n X j=1
|aij |,
sloupcovou normu matice: kAk1 = max j
a euklidovskou normu matice:
n X i=1
|aij |
v uX u n kAk2 = t |aij |2 . i,j=1
Nyn´ı uvaˇzujme syst´em line´arn´ıch rovnic
Ax = b, kde A je regul´arn´ı matice. Z´akladn´ı princip iteraˇcn´ıch metod je ve vyj´adˇren´ı tohoto syst´emu ve tvaru x = T x + g. Definice 3.1.2. Necht’ x0 ∈ Rn je poˇc´ateˇcn´ı aproximace. Pak posloupnost {xk }∞ k=0 urˇcen´a vztahem xk+1 = T xk + g, k = 0, 1, . . . se naz´yv´a iteraˇcn´ı posloupnost a matice T iteraˇcn´ı matice. Oznaˇcme x∗ ˇreˇsen´ı soustavy Ax = b. To bude pr´avˇe tehdy, kdyˇz x∗ je ˇreˇsen´ım syst´emu x = T x + g, tj. x∗ = (I − T )−1 g, je-li (I − T ) regul´arn´ı. Vˇ eta 3.1.1. H je konvergentn´ı matice ⇔ limk→∞ kH k k = 0 pro nˇejakou pˇridruˇzenou maticovou normu ⇔ ρ(H) < 1 ⇔ limk→∞ H k x = o pro x ∈ Rn . D˚ u k a z. viz [4] Lemma 3.1.1. Necht’ ρ(T ) < 1. Pak (I-T) je regul´arn´ı a plat´ı (I − T )−1 = I + T + T 2 + . . . . D˚ u k a z. viz [4] Vˇ eta 3.1.2. Iteraˇcn´ı posloupnost {xk }∞ cen´a procesem xk+1 = T xk + g konverk=0 urˇ guje pro kaˇzdou poˇc´ateˇcn´ı aproximaci x0 ∈ Rn pr´avˇe tehdy, kdyˇz ρ(T ) < 1, pˇriˇcemˇz plat´ı lim xk = x∗ , x∗ = T x∗ + g. k→∞
D˚ u k a z. Necht’ x0 ∈ Rn libovolnˇe, pak xk+1 = T xk + g = T (T xk−1 + g) + g = . . . = T k+1 x0 + (T k + . . . + T + I)g. 37
Necht’ ρ(T ) < 1, pak T je konvergentn´ı a (I − T ) regul´arn´ı. Odtud plyne, ˇze lim xk+1 = lim T k+1 x0 + lim (T k + . . . + I)g = o + (I − T )−1 g = x∗ .
k→∞
k→∞
∗
k→∞
∗
Nyn´ı pro x = T x + g je x∗ − xk = T (x∗ − xk−1 ) = . . . = T k (x∗ − x0 ). Odtud lim (x∗ − xk ) = lim T k (x∗ − x0 ) = o.
k→∞
k→∞
Je-li z ∈ R libovoln´ y vektor a poloˇz´ıme-li x0 = x∗ − z, pak n
lim T k z = lim T k (x∗ − (x∗ − z)) = o,
k→∞
k→∞
coˇz implikuje ρ(T ) < 1. Pro odhad pˇresnosti a poˇctu iteraˇcn´ıch krok˚ u je moˇzn´e pouˇz´ıt n´asleduj´ıc´ıch vztah˚ u: 1. Je-li x∗ hledan´e ˇreˇsen´ı, potom kx∗ − xk k ≤ 2. je-li x∗ hledan´e ˇreˇsen´ı, potom kx∗ − xk k ≤
3.2
kT k kxk 1−kT k kT kk kx1 1−kT k
− xk−1 k, − x0 k.
Jacobiho metoda (metoda prost´ e iterace)
Matici syst´emu Ax = b pˇrepiˇsme ve tvaru A = D − L − U, kde D je diagon´aln´ı, L doln´ı troj´ uheln´ıkov´a a U horn´ı troj´ uheln´ıkov´a matice. L a U jsou tedy matice s nulami na diagon´ale. Pak soustavu Ax = b m˚ uˇzeme pˇrepsat (D − L − U )x = b Dx = (L + U )x + b. Je-li matice D regul´arn´ı, tj. aii 6= 0, i = 1, . . . , n, je x = D−1 (L + U )x + D−1 b. Oznaˇc´ıme-li TJ = D−1 (L + U ) a gJ = D−1 b, pak proces xk+1 = TJ xk + gJ naz´ yv´ame Jacobiho iteraˇcn´ı metoda. Protoˇze iteraˇcn´ı matice TJ m´a zˇrejmˇe diagon´aln´ı prvky nulov´e, m˚ uˇzeme metodu zapsat po sloˇzk´ach: xk+1 i
=−
n X aij j=1 j6=i
aii
xkj +
bi , aii
i = 1, . . . , n,
k ≥ 0.
Z Vˇety 3.1.2 plyne nutn´a a postaˇcuj´ıc´ı podm´ınka pro konvergenci metody ρ(TJ ) < 1. 38
Vˇ eta 3.2.1. Je-li A ryze ˇr´adkovˇe (resp. sloupcovˇe) diagon´ alnˇe dominantn´ı, pak Jacobiho iteraˇcn´ı metoda konverguje pro libovolnou poˇc´ateˇcn´ı aproximaci. D˚ u k a z. Necht’ TJ = D−1 (L + U ) je iteraˇcn´ı matice Jacobiho metody. Pak kTJ k∞
n X
n X aij . = max |tij | = max 1≤i≤n 1≤i≤n aii j=1 j=1 j6=i
Pro ryze ˇr´adkovˇe diagon´alnˇe dominantn´ı matici A plat´ı |aii | >
n X j=1 j6=i
|aij | a tedy
n X aij < 1, aii j=1
i = 1, . . . , n.
Odtud plyne kTJ k < 1 a Jacobiho iteraˇcn´ı metoda konverguje. Nyn´ı mˇejme syst´em s ryze sloupcovˇe diagon´alnˇe dominantn´ı matic´ı A. Z prvn´ı ˇc´asti d˚ ukazu v´ıme, ˇze pro AT Jacobiho metoda konverguje a tedy ρ(D−1 (L + U )) < 1. Matice X = D−1 (LT + U T ) m´a stejn´a vlastn´ı ˇc´ısla jako X T = (L + U )D−1 a tedy tak´e jako matice podobn´a X T D−1 X T D = D−1 (L + U )D−1 D = D−1 (L + U ) = TJ . Odtud plyne, ˇze ρ(TJ ) < 1 a metoda konverguje.
3.3
Gauss-Seidelova metoda
Tato metoda je urˇcitou modifikac´ı Jacobiho metody. Vyuˇzijeme tedy tak´e rozkladu matice soustavy A = D − L − U , kde D, L a U jsou matice stejn´e jako v pˇredchoz´ı kapitole. Soustavu Ax = b pˇrep´ıˇseme na tvar (D − L − U )x = b (D − L)x = U x + b. Je-li opˇet aii 6= 0, i = 1, . . . , n, je matice (D − L) regul´arn´ı a x = (D − L)−1 U x + (D − L)−1 b. Oznaˇcme TGS = (D − L)−1 U a gGS = (D − L)−1 b. Pak proces xk+1 = TGS xk + gGS se naz´ yv´a Gauss-Seidelova iteraˇcn´ı metoda. Z Vˇety 3.1.2 dostaneme nutnou a postaˇcuj´ıc´ı podm´ınku konvergence ρ(TGS ) < 1. 39
Z maticov´eho z´apisu plyne, ˇze pro v´ ypoˇcet (k + 1)-n´ı iterace i-t´e sloˇzky, tj. xk+1 , i vyuˇzijeme jiˇz vypoˇc´ıtan´e (k +1)-n´ı iterace pˇredchoz´ıch sloˇzek. Z´apis metody po sloˇzk´ach tedy bude vypadat takto: xk+1 i
=−
i−1 X aij j=1
aii
xk+1 j
n X aij k bi − xj + , a aii j=i+1 ii
i = 1, . . . , n.
Pˇ r´ıklad 3.3.1. Ukaˇzte, ˇze pro soustavu 3x1 + 2x2 + 2x3 = 1 2x1 + 3x2 + 2x3 = 0 2x1 + 2x2 + 3x3 = −1 Jacobiho iteraˇcn´ı metoda nekonverguje, ale Gauss-Seidelova metoda ano. ˇ sen´ı: Nejdˇr´ıve matici soutavy A rozloˇz´ıme takto: Reˇ A = D − L − U, kde
3 0 0 D = 0 3 0 , 0 0 3
0 0 0 L = −2 0 0 −2 −2 0
Pro Jacobiho iteraˇcn´ı matici plat´ı:
0 −2 −2 a U = 0 0 −2 . 0 0 0
TJ = D−1 (L + U ), tj.
0 −2 −2 0 − 32 − 23 0 0 TJ = 0 31 0 −2 0 −2 = − 32 0 − 23 . − 32 − 23 0 0 0 13 −2 −2 0
1 3
Vlastn´ı ˇc´ısla t´eto matice jsou ˇreˇsen´ım rovnice
det(TJ − λI) = 0 a jsou to 2 , 3 4 = − . 3
λ1,2 = λ3
Spektr´aln´ı polomˇer je roven nejvˇetˇs´ımu (v absolutn´ı hodnotˇe) vlastn´ımu ˇc´ıslu, coˇz je zde 4 ρ(TJ ) = . 3 40
Nutn´a podm´ınka pro konvergenci tedy nen´ı splnˇena. Pro Gauss-Seidelovu metodu plat´ı: TGS = (D − L)−1 U, tj.
1 2
0
1 TGS = − 29 3 2 − 27 − 92
0 0 −2 −2 0 −1 −1 0 0 0 −2 = 0 49 − 92 . 4 16 0 27 − 13 0 0 0 27
Vlastn´ı ˇc´ısla opˇet dostaneme ˇreˇsen´ım rovnice
det(TGS − λI) = 0 a jsou to λ1 = 0, √ 14 2 5 ± i. λ2,3 = 27 27 Spektr´aln´ı polomˇer je tedy
√ 2 6 ρ(TGS ) = < 1. 9 To znamen´a, ˇze je splnˇena nutn´a a postaˇcuj´ıc´ı podm´ınka pro konvergenci. Pro p˚ uvodn´ı soustavu rovnic tedy Jacobiho metoda nekonverguje, ale GaussSeidelova metoda ano. Konvergence Gauss-Seidelovy metody je vidˇet i z Obr´azku 1 v Pˇr´ıloze. Pˇ r´ıklad 3.3.2. Ukaˇzte, ˇze pro soustavu x1 + x2 = 1 2(1 − ε)x1 + x2 + x3 = 2
x3 + x4 = −1
−(1 − ε)2 x1 + x4 = 5
Jacobiho iteraˇcn´ı metoda konverguje, ale Gauss-Seidelova metoda ne. (0 < ε < 0, 1) ˇ sen´ı: Nejdˇr´ıve rozloˇz´ıme matici soutavy A: Reˇ A = D − L − U, kde
D=
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
= I,
L= 41
0 −2(1 − ε) 0 (1 − ε)2
0 0 0 0
0 0 0 0
0 0 0 0
a
U =
0 −1 0 0 0 0 −1 0 0 0 0 −1 0 0 0 0
Pro Jacobiho iteraˇcn´ı matici plat´ı:
.
TJ = D−1 (L + U ), tj. v tomto pˇr´ıpadˇe TJ = I(L + U ) = L + U. Odtud je
TJ =
Vlastn´ı ˇc´ısla t´eto matice jsou
0 −1 0 0 2(1 − ε) 0 −1 0 0 0 0 −1 2 (1 − ε) 0 0 0
.
√
1 − ε, √ = − 1 − ε.
λ1,2 = λ3,4 a spetkr´aln´ı polomˇer je potom
ρ(TJ ) =
√
1 − ε.
Pro 0 < ε < 0, 1 je tedy splnˇena nutn´a a postaˇcuj´ıc´ı podm´ınka konvergence ρ(TJ ) < 1. Pro Gauss-Seidelovu metodu plat´ı: TGS = (D − L)−1 U, tj. TGS
=
1 −2(1 − ε) 0 (1 − ε)2 0 0 = 0 0
0 −1 0 0 0 0 −1 0 0 0 0 −1 0 0 0 0 −1 0 0 2(1 − ε) −1 0 . 0 0 −1 −(1 − ε)2 0 0 0 1 0 0
0 0 1 0
0 0 0 1
=
V´ ypoˇctem pomoc´ı syst´emu Matlab snadno zjist´ıme, ˇze spektr´aln´ı polomˇer ρ(TGS ) > 1. To znamen´a, ˇze nen´ı splnˇena nutn´a podm´ınka pro konvergenci. Pro p˚ uvodn´ı soustavu rovnic tedy Jacobiho metoda konverguje, ale Gauss-Seidelova metoda ne. Divergence Gauss-Seidelovy metody je vidˇet i z Obr´azku 2 v Pˇr´ıloze. 42
3.4
Superrelaxaˇ cn´ı metoda
K ˇreˇsen´ı soustavy Ax = b m˚ uˇzeme pouˇz´ıt obecnˇejˇs´ı iteraˇcn´ı metody z´avisej´ıc´ı na re´aln´em parametru ω, tzv. superrelaxaˇcn´ı metody (nˇekdy kr´atce oznaˇcujeme SOR metody) (D − ωL)xk+1 = ((1 − ω)D + ωU )xk + ωb, kde D je diagon´aln´ı matice, L doln´ı a U horn´ı troj´ uheln´ıkov´a matice a z´aroveˇ n plat´ı A = D − L − U. Vezmeme-li parametr ω = 1, dostaneme Gauss-Seidelovu metodu. Vˇ eta 3.4.1. Nutn´ a podm´ıka pro to, aby superrelaxaˇcn´ı metoda konvergovala pro kaˇzd´y poˇc´ateˇcn´ı vektor x0 a kaˇzdou pravou stranu b je, aby 0 < ω < 2. D˚ u k a z. Oznaˇcme λ1 , λ2 , . . ., λn vlastn´ı ˇc´ısla matice Tω = (D − ωL)−1 ((1 − ω)D + ωU ). Spetkr´aln´ı polomˇer Tω je roven maxi |λi | a, protoˇze matice D − ωL a (1 − ω)D + ωU jsou troj´ uheln´ıkov´e, jej´ı determinant det((D − ωL)−1 det((1 − ω)D + ωU ) = (detD)−1 (1 − ω)n detD = (1 − ω)n . Pak tedy plat´ı λ1 λ2 . . . λn = (1 − ω)n a ρn (Tω ) ≥ |λ1 ||λ2 | . . . |λn | = |1 − ω|n , z ˇcehoˇz vypl´ yv´a ρ(Tω ) ≥ |1 − ω|. Aby iteraˇcn´ı metoda konvergovala, mus´ı b´ yt ρ(Tω ) < 1 a tedy |1 − ω| < 1 neboli 0 < ω < 2.
Vˇ eta 3.4.2. Necht’ matice A je pozitivnˇe definitn´ı a A = D − L − U , U = L∗ , kde D je tak´e pozitivnˇe definitn´ı. Pak superrelaxaˇcn´ı metoda konverguje pro kaˇzd´e 0 < ω < 2, kaˇzdou poˇc´ateˇcn´ı aproximaci x0 a vektor b. D˚ u k a z. Necht’ 0 < ω < 2. Oznaˇcme L1 = (1 − ω)D + ωL,
D1 = (2 − ω)D. 43
Matice D1 je tedy pozitivnˇe definitn´ı a matice A1 = D1 − L1 − L∗1 = ω(D − L − L∗ ) = ωA je tak´e pozitivnˇe definitn´ı. Z tvrzen´ı [soustava Cx = y, C = D − L − L∗ pozitivnˇe definitn´ı, kde D je tak´e pozitivnˇe definitn´ı, pak matice (D−L) je regul´arn´ı a tzv. zobecnˇen´a Gauss-Seidelova metoda xk+1 = (D − L)−1 L∗ xk + (D − L)−1 y konverguje (viz [3])] plyne, ˇze ρ((D1 − L1 )−1 L∗1 ) < 1, ale (D1 − L1 )−1 L∗1 = (D − ωL)−1 ((1 − ω)D + ωL∗ ), takˇze SOR metoda konverguje. Zvol´ıme-li za m´ıru rychlosti konvergence spektr´aln´ı polomˇer pˇr´ısluˇsn´e matice, najdeme tak v nˇekter´ ych pˇr´ıpadech optim´aln´ı parametr SOR metody. Protoˇze plat´ı, ˇze rychlost konvergence metody je t´ım vˇetˇs´ı, ˇc´ım menˇs´ı je odpov´ıdaj´ıc´ı spektr´aln´ı polomˇer matice, dostaneme n´asleduj´ıc´ı nerovnost. Oznaˇcme ω0 optim´aln´ı parametr metody a pro ten plat´ı ρ((D − ω0 L)−1 ((1 − ω0 )D + ω0 U )) ≤ ρ (D − ωL)−1 ((1 − ωD) + ωU ) pro vˇsechna ω splˇ nuj´ıc´ı 0 < ω < 2. ˇ Nyn´ı zavedeme pojem konzistentn´ıho oˇc´ıslov´an´ı mnoˇziny uzl˚ u grafu. Rekneme, ˇze oˇc´ıslov´an´ı O(i) mnoˇziny uzl˚ u grafu G je konzistentn´ı, jestliˇze existuje dalˇs´ı oˇc´ıslov´an´ı W (i) mnoˇziny uzl˚ u takov´e, ˇze pro kaˇzdou hranu (i, k) plat´ı O(i) < O(k) ⇔ W (i) = W (k) − 1, O(i) > O(k) ⇔ W (i) = W (k) + 1. D´ale ˇrekneme, ˇze symetrick´a nebo hermitovsk´a matice m´a vlastnost X, pr´avˇe kdyˇz ji lze permutac´ı ˇr´adk˚ u a sloupc˚ u pˇrev´est na tvar ! D1 −U , −L D2 kde D1 a D2 jsou diagon´aln´ı matice ˇr´adu alespoˇ n 1. Nyn´ı jiˇz m˚ uˇzeme vyslovit d˚ uleˇzitou vˇetu. Oznaˇcen´ı: ρJ - spektr´aln´ı polomˇer Jacobiho iteraˇcn´ı metody, ρGS - spektr´aln´ı polomˇer Gauss-Seidelovy iteraˇcn´ı metody, ρω - spektr´aln´ı polomˇer SOR metody. Vˇ eta 3.4.3. Necht’ matice A soustavy Ax = b je pozitivnˇe definitn´ı a m´a vlastnost X. Pak kaˇzd´e konzistentn´ı oˇc´ıslov´an´ı uzl˚ u grafu G(A) odpov´ıd´a takov´emu poˇrad´ı rovnic soustavy, pˇri kter´em m´a Gauss-Seidelova metoda spektr´ aln´ı polomˇer ρGS = ρ2J , kde ρJ je spektr´ aln´ı polomˇer pˇr´ısluˇsn´e Jacobiho metody. Pˇritom plat´ı ρJ < 1. 44
D˚ u k a z. viz. [3] Obecnˇeji m´a pˇri stejn´em poˇrad´ı superrelaxaˇcn´ı metoda s parametrem ω spektr´aln´ı polomˇer q 1 ρω = {ωρJ + ω 2 ρ2J + 4(1 − ω)}2 , pro ω ≤ ω0 , 4 ρω = ω − 1, pro ω > ω0 , kde ω0 =
1+
√2
1−ρ2J
. Tento polomˇer je nejmenˇs´ı pro ω = ω0 , a pak ρω0
p 1 − ρ2J p = ω0 − 1 = . 1 + 1 − ρ2J 1−
Podle posledn´ı vˇety tedy spektr´aln´ı polomˇer Gauss-Seidelovy ani SOR metody nez´avis´ı na volbˇe konzistentn´ıho uspoˇr´ad´an´ı, jichˇz m˚ uˇze b´ yt nˇekolik. Vˇ eta 3.4.4. Necht’ matice A soustavy Ax = b je monot´onn´ı a takov´ a, ˇze A = M −N , kde M je monot´onn´ı a N ≥ 0. Pak iteraˇcn´ı metoda M xk+1 = N xk + b
(3.1)
konverguje pro kaˇzd´y poˇc´ateˇcn´ı vektor x0 a kaˇzdou pravou stranu b. Matice A se naz´ yv´a monot´onn´ı, jestliˇze plat´ı A−1 ≥ 0 a pro kaˇzd´e dva vektory, splˇ nuj´ıc´ı Ax ≥ Ay, plat´ı x ≥ y. D˚ u k a z. (vˇety) Protoˇze (3.1) je ekvivalentn´ı iteraci xk+1 = M −1 N xk + M −1 b, rozhoduje o konvergenci spektr´aln´ı polomˇer ρ(M −1 N ). Protoˇze M −1 ≥ 0, je M −1 N ≥ 0. Pak existuje vektor z 6= 0, z ≥ 0, tak, ˇze M −1 N z = ρz. D´ale je A−1 N = A−1 M M −1 N = A−1 (A + N )M −1 N = (I + A−1 N )M −1 N, odkud A−1 N z = (I + A−1 N )M −1 N z = (I + A−1 N )ρz = ρz + A−1 N ρz. Protoˇze A−1 ≥ 0 a A−1 N ≥ 0, pro ρ ≥ 1 by platilo (ρ − 1)A−1 N z + ρz = 0 a to nen´ı moˇzn´e. Proto je ρ < 1. Nˇekdy lze pouˇz´ıt i tzv. symetrickou superrelaxaˇcn´ı metodu. V t´eto metodˇe je jeden krok vlastnˇe sloˇzen´ y ze dvou ”p˚ ulkrok˚ u”. Opˇet vych´az´ıme ze z´apisu matice A = D − L − U . Pak 1
(D − ωL)xk+ 2 = ((1 − ω)D + ωU )xk + ωb, 45
1
(D − ωU )xk+1 = ((1 − ω)D + ωL)xk+ 2 + ωb. Matice metody, pomoc´ı kter´e poˇc´ıt´ame xk+1 z xk je tvaru (D − ωU )−1 ((1 − ω)D + ωL)(D − ωL)−1 ((1 − ω)D + ωU ). O konvergenci a jej´ı rychlosti pak opˇet rozhoduje spektr´aln´ı polomˇer t´eto matice. Pˇ r´ıklad 3.4.1. Zjistˇete spektr´ aln´ı polomˇer superrelaxaˇcn´ı metody uˇzit´e na soustavu rovnic 3x1 + 2x2 + 2x3 = 1 2x1 + 3x2 + 2x3 = 0 2x1 + 2x2 + 3x3 = −1 pro parametr a) ω = 1, b) ω = 1, 1, c) ω = 0, 9, d) ω = 0, 926. ˇ sen´ı: Nejdˇr´ıve Reˇ 3 0 D= 0 3 0 0
urˇc´ıme matice D, L a U 0 0 0 0 , L = −2 0 3 −2 −2
tak, aby A = D − L − U . Tedy 0 0 −2 −2 0 , U = 0 0 −2 . 0 0 0 0
Pak tyto matice dosad´ıme do matice superrelaxaˇcn´ı metody (D − ωL)−1 ((1 − ω)L + ωU ),
urˇc´ıme vlastn´ı ˇc´ısla a z nich spektr´aln´ı polomˇer. Dalˇs´ı v´ ypoˇcty jsou provedeny ve v´ ypoˇcetn´ım syst´emu Matlab. a) Pro ω = 1 bude matice metody 0 −0, 6667 −0, 6667 0 0, 4444 −0, 2222 , 0 0, 1481 0, 5926
jej´ı vlastn´ı ˇc´ısla λ1 = 0, λ2 = 0, 5185 + 0, 1656i, λ3 = 0, 5185 − 0, 1656i. Odtud je spektr´aln´ı polomˇer ρ = 0, 5443. b) Pro ω = 1, 1 bude matice metody −0, 1 −0, 7333 −0, 7333 0, 0733 0, 4378 −0, 1956 , 0, 0196 0, 2167 0, 5812
jej´ı vlastn´ı ˇc´ısla λ1 = −0, 0038, λ2 = 0, 4614 + 0, 2316i, λ3 = 0, 4614 − 0, 2316i. Odtud je spektr´aln´ı polomˇer ρ = 0, 5162. c) Pro ω = 0, 9 bude matice metody 0, 1 −0, 6 −0, 6 −0, 06 0, 46 −0, 24 , −0, 024 0, 084 0, 604 46
jej´ı vlastn´ı ˇc´ısla λ1 = 0, 0029, λ2 = 0, 5806 + 0, 1167i, λ3 = 0, 5806 − 0, 1167i. Odtud je spektr´aln´ı polomˇer ρ = 0, 5922. d) Pro ω = 0, 926 bude matice metody 0, 074 −0, 6173 −0, 6173 −0, 0457 0, 4551 −0, 2362 , −0, 0175 0, 1002 0, 6009
jej´ı vlastn´ı ˇc´ısla λ1 = 0, 0012, λ2 = 0, 5644 + 0, 1279i, λ3 = 0, 5644 − 0, 1279i. Odtud je spektr´aln´ı polomˇer ρ = 0, 5787. Konvergenci aproximac´ı pro jednotliv´e pˇr´ıpady najdeme na Obr. 3 - Obr. 6 v Pˇr´ıloze.
47
Kapitola 4 Singul´ arn´ı line´ arn´ı syst´ emy 4.1
Semikonvergentn´ı matice
Pro ˇreˇsen´ı singul´arn´ıch syst´em˚ u line´arn´ıch rovnic se pouˇz´ıv´a tzv. semiiteraˇcn´ı metoda. K zaveden´ı konvergenˇcn´ıch krit´eri´ı t´eto metody budeme potˇrebovat definovat pojem semikonvergentn´ı matice. D´ale uvedeme nˇekolik vˇet t´ ykaj´ıc´ıch se tˇechto matic. ˇ Definice 4.1.1. Ctvercov´ a matice H se naz´yv´a semikonvergentn´ı, kdyˇz lim H j
j→∞
existuje. Pˇripomeˇ nme, ˇze H je konvergentn´ı tehdy a jen tehdy, je-li ρ(H) < 1. Pro semikonvergentn´ı matice plat´ı obdobn´e tvrzen´ı. Vˇ eta 4.1.1. Necht’ H je ˇctvercov´a matice. Pak H je semikonvergentn´ı pr´avˇe tehdy, kdyˇz plat´ı n´ asleduj´ıc´ı podm´ınky: 1. ρ(H) ≤ 1, 2. ρ(H) = 1, pak element´arn´ı dˇelitel´e asociovan´ı s vlastn´ım ˇc´ıslem 1 matice H jsou line´arn´ı, tj. h(I − H)2 = h(I − H), 3. ρ(H) = 1, pak λ ∈ σ(H), |λ| = 1 implikuje λ = 1. D˚ u k a z. viz [2] Definice 4.1.2. Necht’ A je ˇctvercov´a matice s indexem k. Drazinova inverze matice A je matice AD , kter´ a splˇ nuje podm´ınky AD AAD = AD ,
AAD = AD A, 48
Ak = AD Ak+1 .
Pro Drazinovu inverzi tak´e plat´ı ( y, Ay = x, x ∈ R(Ak ) AD = 0, Ak = 0.
(4.1)
Lemma 4.1.1. Necht’ H ∈ Cn,n . Pak H je semikonvergentn´ı pr´avˇe tehdy, kdyˇz existuje regul´arn´ı matice P tak, ˇze ! I 0 H=P P −1 , 0 K kde I chyb´ı, kdyˇz 1 6∈ σ(H) a kde ρ(K) < 1 nebo K chyb´ı. D˚ u k a z. viz [2] Lemma 4.1.2. Necht’ H je semikonvergentn´ı. Pak lim H i = I − E,
i→∞
kde E = (I − H)(I − H)D . D˚ u k a z. Plyne pˇr´ımo z Lemmatu 4.1.1 a z charakterizace Drazinovy inverze (4.1).
Nyn´ı jeˇstˇe uvedeme nˇekolik doplˇ nuj´ıc´ıch tvrzen´ı o semikonvergentn´ıch matic´ıch. ˇ adkovˇe stochastick´a matice je semikonvergentn´ı pr´avˇe tehdy, kdyˇz Lemma 4.1.3. R´ |λ| = 1 ⇒ λ = 1. D˚ u k a z. viz [2] D˚ usledek 4.1.1. Necht’ S je ˇr´adkovˇe stochastick´a matice, pro kterou plat´ı 0 < sii < 1,
i = 1, . . . , n.
Pak S je semikonvergentn´ı. D˚ u k a z. Plyne z tvrzen´ı, ˇze je-li S ˇr´adkovˇe stochastick´a matice, pro kterou plat´ı 0 < sii < 1, i = 1, . . . , n, pak |λ| = 1 ⇒ λ = 1, tj. λ = 1 je vlastn´ı hodnota s absolutn´ı hodnotou 1. D˚ usledek 4.1.2. Necht’ S1 , S2 jsou ˇr´adkovˇe stochastick´e matice splˇ nuj´ıc´ı pˇredpoklad pˇredchoz´ıho D˚ usledku. Pak matice ! 0 −S1 T = 0 S2 S1 je semikonvergentn´ı. D˚ u k a z. Zˇrejm´ y. Jeˇstˇe budeme potˇrebovat definici indexu matice. Definice 4.1.3. Index ˇctvercov´e matice A je nejmenˇs´ı nez´ aporn´e ˇc´ıslo k takov´e, ˇze h(Ak ) = h(Ak+1 ). 49
4.2
Semiiteraˇ cn´ı metoda
Nyn´ı jiˇz m˚ uˇzeme rozˇs´ıˇrit konvergenci iteraˇcn´ıch metod na singul´arn´ı pˇr´ıpad. Uvaˇzujme tedy konzistentn´ı line´arn´ı syst´em Ax = b se singul´arn´ı matic´ı A, kter´ y je vˇsak ˇreˇsiteln´ y. Vˇ eta 4.2.1. Necht’ A = M − N ∈ Cn,n s M regul´arn´ı. Pak pro H = M −1 N a c = M −1 b iteraˇcn´ı metoda xk+1 = Hxk + c (4.2) konverguje k ˇreˇsen´ı x∗ soustavy Ax = b, pro kaˇzd´e x0 pr´avˇe tehdy, kdyˇz H je semikonvergentn´ı. V tomto pˇr´ıpadˇe lim xk = (I − H)D c + (I − E)x0 ,
k→∞
(4.3)
kde E = (I − H)(I − H)D . D˚ u k a z. Pˇredpokl´adejme, ˇze x∗ je ˇreˇsen´ı soustavy Ax = b. Pak x∗ = Hx∗ + c. Nyn´ı odeˇcteme-li tuto rovnici od (4.2), dostaneme chybu rovnice xk+1 − x∗ = H(xk − x∗ ) = . . . = H k+1 (x0 − x∗ ). Odtud posloupnost {xi } konverguje k nˇejak´emu ˇreˇsen´ı syst´emu Ax = b pr´avˇe tehdy, kdyˇz posloupnost {xi − x∗ } konverguje k vektoru nulov´e vzd´alenosti od A. Ale to je pr´avˇe kdyˇz lim H i (x0 − x∗ ) (4.4) i→∞
existuje. Pak z Lemmatu 4.1.1 plyne, ˇze (4.4) plat´ı pr´avˇe, kdyˇz H je semikonvergentn´ı. Nav´ıc, je-li H semikonvergentn´ı, pak identita (4.3) plyne z Lemmatu 4.1.2. Oznaˇcme r vektor rezidu´ı. Asociovan´e norm´aln´ı rovnice syst´emu Ax = b jsou tvaru r + Ax = b, kde AT r = o,
nebo − AT r = o.
Pak iteraˇcn´ı proces rk+1 = −Axk + b
xk+1 = AT rk+1 + xk je zˇrejmˇe ekvivalentn´ı xk+1 = Hxk + c,
H = I − AT A, 50
c = AT b.
V´ yhoda je v tom, ˇze se vyhneme poˇc´ıt´an´ı AT A. D˚ uleˇzit´ ym n´astrojem pro vyˇsetˇrov´an´ı konvergence iteraˇcn´ıch metod pro singul´arn´ı syst´emy je pouˇzit´ı kvadratick´ ych forem a pozitivn´ı definitnosti. D˚ ukazy konvergenˇcn´ıch v´ ysledk˚ u pro singul´arn´ı pˇr´ıpad jsou komplikovanˇejˇs´ı, ale obdobn´e jako pro regul´arn´ı pˇr´ıpad. Nejsou zde proto uvedeny. Pˇripomeˇ nme jeˇstˇe, ˇze A∗ znaˇc´ı matici transponovanou a komplexnˇe sdruˇzenou k matici A. Vˇ eta 4.2.2. Necht’ A ∈ Cn,n je hermitovsk´ a (tj. A∗ = A). Necht’ A = M − N s M regul´arn´ı a necht’ H = M −1 N . Pˇredpokl´adejme, ˇze M ∗ + N je pozitivnˇe definitn´ı. Pak H je semikonvergentn´ı pr´avˇe tehdy, kdyˇz A je pozitivnˇe semidefinitn´ı. Nakonec jiˇz m˚ uˇzeme na singul´arn´ı pˇr´ıpad rozˇs´ıˇrit konkr´etn´ı iteraˇcn´ı metody. D˚ usledek 4.2.1. Necht’ A ∈ Cn,n m˚ uˇzeme rozloˇzit na A=D−L−U a pˇredpokl´adejme, ˇze A je hermitovsk´ a a D pozitivnˇe definitn´ı. Potom blokov´a Jacobiho iteraˇcn´ı matice TJ je semikonvergentn´ı pr´avˇe tehdy, kdyˇz 2D − A a A jsou pozitivnˇe definitn´ı matice. D˚ usledek 4.2.2. Mˇejme stejn´e pˇredpoklady jako v pˇredchoz´ım D˚ usledku. Pak blokov´ a superrelaxaˇcn´ı iteraˇcn´ı matice (D − ωL)−1 ((1 − ω)D + ωU ) je semikonvergentn´ı pro vˇsechna 0 < ω < 2 pr´avˇe tehdy, kdyˇz A je pozitivnˇe definitn´ı. Pˇ r´ıklad 4.2.1. Semiiteraˇcn´ı metodou ˇreˇste syst´em x − 2y = 1
2x − 4y = 2. ˇ sen´ı: Matice soustavy Reˇ 1 −2 2 −4
A=
!
je singul´arn´ı, ale soustava je ˇreˇsiteln´a. Pro ˇreˇsen´ı t´eto soustavy semiiteraˇcn´ı metodou potˇrebujeme rozloˇzit matice A na matici N a regul´arn´ı matici M tak, ˇze A = M −N . Pak tedy ! ! 0 6 −1 8 M= a N= . 8 0 6 4 Nyn´ı T =M
−1
N=
g=M
−1
0 1 6
b=
! −1 8 = 0 6 4 ! ! 1 0 18 = 1 0 2 6 1 8
!
51
3 4
− 16 !
1 4 1 6
1 2 4 3
.
!
,
Iteraˇcn´ı proces pak vypad´a takto: ! 3 xk+1 4 = k+1 y − 61 xk+1 y k+1
!
1 2 4 3
!
xk yk
!
+
3 k x + 12 y k + 41 4 − 16 xk + 43 y k + 61
=
1 4 1 6
!
.
x3 y3
!
!
Jako poˇc´ateˇcn´ı aproximaci zvol´ıme x0 y0 a potom dost´av´ame: ! x1 = y1
2 8 3
!
,
x2 y2
!
!
=
=
1 2
37 12 61 18
! !
,
=
613 144 901 216
!
Vˇsechny aproximace leˇz´ı na jedn´e pˇr´ımce, jak ukazuje Obr. 7 v Pˇr´ıloze. Pˇr´ımka 1 1 y = x− 2 2 je pˇr´ımka zadan´a soustavou a aproximace leˇz´ı na pˇr´ımce 4 2 y = x+ . 3 3
52
.
Z´ avˇ er C´ılem t´eto pr´ace bylo popsat r˚ uzn´e typy speci´aln´ıch matic a jejich vlastnosti a zab´ yvat se ˇreˇsen´ım soustav line´arn´ıch rovnic s tˇemito maticemi. Tzv. iteraˇcn´ı metody jsou d˚ uleˇzit´ ym prostˇredkem ˇreˇsen´ı takov´ ych soustav. Metody uveden´e v t´eto pr´aci lze vyuˇz´ıt k ˇreˇsen´ı probl´em˚ u napˇr. matematick´e fyziky nebo ekonomie, kter´e ˇcasto vedou na soustavy vysok´ ych ˇr´ad˚ u. V souˇcasn´e dobˇe se v´ yzkum vˇenuje i semiiteraˇcn´ı metodˇe, jej´ıˇz z´aklady jsou zde uvedeny. Tyto probl´emy je bez vyuˇzit´ı poˇc´ıtaˇce t´emˇeˇr nemoˇzn´e ˇreˇsit. Programy vytvoˇren´e v t´eto diplomov´e pr´aci mohou b´ yt tedy pouˇzity v praxi pro jednoduˇsˇs´ı v´ ypoˇcet. Douf´am, ˇze tato pr´ace pom˚ uˇze dalˇs´ım student˚ um ve studiu soustav se speci´aln´ımi maticemi.
53
Literetura ˇ [1] BASTINEC J., DIBL´IK J.: Matematika IV. VUT v Brnˇe, Brno 1991. [2] BERMAN A., PLEMMONS R. J.: Nonnegative Matrices in the Mathematical Sciences. SIAM, Philadelphia 1994. [3] FIEDLER M.: Speci´aln´ı matice a jejich pouˇzit´ı v numerick´e matematice. SNTL, Praha 1981. ´ I.: Numerick´e metody. MU v Brnˇe, Brno 1999. [4] HOROVA ´ I., BUD´IKOVA ´ M.: Linear smoothers and additive models. MU [5] HOROVA v Brnˇe, Brno 1998. ˇ [6] MARCUK G. I.: Metody numerick´e matematiky. Academia, Praha 1987. [7] M´IKA S.: Numerick´e metody algebry. SNTL, Praha 1982. ´ ˇ sen´ı syst´em˚ [8] NAVRATIL L.: Reˇ u line´arn´ıch rovnic semiiteraˇcn´ı metodou. Diplomov´a pr´ace, Brno 1998. [9] RALSTON A.: Z´aklady numerick´e matematiky. Academia, Praha 1973. ˇ F.: Line´arn´ı algebra zamˇeˇren´a na numerickou anal´ [10] SIK yzu. MU v Brnˇe, Brno 1998.
54
Pˇ r´ıloha V pˇr´ıloze jsou uvedeny nˇekter´e programy pro v´ ypoˇcetn´ı syst´em Matlab. Prvn´ı poˇc´ıt´a ˇreˇsen´ı syst´emu line´arn´ıch rovnic Gauss-Seidelovou metodou a druh´ y metodou superrelaxaˇcn´ı pro r˚ uzn´e parametry ω. Jako vstupn´ı data tedy potˇrebujeme ˇctvercovou matici soustavy A, vektor b, vektor poˇc´ateˇcn´ıch aproximac´ı x0 a poˇcet iteraˇcn´ıch krok˚ u, kter´e se maj´ı poˇc´ıtat, tj. n. Pro SOR metodu jeˇstˇe mus´ıme zadat, pro kter´e ω m´a b´ yt v´ ypoˇcet proveden. Po proveden´ı v´ ypoˇctu z´ısk´ame n-tou aproximaci ˇreˇsen´ı a obr´azek. Z toho snadno zjist´ıme, zda m´ame jiˇz dostateˇcnˇe pˇresn´e ˇreˇsen´ı nebo je nutn´e zv´ yˇsit poˇcet iterac´ı. V dalˇs´ı ˇc´asti pˇr´ılohy jsou obr´azky z´ıskan´e z tˇechto program˚ u pro nˇekter´e pˇr´ıklady uvedn´e v jednotliv´ ych kapitol´ach. Vˇsechny obr´azky jsou vykresleny pro z ∈ (1, n).
55
Gauss-Seidelova metoda function [x] = gs(A, b, x0, n) s = size(A); L = −tril(A); U = −triu(A); for i = 1 : s(1); U (i, i) = 0; L(i, i) = 0; end D = A + L + U; T = inv(D − L) ∗ U ; c = inv(D − L) ∗ b; P = zeros(n, s(1)); P (1, :) = x0′ ; for i = 2 : n; P (i, :) = (T ∗ P (i − 1, :)′ + c)′ ; end; z = 1 : n; for i = 1 : s(1); plot(z, P (:, i)′ ); hold on; end; hold off; x = P (n, :)′
56
Superrelaxaˇ cn´ı metoda function [x] = sor(A, b, omg, x0, n) s = size(A); L = −tril(A); U = −triu(A); for i = 1 : s(1); U (i, i) = 0; L(i, i) = 0; end D = A + L + U; T = inv(D − omg ∗ L) ∗ ((1 − omg) ∗ D + omg ∗ U ); c = inv(D − omg ∗ L) ∗ (omg ∗ b); P = zeros(n, s(1)); P (1, :) = x0′ ; for i = 2 : n; P (i, :) = (T ∗ P (i − 1, :)′ + c)′ ; end; z = 1 : n; for i = 1 : s(1); plot(z, P (:, i)′ ); hold on; end; hold off; x = P (n, :)′
57
