Numerické metody pro nalezení

Masarykova univerzita Brno Fakulta pˇr´ırodovˇedeck´a Katedra aplikovan´e matematiky

Numerick´ e metody pro nalezen´ı vlastn´ıch ˇ c´ısel matic Diplomov´a pr´ace

kvˇeten 2006

Alena Baˇstincov´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 20. kvˇetna 2006

Obsah 1 Z´ akladn´ı kapitola

7

2 Typy metod pro hled´ an´ı vlastn´ıch ˇ c´ısel

8

3 Klasick´ e metody urˇ cen´ı koeficient˚ u charakteristick´ eho polynomu 3.1 Krylovova metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Faddˇejevova-Leverrierova metoda . . . . . . . . . . . . . . . . . . . . . . . .

10 10 11

4 Poloha a odhad vlastn´ıch ˇ c´ısel 4.1 Gerˇsgorinovy vˇety . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13 13

5 Metody v´ ypoˇ ctu dominantn´ıho 5.1 Mocninn´a metoda . . . . . . . 5.2 Metoda Rayleighova pod´ılu . 5.3 V´ ypoˇcet dalˇs´ıch vlastn´ıch ˇc´ısel

17 17 20 22

vlastn´ıho ˇ c´ısla . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . mocninnou metodou . . . . . . . . . . . . . .

6 Metody pro v´ ypoˇ cet vlastn´ıch ˇ c´ısel a vlastn´ıch vektor˚ u symetrick´ ych matic 6.1 Jacobiho metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Householderova matice zrcadlen´ı . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Givensova-Householderova metoda . . . . . . . . . . . . . . . . . . . . . . . 6.3.1 Householderova metoda . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.2 Givensova metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 QR-rozklad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Konstrukce QR-rozkladu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5.1 QR-rozklad pomoc´ı Gram-Schmidtova algoritmu . . . . . . . . . . . . 6.5.2 QR-rozklad pomoc´ı Householderovy matice . . . . . . . . . . . . . . 6.5.3 QR-rozklad pomoc´ı Givensovy matice . . . . . . . . . . . . . . . . . . 6.5.4 Srovn´an´ı algoritm˚ u . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5.5 QR-rozklad a vlastn´ı ˇc´ısla matice A – QR-algoritmus . . . . . . . . .

24 24 32 35 35 37 40 40 40 42 46 48 49

7 Podm´ınˇ enost probl´ emu vlastn´ıch ˇ c´ısel 7.1 Glob´aln´ı ˇc´ıslo podm´ınˇenosti . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Odhad chyby vypoˇc´ıtan´eho vlastn´ıho ˇc´ısla . . . . . . . . . . . . . . . . . . . 7.3 Relativn´ı chyba vypoˇc´ıtan´eho vlastn´ıho ˇc´ısla . . . . . . . . . . . . . . . . . .

51 51 52 53

4

´ Uvod C´ılem m´e diplomov´e pr´ace je popsat numerick´e metody pro nalezen´ı vlastn´ıch ˇc´ısel matic. Vlastn´ı ˇc´ısla a vlastn´ı vektory maj´ı velmi ˇsirok´e spektrum aplikac´ı, napˇr´ıklad se pouˇz´ıvaj´ı pˇri hled´an´ı ˇreˇsen´ı diferenci´aln´ıch rovnic a jejich soustav a to jak u obyˇcejn´ ych difarenci´aln´ıch rovnic, tak u parci´aln´ıch diferenci´aln´ıch rovnic a jejich soustav. Tot´eˇz plat´ı i pro diferenˇcn´ı rovnice a jejich soustavy. Mnoh´e technick´e probl´emy se daj´ı popsat pomoc´ı diferenci´aln´ıch nebo diferenˇcn´ıch rovnic a jejich soustav, jako napˇr´ıklad popis obvod˚ u v elektrotechnice.Pokud m´a obvod vˇetˇs´ı poˇcet prvk˚ u, dost´av´ame soustavu diferenci´aln´ıch rovnic vyˇsˇs´ıho ˇr´adu. Pro jejich ˇreˇsen´ı potˇrebujeme zn´at vlastn´ı ˇc´ısla matice soustavy. Odtud je zˇrejm´a duleˇzitost u ´lohy o nalezen´ı vlastn´ıch ˇc´ısel matice. Pˇr´ım´e metody hled´an´ı vlastn´ıch ˇc´ısel jsou mnohdy neefektivn´ı a proto je nutn´e ˇreˇsit tuto u ´lohu numericky. Pˇri numerick´em ˇreˇsen´ı se sice dopouˇst´ıme urˇcit´e chyby, ale souˇcasnˇe se dostaneme k ˇreˇsen´ı, alespoˇ n pˇribliˇzn´emu, v relativnˇe kratˇs´ım ˇcase s poˇzadovanou pˇresnost´ı. Ve sv´e pr´aci nejdˇr´ıve definuji z´akladn´ı pojmy, nevˇenuji se pˇr´ım´ ym metod´am v´ ypoˇctu vlastn´ıch ˇc´ısel a zab´ yv´am se numerick´ ymi metodami jejich urˇcen´ı. Postupnˇe uv´ad´ım ˇradu zp˚ usob˚ u nalezen´ı vlastn´ıch ˇc´ısel a jim pˇr´ısluˇsn´ ym vlastn´ım vektor˚ um. Nejdˇr´ıve uv´ad´ım klasick´e metody urˇcen´ı koˇren˚ u charakteristick´eho polynomu, d´ale se vˇenuji odhadu polohy vlastn´ıch ˇc´ısel, pot´e n´asleduj´ı metody v´ ypoˇctu dominantn´ıho vlastn´ıho ˇc´ısla. Nejv´ıce m´ısta vˇenuji metod´am pro v´ ypoˇcet vlastn´ıch ˇc´ısel symetrick´ ych matic. Z´avˇereˇcn´a kapitola je vˇenov´ana probl´emu podm´ınˇenosti vlastn´ıch ˇc´ısel. Nevˇenovala jsem se rozboru jednotliv´ ych algoritm˚ u pˇri jejich zpracov´an´ı na poˇc´ıtaˇci, protoˇze tato problematika z´avis´ı na volbˇe programovac´ıho jazyka a softwarov´em vybaven´ı poˇc´ıtaˇce.

5

Oznaˇ cen´ı N Z R C pn (x) Am,n A = (aij ) I ei O o det A = |A| A−1 hod (A) tr(A) ρ(A) AH (Rn , +, .) dim P < ·, · > kxk kxk2 kAk kAk2 kAk∞ 

mnoˇzina pˇrirozen´ ych ˇc´ısel mnoˇzina cel´ ych ˇc´ısel mnoˇzina re´aln´ ych ˇc´ısel mnoˇzina komplexn´ıch ˇc´ısel polynom n-t´eho stupnˇe promˇenn´e x matice typu m, n (s m ˇr´adky a n sloupci) matice s prvky aij jednotkov´a matice jednotkov´ y vektor s 1 na i-t´em m´ıstˇe nulov´a matice nulov´ y vektor determinant matice A matice inverzn´ı k matici A hodnost matice A stopa matice A spektr´aln´ı polomˇer matice A ¯T matice hermitovsky sdruˇzen´a, tj.AH = A vektorov´ y prostor vˇsech uspˇr´adan´ ych n-tic dimenze prostoru P . standardn´ı skal´arn´ı souˇcin norma vektoru x eukleidovsk´a norma vektoru x norma matice A euklidovsk´a norma matice A krychlov´a norma matice A konec d˚ ukazu

6

Kapitola 1 Z´ akladn´ı kapitola Definice 1.0.1. Necht’ A je ˇctvercov´a matice ˇr´adu n. Jej´ı vlastn´ı ˇc´ısla λ1 , . . . , λn jsou koˇreny rovnice det(A − λI) = 0,

zvan´e charakteristick´a rovnice. Ke kaˇzd´emu vlastn´ımu ˇc´ıslu λi existuje aspoˇ n jedno nenulov´e (1) (2) (n) T ˇreˇsen´ı soustavy rovnic Ax = λi x. Toto ˇreˇsen´ı xi , kde xi = (xi , xi , . . . , xi ), nazveme prav´ym vlastn´ım vektorem matice A. (Vˇsude v dalˇs´ım bude pojem vlastn´ı vektor znaˇcit v´ yhradnˇe prav´ y vlastn´ı vektor.) Lev´y vlastn´ı vektor yi odpov´ıdaj´ıc´ı vlastn´ımu ˇc´ıslu λi je ˇreˇsen´ım rovnice yT A = λi yT . Lev´ y vlastn´ı vektor matice A je tedy vlastn´ım vektorem T transponovan´e matice A a snadno lze uk´azat , ˇze odpov´ıd´a-li lev´ y vlastn´ı vektor yk vlastn´ımu ˇc´ıslu λk a prav´ y vlastn´ı vektor xi vlastn´ımu ˇc´ıslu λi a plat´ı λk 6= λi jsou vektory yk a xi ortogon´aln´ı. (Ve vˇetˇsinˇe d´ale uveden´ ych pˇr´ıklad˚ u se budou vyskytovat re´aln´e matice , budeme pˇredpokl´adat , pokud nebude ˇreˇceno jinak, ˇze matice A je re´aln´a. Mnoh´e vˇety budou vˇsak platit i pro komplexn´ı matice nebo budeme-li pˇredpokl´adat symetrii, pro hermitovsk´e matice, (d˚ ukazy n´asleduj´ıc´ıch vˇet viz. [8]). Vˇ eta 1.0.1. Jsou-li λ1 , . . . , λn vlastn´ı ˇc´ısla matice A, m´a matice Ak vlastn´ı ˇc´ısla λk1 , . . . , λkn . Obecnˇeji, je-li p(x) libovoln´y polynom, m´a matice p(A) vlastn´ı ˇc´ısla p(λ1 ), . . . , p(λn ). Vˇ eta 1.0.2. Je-li matice A re´aln´a a symetrick´a, jsou vˇsechna jej´ı vlastn´ı ˇc´ısla a vˇsechny pˇr´ısluˇsn´e vlastn´ı vektory re´aln´e. Kromˇe toho vlastn´ı vektory pˇr´ısluˇsn´e r˚ uzn´ym vlastn´ım ˇc´ısl˚ um jsou ortogon´aln´ı a lev´y vlastn´ı vektor a prav´y vlastn´ı vektor pˇr´ısluˇsn´e t´emuˇz vlastn´ımu ˇc´ıslu jsou si rovny. Vˇ eta 1.0.3. Podobnostn´ı transformace PAP−1 nemˇen´ı vlastn´ı ˇc´ısla matice A. Vˇ eta 1.0.4. (Cayley-Hamilton) Necht’ je f (λ) = det(A−λI) = 0 charakteristick´a rovnice matice A. Pak plat´ı f (A) = 0. Vˇ eta 1.0.5. Vlastn´ı ˇc´ısla horn´ı (doln´ı) troj´ uheln´ıkov´e matice jsou prvky na jej´ı diagon´ale. Vˇ eta 1.0.6. Libovoln´a matice A je podobn´a diagon´ aln´ı matici D pr´avˇe tehdy, kdyˇz m´a matice Akompletn´ı soubor n line´arnˇe nez´ avisl´ych vlastn´ıch vektor˚ u. 7

Kapitola 2 Typy metod pro hled´ an´ı vlastn´ıch ˇ c´ısel Podle z´akladn´ı definice v´ıme, ˇze vlastn´ı ˇc´ısla dan´e matice jsou koˇreny jej´ıho charakteristick´eho polynomu. Z algebraick´e teorie uˇzeme algebraicky √ v´ıme, ˇze koˇreny polynomu stupnˇe n > 4 nem˚ (tj. pomoc´ı operac´ı ±, ×, ÷, ) vyj´adˇrit ve tvaru vzorce. Proto se obecnˇe nedaj´ı z´ıskat vlastn´ı ˇc´ısla pˇresnˇe ( aˇz na zaokrouhlovac´ı chyby) po koneˇcn´em poˇctu operac´ı. K ˇreˇsen´ı naˇseho probl´emu m˚ uˇzeme pˇristupovat v´ıce zp˚ usoby. 1. Pouˇzijeme-li libovolnou metodu na hled´an´ı koˇren˚ u charakteristick´eho polynomu p(λ). Pro jednoduch´ y koˇren m˚ uˇzeme pouˇz´ıt Newtonovu metodu ci+1 = ci − p(ci )/p′ (ci )

i = 1, 2, . . . ,

pˇri vhodn´e volbˇe poˇc´ateˇcn´ı aproximace c0 , metodu seˇcen, metodu p˚ ulen´ı intervalu atd. Modifikovan´a Newtonova metoda se d´a pouˇz´ıt i na hled´an´ı n´asobn´ ych koˇren˚ u. V pˇr´ıpadˇe komplexnˇe sdruˇzen´e dvojice koˇren˚ u m˚ uˇzeme pouˇz´ıt napˇr. Bairstowovu metodu. Hled´an´ı velk´eho poˇctu koˇren˚ u t´ımto zp˚ usobem je vˇsak dost n´aroˇcn´e a probl´em b´ yv´a nestabiln´ı. 2. Z´ısk´an´ı vlastn´ıch ˇc´ısel bez znalost´ı charakteristick´eho polynomu, pˇri vyuˇz´ıv´an´ı vlastnost´ı podobn´ ych matic. C´ılem je naj´ıt podobnou matici v jednoduˇsˇs´ım tvaru, ze kter´eho se d´a vlastn´ı ˇc´ıslo urˇcit (napˇr´ıklad z diagon´aln´ı nebo troj´ uheln´ıkov´e matice). Takovou matici (nˇekdy jen nˇekter´e jej´ı vlastn´ı ˇc´ıslo) m˚ uˇzeme z´ıskat jako limitu posloupnosti podobnostn´ıch transformac´ı. V´ ybˇer tˇechto transformac´ı b´ yv´a zaloˇzen na speci´aln´ıch vlastnostech matic a jejich vlastn´ıch vektor˚ u. 3. Neline´arn´ı pˇr´ıstup, vlastn´ı probl´em (A − λI)x = 0 uvaˇzujeme jako soustavu n rovnic pro ych x1 , ..., xn , λ, kterou dopln´ıme P n 2+ 1 nezn´am´ normovanou podm´ınkou napˇr´ıklad xi = 1 na soustavu n + 1 neline´arn´ıch rovnic. Tato soustava se d´a ˇreˇsit napˇr´ıklad Newtonovou metodou. Pˇritom se vˇsak nevyuˇz´ıvaj´ı algebraick´e vlastnosti soustavy, kter´e m˚ uˇzou v´ ypoˇcet znaˇcnˇe ulehˇcit. Proto je tento postup znaˇcnˇe neefektivn´ı. Pozn´ amka 2.0.1. Pod pojmem u ´pln´y probl´em vlastn´ıch ˇc´ısel se rozum´ı u ´loha naj´ıt vˇsechna vlastn´ı ˇc´ısla a pˇr´ıpadnˇe i pˇr´ısluˇsn´e vlastn´ı vektory. 8

Pojem ˇc´asteˇcn´y probl´em vlastn´ıch ˇc´ısel znamen´a naj´ıt jedno nebo v´ıce vlastn´ıch ˇc´ısel spolu s pˇr´ısluˇsn´ ymi vlastn´ımi vektory. ´ y a ˇc´asteˇcn´ Upln´ y probl´em vystupuj´ı jako naprosto odliˇsn´e u ´lohy nejen oborem aplikac´ı, ale i metodami ˇreˇsen´ı. ˇ sen´ı u Reˇ ´pln´eho probl´emu je n´aroˇcnˇejˇs´ı. Neexistuje univerz´aln´ı algoritmus, kter´ y by byl stejnˇe efektivn´ı pro vˇsechny typy matic.

9

Kapitola 3 Klasick´ e metody urˇ cen´ı koeficient˚ u charakteristick´ eho polynomu Dˇr´ıve se vˇetˇsina metod na v´ ypoˇcet vlastn´ıch ˇc´ısel zakl´adala pr´avˇe na v´ ypoˇctu koeficient˚ u charakteristick´eho polynomu. Jejich v´ ypoˇcet pomoc´ı souˇctu hlavn´ıch minor˚ u je vˇsak nerentabiln´ı. Existuj´ı mnohem jednoduˇsˇs´ı metody na urˇcen´ı koeficient˚ u, kter´e maj´ı stejn´ y charakter (tj. pˇri v´ ypoˇctu bez zaokrouhlov´an´ı z´ısk´ame po koneˇcn´em poˇctu krok˚ u pˇresn´e koeficienty). Zaokrouhlovac´ı chyby vˇsak m˚ uˇzou vypoˇc´ıtan´e koeficienty hodnˇe odd´alit od jejich pˇresn´ ych hodnot. Proto se tyto metody moc nepouˇz´ıvaj´ı.

3.1

Krylovova metoda

Charakteristickou rovnici m˚ uˇzeme zapsat ve tvaru n

p(λ) = λ +

n−1 X

bi λi = 0.

i=0

Z Cayleyovy − Hamiltonovy vˇety plyne n

A +

n−1 X

bi Ai = 0.

i=0

Tedy pro kaˇzd´ y vektor y plat´ı An y +

n−1 X i=0

bi Ai y = O.

(3.1)

Rovnice (3.1) je soustava n line´arn´ıch rovnic pro n nezn´am´ ych b0 , . . . , bn−1 . Pozn´ amka 3.1.1. K v´ ypoˇctu vektoru Ai y podle rovnice Ai y = A(Ai−1 y) je tˇreba n2 n´asoben´ı, takˇze k sestaven´ı soustavy (3.1) je tˇreba ˇr´adovˇe n3 operac´ı.

10

3.2

Faddˇ ejevova-Leverrierova metoda

Metoda se op´ır´a o fakt, ˇze souˇcet vlastn´ıch ˇc´ısel libovoln´e matice je roven jej´ı stopˇe. Algoritmus Faddˇejˇevovy-Leverrierovy metody poˇc´ıt´a jednoduch´ ym zp˚ usobem koˇreny charakteristick´e rovnice. Algoritmus 1. Je d´ana matice A ˇr´adu n. Krok 1: Poloˇzme B1 = A pak p1 = tr(B1 ) Krok 2: B2 = A(B1 − p1 I) a p2 = 12 tr(B2 ) .. . Krok n: Bn = A(Bn−1 − pn−1 I) a pn = n1 tr(Bn ) Krok n+1: Charakteristick´y polynom je ve tvaru p(λ) = λn − p1 λn−1 − . . . − pn−1 λ − pn . Pozn´ amka 3.2.1. Pro inverzn´ı matici A−1 plat´ı A−1 =

1 (Bn−1 − pn−1 I). pn

Pozn´ amka 3.2.2. D˚ ukazy konvergence popsan´ ych metod v t´eto kapitole a anal´ yzu chyb m˚ uˇzeme naj´ıt v literatuˇre, viz.[1],[10]. Pˇ r´ıklad 3.2.1. Najdˇete koeficienty charakteristick´eho polynomu uˇzit´ım F.-L. metody pro matici   8 −1 3 −1 −1 6 2 0  . A= 3 2 9 1 −1 0 1 7 B1 = A ⇒ tr(B1 ) = 30 ⇒ p1 = 30,   −165 22 −42 18  22 −139 −33 3   B2 = A(B1 − 30I) =   −42 −33 −175 −17  , 18 3 −17 −159 1 1 ⇒ p2 = tr(B2 ) = (−638) = −319, 2 2  1066 −106 146 −106 992 132 B3 = A(B2 + 319I) =   146 132 1087 70 −34 67

 −70 −34  , −67  1085

1 1 ⇒ p3 = tr(B3 ) = 4230 = 1470, 3 3   −2138 0 0 0  0 −2138 0 0   B4 = A(B3 − 1410I) =   0 0 −2138 0  0 0 0 −2138 11

1 1 ⇒ p4 = tr(B4 ) = (−8552) = −2138, 2 4 ⇒ p(λ) = λ4 − 30λ3 + 319λ2 − 1410λ + 2138. Pozn´ amka 3.2.3. F.-L. metoda je i pˇres jednoduch´ y algoritmus m´enˇe v´ yhodn´a neˇz Krylovova metoda, protoˇze vyˇzaduje skuteˇcnˇe poˇc´ıtat matice Ak pro k = 1, . . . , n.

12

Kapitola 4 Poloha a odhad vlastn´ıch ˇ c´ısel 4.1

Gerˇ sgorinovy vˇ ety

Pˇresn´a znalost vlastn´ıch ˇc´ısel dan´e matice n´as v nˇekter´ ych praktick´ ych aplikac´ıch nemus´ı zaj´ımat a staˇc´ı zn´at polohu vlastn´ıch ˇc´ısel v urˇcit´ ych oblastech komplexn´ı roviny. Tyto informace m˚ uˇzeme z´ıskat i bez pˇr´ım´ ych v´ ypoˇct˚ u vlastn´ıch ˇc´ısel dan´e matice. K nalezen´ı polohy vlastn´ıch ˇc´ısel lze pouˇz´ıt n´asleduj´ıc´ı vˇetu. Vˇ eta 4.1.1. Gerˇ sgorinova vˇ eta ’ Necht A = {aij } je ˇctvercov´a matice ˇr´adu n. Definujme ri :=

n X

j=1,j6=i

|aij |,

i = 1, . . . , n.

(4.1)

Potom kaˇzd´e vlastn´ı ˇc´ıslo λ matice A splˇ nuje aspoˇ n jednu z n´ asleduj´ıc´ıch nerovnost´ı |λ − aii | ≤ ri ,

i = 1, . . . , n.

(4.2)

Jin´ymi slovy, vˇsechna vlastn´ı ˇc´ısla matice A leˇz´ı v oblasti K=

n [

Ri ,

(4.3)

i=1

kde Ri jsou kruhy o polomˇeru ri a stˇredu aii . D˚ ukaz. Necht’ λ je vlastn´ı ˇc´ıslo matice A a x je vlastn´ı vektor odpov´ıdaj´ıci vlastn´ımu ˇc´ıslu λ. Potom ze vztahu Ax = λx nebo ze vztahu (A − λI) = 0 dostaneme (λ − aii )xi =

n X

aij xj ,

i = 1, . . . , n

j=1,j6=i

kde xi je i-t´ y prvek vektoru x. ’ Necht xk je nejvˇetˇs´ı prvek vektoru x (v absolutn´ı hodnotˇe). Protoˇze |xj |/|xk | ≤ 1 pro j 6= k, je |λ − akk | ≤

n X j=1

|akj |(|xj |/|xk |) ≤

Tedy λ leˇz´ı v kruhu {λ : |λ − akk | ≤ rk }. 13

n X

j=1,j6=k

|akj |.

(4.4)

Definice 4.1.1. Kruhy Ri := {z : |z − aii | ≤ ri }, i = 1, . . . , n, se naz´ yvaj´ı Gerˇsgorinovy kruhy v komplexn´ı rovinˇe. Pozn´ amka 4.1.1. Vˇeta n´am nezaruˇcuje, ˇze v kaˇzd´em kruhu bude nˇejak´e vlastn´ı ˇc´ıslo, pouze n´am ˇr´ık´a, ˇze vlastn´ı ˇc´ısla matice A leˇz´ı ve sjednocen´ı Gerˇsgorinov´ ych kruh˚ u. N´asleduj´ıc´ı vˇeta polohu vlastn´ıch ˇc´ısel upˇresˇ nuje. Vˇ eta 4.1.2. Gerˇ sgorinova zobecnˇ en´ a vˇ eta Necht’ r Gerˇsgorinov´ych kruh˚ u je disjunktn´ıch. Pak pr´avˇe r vlastn´ıch ˇc´ısel matice A leˇz´ı ve sjednocen´ı tˇechto kruhu. D˚ ukaz. V d˚ ukazu t´eto vˇety se pouˇz´ıva vlastnost´ı z komplexn´ı anal´ yzy, viz [2]. Pozn´ amka 4.1.2. Urˇcen´ı polohy vlastn´ıho ˇc´ısla dan´e matice pomoc´ı Gerˇsgorinov´ ych vˇet je pomˇernˇe jednoduch´e. Pro zaj´ımavost uvedeme jeˇstˇe jednu vˇetu, kter´a sice tak´e urˇcuje polohu vlastn´ıch ˇc´ısel, ale jej´ı pouˇzit´ı je uˇz sloˇzitˇejs´ı a v urˇcit´ ych pˇr´ıkladech nepraktick´e. Vˇ eta 4.1.3. Necht’ A je ˇctvercov´a (obecnˇe komplexn´ı) matice n-t´eho ˇr´adu, necht’ α je (komplexn´ı) ˇc´ıslo, pro kter´e stopa matice tr((αI − A)−1 ) 6= 0. Pak v kaˇzd´em uzavˇren´em kruhu obsahuj´ıcim ˇc´ıslo α a α ˜ , kde n , α ˜ =α− tr((αI − A)−1 ) leˇz´ı alespoˇ n jedno vlastn´ı ˇc´ıslo matice A. n (α − α ˜) Definujme r = , pak v kruhu o stˇredu a polomˇeru r leˇz´ı alespoˇ n −1 2tr((αI − A) ) 2 jedno vlastn´ı ˇc´ıslo matice A. Pozn´ amka 4.1.3. Tato vˇeta nen´ı obecnˇe zn´ama a vypl´ yva z vˇet o koˇrenech polynomi´aln´ı rovnice.D˚ ukaz viz.[9] Pˇ r´ıklad 4.1.1. Uˇzit´ım Gerˇsgorinov´ych vˇet urˇcete pˇribliˇznou polohu vlastn´ıch ˇc´ısel komlexn´ı matice   1 −1/2 1/4 −1/4  1/4 1 + 2i 0 1/4    −1/2 1/4 −1 1/2  1/4 −1/2 1/2 −2 − 2i ˇ sen´ı 1. Reˇ P r1 = Pni=1,i6=1 |a1i | = 1/2 + 1/4 + 1/2 = 1 r2 = Pni=1,i6=2 |a2i | = 1/4 + 0 + 1/4 = 1/2 r3 = Pni=1,i6=3 |a3i | = 1/2 + 1/4 + 1/2 = 5/4 r4 = ni=1,i6=4 |a4i | = 1/4 + 1/2 + 1/2 = 5/4 R1 R2 R3 R4

= {z = {z = {z = {z

: |z − 2| ≤ 1} : |z − 1 − 2i| ≤ 1/2} : |z + 1| ≤ 5/4} : |z + 2 + 2i| ≤ 5/4}

14

R2 2

R3

R1

1

-2

-1

1

0

2

-1

-2

R4 Obr´azek 4.1: Gerˇsgorinovy kruhy

Podle Gerˇsgorinov´ych vˇet tedy S leˇz´ı jedno vlastn´ı ˇc´ıslo v kruhu R1 , jedno v kruhu R2 a zbyl´a dvˇe ve sjednocen´ı kruh˚ u R3 R4 . viz obr(4.1). Uved’me pˇresnou hodnotu vlastn´ıch ˇc´ısel: λ1 = 1.9285 − i0.0446 λ2 = 1.0063 + i2.0678 λ3 = −0.9079 − i0.0855 λ4 = −2.0269 − i1.9377 coˇz pˇresnˇe odpov´ıd´a poloze urˇcen´e pomoc´ı Gerˇsgorinov´ych kruh˚ u. Pozn´amky ke Gerˇsgorinovˇe vˇetˇe 1. Ze vztahu (4.4) pro maxim´aln´ı souˇradnici |xi | m˚ uˇzeme z´ıskat odhad X X |λi | ≥ |aii | − |aij | ≥ min(|akk | − |akj |) k

j6=i

a

min |λi | ≥ (|aii | − i

j6=i

X j6=i

|aij |).

Pro matici s pˇrevl´adaj´ıc´ı diagon´alou plat´ı X X 0 < min(|aii | − |aij |) ≤ |λi | ≤ max |aij | = ||A||∞ i

i

j6=i

15

j6=i

2. K matici A m˚ uˇzeme pomoc´ı jednoduch´e podobnostn´ı transformace D−1 AD = B (D je diagon´aln´ı) z´ıskat podobnou matici B, kter´a m´a jin´e Gerˇsgorinovy kruhy. Potom vˇsechna vlastn´ı ˇc´ısla leˇz´ı v oblasti KA ∩ KB . C´ılem tˇechto transformac´ı je rozklad oblasti K na souvisl´e komponenty, pˇr´ıpadn´a izolace jednoho kruhu, ve kter´em pak m˚ uˇzeme zaruˇcit existenci pr´avˇe jednoho vlastn´ıho ˇc´ısla. 3. Pokud det(λI − A) = det(λI − AT ), m˚ uˇzeme vytvoˇrit Gerˇsgorinovy kruhy i pro matici T A a z´ıskat oblast KA ∩ KAT .,ve kter´e vlastn´ı ˇc´ısla leˇz´ı.     3 2 3 −2 Pˇ r´ıklad 4.1.2. Matice A1 = resp. A2 = maj´ı stejn´e oblasti KA1 = 1 1 −1 1 KA2 := KA . Na obr.2 vid´ıme, ˇze v pˇr´ıpadˇe matice A2 , ˇza´dn´y z mal´ych kruh˚ u neobsahuje vlastn´ı ˇc´ıslo.

Obr´azek 4.2:

Pouˇzit´e znaˇcen´ı: Hranice oblasti KA je znaˇcena pˇreruˇsovanˇe Hranice oblasti KAT je znaˇcena plnou ˇcarou ˇ Sedou barvou je znaˇcena hranice √ oblasti KA ∩ KAT ⋆ vlastn´ı ˇc´ısla A1 λ1,2 = 2 ± 3 • vlatn´ı ˇc´ısla A2 λ1,2 = 2 ± i

16

Kapitola 5 Metody v´ ypoˇ ctu dominantn´ıho vlastn´ıho ˇ c´ısla ´ Umluva: Oˇc´ıslujeme-li vlastn´ı ˇc´ısla dan´e matice A tak, aby platilo |λ1 | ≥ |λ2 | ≥ . . . ≥ |λn | (kaˇzd´e ˇc´ıslo p´ıˇseme tolikr´at, kolik ˇcin´ı jeho n´asobnost), pak budeme vlastn´ı ˇc´ıslo λ1 naz´ yvat dominantn´ı vlastn´ı ˇc´ıslo.

5.1

Mocninn´ a metoda

Mocninn´a metoda je nejˇcastˇeji pouˇz´ıvanou metodou pro nalezen´ı dominantn´ıho vlastn´ıho ˇc´ısla a pˇr´ısluˇsn´eho vlastn´ıho vektoru dan´e matice. Metoda je obzvlaˇstˇe vhodn´a pro ˇr´ıdk´e matice, protoˇze spoˇc´ıv´a pouze v n´asoben´ı sloupcov´ ych vektor˚ u dan´e matice. Z´akladn´ı pˇredpoklad k uˇzit´ı t´eto metody je, ˇze dan´a matice m´a dominantn´ı vlastn´ı ˇc´ıslo λ1 a ˇze nem´a neline´arn´ı element´arn´ı dˇelitele, tj. ˇze existuje n line´arnˇe nez´avisl´ ych vlastn´ıch vektor˚ u t´eto matice, kde n je ˇr´ad matice. Konstrukce: Necht’ x je libovoln´ y vektor, x ∈ Rn , za pˇredpokladu, ˇze {v1 , . . . , vn } je mnoˇzina line´arnˇe nez´avisl´ ych vlastn´ıch vektor˚ u, m˚ uˇzeme vektor x vyj´adˇrit jako line´arn´ı kombinaci vektor˚ u vi , i = 1, . . . , n x=

n X

αi vi .

(5.1)

i=1

N´asoben´ım obou stran rovnice (5.1) maticemi A, A2 , . . . , Ak dostaneme syst´em rovnic Ax =

n X

αi Avi =

A2 x =

αi A2 vi =

A x=

n X

n X

αi λ2i vi ,

i=1

i=1

k

αi λi vi ,

i=1

i=1

n X

n X

k

αi A vi =

n X i=1

i=1

17

.. . αi λki vi .

(5.2)

Pro λk1 , kter´e jsme vypoˇc´ıtali ze syst´emu (5.2), dost´av´ame k

A x=

λk1

n X

αi (

i=1

λi k ) vi . λ1

Z pˇredpokladu, ˇze λ1 je dominantn´ı vlastn´ı ˇc´ıslo a tedy |λ1 | > |λj | j = 2, . . . , n, plyne, ˇze lim (

k→∞

λj k ) =0 λ1

a tedy lim Ak x = lim λk1 α1 v1 .

k→∞

k→∞

(5.3)

Tento postup bude konvergovat k nule, jestliˇze |λ1 | < 1 a divergovat, jestliˇze |λ1 | ≥ 1, ovˇsem za pˇredpokladu, ˇze α1 6= 0.

Pozn´ amka 5.1.1. Popsan´a konstrukce je i d˚ ukazem n´asleduj´ıc´ı vˇety.

Vˇ eta 5.1.1. Von Mises Jestliˇze matice A m´a n line´arnˇe nez´ avisl´ych vektor˚ u a je-li vlastn´ı ˇc´ıslo λ1 dominantn´ı a pro vektor x0 ∈ Rn plat´ı, ˇze hx0 , v1 i = 6 0 .Pak lim (

k→∞

Ak x0 ) = α1 v1 . λk1

(5.4)

D˚ usledek 5.1.1. Je-li y libovoln´y vektor, kter´y nen´ı ortogon´aln´ı k vlastn´ımu vektoru v1 , plyne z vˇety 5.1.1,ˇze λ1 = lim ( k→∞

yT xk+1 ), yT xk

kde xk+1 = Axk = Ak x0 . ˇ ısla yT xk+1 = yT Axk se naz´ Definice 5.1.1. C´ yvaj´ı Schwarzov´ymi konstantami. Algoritmus 2. Je zad´ana matice A Krok 1: Zvol´ıme x0 Krok 2: Pouˇzijeme iteraˇcn´ı formuli xd k+1 = Axk

Krok 3:

xk+1 =

xd k+1 (j) cj |} ⇒ λ1 = maxj=1,...,n {|x n d (j) max{|xk+1 |}

.. . d (j) (j) Krok n: Zastaven´ı v´ypoˇctu po n kroc´ıch ⇒ λ1 = maxj=1,...,n {|xn |} (k+1) (k) nebo zastaven´ı v´ypoˇctu pro |λ1 − λ1 | < δ. 18

Pozn´ amka 5.1.2. Nejˇcastˇejˇs´ı volbou poˇc´ateˇcn´ıho vektoru x0 je vektor x0 = (1, . . . , 1)T . Pˇ r´ıklad 5.1.1. Najdˇete dominantn´ı vlastn´ı ˇc´ıslo  4 2 3 3 0 4  A= 1 2 5 2 6 0 3 6 5 ˇ sen´ı 2. Zvol´ıme x0 = (1, 1, 1, 1, 1)T Reˇ   13 11    xb1 = Ax0 =  12 , 11 17 

 9.7647  8.7647    , 9.5882 xb2 = Ax1 =     7.7059  12.3529

(1)

λ1

matice  2 2 1 3  0 4 . 2 1 1 2

  0.7647 0.6471   , 0.7059 = 17 x1 =    0.6471 1

(2)

λ1 = 12.3529

•



 0.7905 0.7095   , 0.7762 x2 =    0.6238 1

• •

x10

xc 11 = Ax10 Vlastn´ı ˇc´ısla matice A jsou

  0.7731 0.6957    = 0.7735 , 0.6125 1

  10.0285  9.0247     =  10.0307 ,  7.9454  12.9722

(11)

λ1

= 12.9722.

λ1 = 12.9722, λ2 = 3.8755, λ3 = −3.0794, λ4 = −0.0297 − i0.0164. Takˇze je vidˇet,ˇze po jeden´acti kroc´ıch jsme dostali pˇresn´e ˇreˇsen´ı zadan´eho pˇr´ıkladu. 19

Pˇ r´ıklad 5.1.2. Pro matici

  1.5 −2 0.4 A =  3 0.86 −0.5 2 1.5 1.5

vˇsak metoda nebude konvergovat, protoˇze ˇc´ıseln´e hodnoty budou oscilovat. λ1 = 2.13746 λ2,3 = 0.86127 ± i2.25118

⇒ |λ2,3 | = 2.66

Absolutn´ı hodnoty vlastn´ıch ˇc´ısel jsou si rovny a tedy mocninn´ a metoda nedok´ aˇze urˇcit dominantn´ı vlastn´ı ˇc´ıslo. Pozn´ amka 5.1.3. Nev´ yhody mocninn´e metody: • odhad chyby • konvergence (obvykle v praxi nev´ıme, zda jsou splnˇeny pˇredpoklady mocninn´e metody) • volba x0 (bude-li vektor x0 takovou line´arn´ı kombinac´ı vlastn´ıch vektor˚ u, ˇze koeficient u vlastn´ıho vektoru odpov´ıdaj´ıc´ıho dominantn´ımu vlastn´ımu ˇc´ıslu bude roven 0, potom mocninn´a metoda nevypoˇcte dominantn´ı vlastn´ı ˇc´ıslo). Pozn´ amka 5.1.4. Rychlost konvergence mocninn´e metody z´avis´ı hlavnˇe na volbˇe vektoru |λ2 | x0 a na velikosti pod´ılu . |λ1 |

5.2

Metoda Rayleighova pod´ılu

Metoda Rayleighova pod´ılu je modifikovanou mocninnou metodou a zamˇeˇruje se na v´ ypoˇcet dominantn´ıho vlastn´ıho ˇc´ısla symetrick´e matice. Pro tuto ˇc´ast tedy budeme vˇzdy pˇredpokl´ad´at, ˇze matice A je symetrick´a. Potom vlastn´ı vektory mus´ı b´ yt ortonorn´aln´ı (tj. vTi vj = 0 pro i 6= j, vTi vi = 1). Odvozen´ı: 1. Zvol´ıme x0 jako line´arn´ı kombinaci vlastn´ıch vektor˚ u x0 =

n X

αi vi .

i=1

2. Sestroj´ıme posloupnost xk = Axk−1 , xk = Ak x0 ,

xk = α1 Ak v1 + . . . + αn Ak vn .

20

3. Plat´ı Avi = λi vi , potom xk = α1 λk1 v1 + α2 λk2 v2 + . . . + αn λkn vn , kde λ1 je dominantn´ı vlastn´ı ˇc´ıslo. 4. Dostaneme xk =

λk1 [α1 v1

+

n X

αi (

i=2

Sumu

Pn

i=2

αi (

λi )vi ]. λ1

λi )vi ] definujme jako wk , wk → o. λ1

5. Analogicky xk+1 6. Vyj´adˇr´ıme souˇcin xTk xk , xTk xk = λk1 [α1 v1 +

n X i=2

n

αi (

n

X X λi T k λi λi 2 )vi ]λ1 [α1 v1 + [α + αi ( )vi ] = λ2k αi2 ( )2k ] = 1 1 λ1 λ1 λ1 i=2 i=2 2 T λ2k 1 [α1 + wk wk ]

a souˇcin xTk xk+1 xTk xk+1

=

λk1 [α1 vT1

+

n X i=2

[α12 λ2k+1 1

+

n X i=2

n

X λi λi αi ( )k+1 vi ] = [α v + αi ( )k vTi ]λk+1 1 1 1 λ1 λ1 i=2

αi2 (

λi 2k+1 2 T ) ] = λ2k 1 [α1 + wk wk+1 ]. λ1

Dost´av´ame →0

xTk Axk xTk Axk+1 = lim k→∞ xT xk k→∞ xTk xk k lim

z }| { T 2 w λ2k+1 (α + k wk+1 ) = 1 2k 2 1 = λ1 . T λ1 (α1 + wk wk+1 ) | {z } →0

Pozn´ amka 5.2.1. Souˇcin wTk wk konverguje k nule pro k → ∞ dvakr´at rychleji neˇz wk k nulov´emu vektoru, z toho vypl´ yv´a, ˇze metoda Raleighova pod´ılu bude rychlejˇs´ı neˇz mocninn´a metoda. Pˇ r´ıklad 5.2.1. Metodou Rayleighova pod´ılu urˇcete dominantn´ı vlastn´ı ˇc´ıslo matice   1 1 0 A = 1 1 1 . 0 1 1 21

ˇ sen´ı 3. x0 = (1 1 Reˇ

1)T   2 x1 = Ax0 = 3 , 2

λ1 =

  5  x2 = Ax1 = 7 , 5

λ1 =

  12 x3 = Ax2 = 17 , 12

λ1 =

(1)

(2)

xT0 x1 = 2.3333, xT0 x0

xT1 x2 = 2.4118, xT1 x1

(3)

xT2 x3 = 2.4142. xT2 x2

Vlastn´ı ˇc´ısla matice A jsou

λ1 = 2.4142, λ2 = 1, λ3 = −0.4142. Tedy uˇz po tˇrech kroc´ıch jsme dostali pˇresn´e ˇreˇsen´ı.

5.3

V´ ypoˇ cet dalˇ s´ıch vlastn´ıch ˇ c´ısel mocninnou metodou

Pokud jiˇz zn´ame vlastn´ı ˇc´ıslo λ1 matice A a k nˇemu pˇr´ısluˇsn´ y vlastn´ı vektor v1 , m˚ uˇzeme vypoˇc´ıtat n´asleduj´ıc´ı vlastn´ı ˇc´ıslo λ2 a vlastn´ı vektor v2 opˇet mocninnou metodou, kterou pouˇzijeme na redukovanou matici. Vˇ eta 5.3.1. O redukci Necht’ λ1 6= 0 je vlastn´ı ˇc´ıslo matice A s vlastn´ım vektorem v1 a vektor x je libovoln´y vektor s vlastnost´ı xT v1 = 1. Potom vlastn´ı ˇc´ısla matice B = A − λ1 v1 xT jsou 0, λ2 , . . . , λn (kde λ1 , λ2 , . . . , λn jsou vlastn´ı ˇc´ısla matice A). D˚ ukaz. Necht’

  λ1 δ1 0 · · · 0  0 λ2 δ2 · · ·  .  . . .   .. .. .. J = V−1 AV =  .. , .  . . . .. .. . . δn−1   .. 0 0 0 . . . λn

je Jordan˚ uv tvar matice, kde δi ∈ {0, 1}, i = 1, . . . , n − 1. Jsou-li v1 , . . . , vn sloupce matice V, potom matice C = V−1 BV m´a tvar C = J − λ1 V−1 v1 xT V = J − λ1 e1 (xT v1 , . . . , xT vn ) = 22

= J − λ1



xT v2 . . . xT vn 0n−1,n−1

1 01,n−1



=

  0 δ1 − λ1 xT v2 −λ1 xT v3 · · · −λ1 xT vn 0  λ2 δ2 ··· 0 .  . ... ...   .. 0 =  ..  .  . . . .. .. ..  .. δn−1  0 0 0 ... λn

coˇz vˇetu dokazuje (vlastn´ı ˇc´ısla jsou na diagon´ale). V´ ybˇ er vektoru x: Vˇeta o redukci zaruˇcuje ˇsirok´ y v´ ybˇer vektoru x. Napˇr.

1. Wielandtova redukce V´ yhoda t´eto metody je v tom, ˇze v kaˇzd´e dalˇs´ı f´azi pracujeme s menˇs´ı matic´ı a prov´ad´ıme m´enˇe v´ ypoˇct˚ u. Poloˇz´ıme x=

1 j T v r λ1 1 j

kde rj je j-t´ y ˇr´adek matice A a v1j 6= 0. Index j vybereme tak, aby odpov´ıdal nejvˇetˇs´ı sloˇzce vektoru x. 2. Hotellingova redukce Zde poloˇz´ıme x = y1 , kde y1 je lev´ y vlastn´ı vektor k λ1 a je normalizov´an, tak, ˇze T plat´ı y1 x = 1. Protoˇze y1 obvykle nezn´ame, pouˇz´ıv´a se tato metoda nejsnadnˇeji u symetrick´ ych matic, v tomto pˇr´ıpadˇe je xi = vi .

23

Kapitola 6 Metody pro v´ ypoˇ cet vlastn´ıch ˇ c´ısel a vlastn´ıch vektor˚ u symetrick´ ych matic 6.1

Jacobiho metoda

Jacobiho metoda m˚ uˇze naj´ıt vˇsechna vlastn´ı ˇc´ısla a jim odpov´ıdaj´ıc´ı vlastn´ı vektory symetrick´e matice A. Metoda je vhodn´a hlavnˇe pro pln´e matice. Necht’ A je symetrick´a, potom existuje ortonorm´aln´ı b´aze sloˇzen´a z vlastn´ıch vektor˚ u A = MT DM λi jsou re´aln´a vlastn´ı ˇc´ısla matice A, D = diag(λ1 , . . . , λn ) a T je ortogon´aln´ı matice. Pˇri prvn´ım kroku Jacobiho metody poloˇz´ıme A = A1 a sestroj´ıme posloupnost {Sk }k≥1 element´arn´ıch ortogon´aln´ıch matic takovou, aby Ak+1 = STk Ak Sk = (S1 . . . Sk )T A(S1 . . . Sk ) k = 1, 2, . . . konverguj´ıc´ı k D. Protoˇze Ak+1 jsou podobn´e matici A, maj´ı stejn´a vlastn´ı ˇc´ısla. Necht’ S je matice tvaru   1 ··· 0 ··· 0 ··· 0 .. .. ..   .. . . . . . . .   0 · · · cos α · · · sin α · · · 0 . .. .. ..  ... S= . . .  ..  0 · · · − sin α · · · cos α · · · 0   .  . . . . . . . . . . . . . . 0 ··· 0 ··· 0 ··· 1 (tzn. matice rovinn´e rotace nebo Givensova transformace)

kde prvky cos α jsou na pozc´ıch (p,p) a (q,q),sin α na pozici (p,q) a − sin α na pozici (q,p). Pak plat´ı vˇeta Vˇ eta 6.1.1. Necht’ p,q jsou pˇrirozen´a ˇc´ısla, 1 ≤ p < q ≤ n, α je re´aln´e ˇc´ıslo, necht’ S je ortogon´aln´ı matice. 24

1. Je-li A = (aij ) symetrick´a, je B = ST AS = (bij ) symetrick´a a n X

b2ij

n X

=

a2ij

i,j=1

i,j=1

2. Je-li apq 6= 0, existuje jedin´e α ∈ h−π/4, 0) ∪ (0, π/4) tak, ˇze bpq = 0, kde α je jedin´e ˇreˇsen´ı rovnice aqq − app 2apq

cotg 2α = leˇz´ıc´ı v t´eto mnoˇzinˇe. Potom

n X i=1

b2ii

=

n X

a2ii + 2a2pq .

i=1

D˚ ukaz. 1. Protoˇze A = SBST a v´ıme, ˇze pro dvˇe matice K,L plat´ı tr(KL) = tr(LK), m´ame

n X

a2ij = tr(AT A) = tr(SBT ST SBST ) =

i,j=2

tr(SBT BST ) = tr(ST SBT B) = tr(BT B) = n X

b2ij .

i,j=2

2. Transformace na pozic´ıch (p,q);(q,q);(p,p);(q,p) m´a tvar         cos α − sin α cos α − sin α app apq bpp bpq · = · sin α cos α aqp aqq sin α cos α bqp bqq     cos α − sin α app cos α − apq sin α apq cos α − aqq sin α · = sin α cos α app sin α + apq cos α apq sin α + aqq cos α a tedy •

bpp = app cos2 α − 2apq sin α cos α + aqq sin2 α app cos2 α + aqq sin2 α − apq sin 2α 25

•

bpq = bqp = app cos α sin α + apq sin2 α + apq cos2 α − aqq sin α cos α = apq cos 2α + 1/2(apq − aqq ) sin 2α

•

bqq = app sin2 α + 2apq sin α cos α + aqq cos2 α app sin2 α + aqq cos2 α + apq sin 2α

Stejnˇe jako v ˇc´asti (1) a2pp + a2qq + 2a2pq = b2pp + b2qq + 2b2pq pro libovoln´e α. Zvol´ıme-li α tak, aby platilo cotg 2α = − je bpq = bqp = 0 a tedy

app − aqq 2apq

b2pp + b2qq = a2pp + a2qq + 2a2pq

ostatn´ı aii = bii pro i 6= p, q. Pozn´ amka 6.1.1. • Pˇri transformaci

A → B = ST · A · S

se mˇen´ı pouze p-t´e a q-t´e ˇr´adky a sloupce, pˇresnˇeji pro libovoln´e α : – bij = aij

pro i 6= p, q

a j 6= p, q

– bpi = bip = api cos α − aqi sin α pro i 6= p, q – bqi = biq = api sin α − aqi cos α pro i 6= p, q – bpp = app cos2 α + aqq sin2 α − apq sin 2α – bqq = app sin2 α + aqq cos2 α + apq sin 2α –

1 bpq = bqp = apq cos 2α + (app − aqq ) sin 2α 2

• Pouˇzijeme-li vztahy mezi goniometrick´ ymi funkcemi, lze prvky matice B vyj´adˇrit pomoc´ı prvk˚ u matice A.

26

Postup v´ ypoˇ ctu: • Nejprve poloˇz´ıme

K=

• Oznaˇc´ıme-li t = tg α je

• D´ale

aqq − app 2apq

(= cotg 2α)

( koˇren t2 + 2Kt − 1 pro K 6= 0 t= 1 pro K = 0

c= √

1 1 + t2

s= √

(= cos α)

t 1 + t2

(= sin α)

• Pro prvky matice B plat´ı vztahy: bpi = bip = c · api − s · aqi

i 6= p, q

bqi = biq = c · aqi + s · api

i 6= p, q

bpi = bip = app − t · apq bpi = bip = aqq + t · apq

Uved’me odvozen´ı napˇr. pro bqq bqq = app sin2 α + aqq (1 − sin2 α) + apq sin 2α = aqq − (aqq + app ) sin2 α + apq sin 2α = aqq + apq (sin 2α − 2 cotg 2α sin2 α).

Protoˇze

−2 cot 2α sin2 α + sin 2α =

sin2 2α − 2 cos2 2α sin2 α 2 sin α cos α

a d´ale ˇcitatel 4 sin2 α cos2 α − 2 sin2 α cos2 α + 2 sin4 α = 2 sin2 α(sin2 α + cos2 α) = 2 sin2 α je bqq = aqq +

sin α apq = aqq + t · apq . cos α

Jeden krok Jacobiho metody: (k) M´ame-li sestrojenou matici Ak = [aij ], vybereme (p,q) tak, aby a(k) p,q 6= 0. Sestroj´ıme Sk jako ve vˇetˇe 6.1.1, urˇc´ıme α ∈ (−π/4, 0) ∪ (0, π/4) tak, aby (k)

cotg 2αk =

(k)

aqq − app

poloˇz´ıme

(k)

2apq

,

(k+1)

Ak+1 = STk ASk = [aij 27

].

Strategie pro volbu (p,q): 1. Klasick´ a Jacobiho metoda: Zvol´ıme (p,q) takov´a , aby platilo (k)

|a(k) pq | = max |aij | i6=j

a (p,q) se mˇen´ı pro r˚ uzn´a k. 2. Cyklick´ a Jacobiho metoda: Nuluj´ı se vˇsechny nediagon´aln´ı prvky cyklickou smyˇckou, napˇr. (p,q) vol´ıme (1, 2) (1, 3)

...

(1, n); (2, 3)

...

(2, n);

...

; (n − 1, n).

Zˇrejmˇe, je-li nˇekter´ y prvek nulov´ y, postupujeme d´ale (tj. vol´ıme αk = 0 nebo Sk = I) 3. Prahov´ a Jacobiho metoda: Postupujeme jako u cyklick´e Jacobiho metody, ale nediagon´aln´ı prvky, kter´e jsou v absolutn´ı hodnotˇe menˇs´ı neˇz ”jist´a” mez, kter´a se zmenˇsuje s kaˇzdou smyˇckou, se neanuluje. Pozn´ amka 6.1.2. Co se t´ yˇce konvergence, uk´aˇzeme myˇslenku d˚ ukazu pro nejjednoduˇsˇs´ı pˇr´ıpad. Oznaˇc´ıme Pn mnoˇzinu vˇsech permutac´ı ˇc´ısel 1, 2, . . . , n. Vˇ eta 6.1.2. Posloupnost matic {Ak }∞ ıskan´ych klasickou Jacobiho metodou je konverk=1 z´ gentn´ı, lim Ak = diag(λs(i) ) k→∞

pro jistou permutaci s ∈ Pn . K d˚ ukazu potˇrebujeme n´asleduj´ıc´ı lemma. Lemma 6.1.1. Bud’ X koneˇcnˇedimenzion´aln´ı normovan´y vektorov´y prostor, {xk } ohraniˇcen´a posloupnost v X, kter´ a m´a pouze koneˇcn´y poˇcet hromadn´ych bod˚ u, necht’ lim ||xk+1 − xk || = 0.

k→∞

Potom je posloupnost {xk } konvergentn´ı. D˚ ukaz. vˇety 6.1.2 (k) Oznaˇcme Ak = [aij ] = Dk + Bk ,

(k)

Dk = diag(aii ).

• Nejprve dok´aˇzeme, ˇze limk→∞ Bk = 0. Oznaˇcme X (k) Ωk = |aij |2 . i6=j

Pak plat´ı 2 Ωk ≤ n(n − 1)|a(k) pq |

28

(k) nebot’ m´ame n(n-1) nediagon´aln´ıch prvk˚ u a ˇc´ıslo |apq | je maxim´aln´ı. D´ale podle vˇety 6.1.1 (k) Ωk+1 = Ωk − 2|aij |2 ,

tedy Ωk+1 ≤ (1 − tj.

2 )Ωk n(n − 1)

lim Ωk = 0.

k→∞

• Nyn´ı dok´aˇzeme, ˇze limk→∞ (Dk+1 − Dk ) = O. Pro diagon´aln´ı prvky matice Ak+1 plat´ı   i 6= p, q, 0, (k+1) (k) (k) aii − aii = −(tg αk )apq , i = p,   (k) (tg αk )apq , i = q. (k)

Protoˇze |αk | ≤ π/4 a limk→∞ apq = 0 je d˚ ukaz proveden.

• Necht’ {Dk′ } je posloupnost, kter´a konverguje k matici D, potom tak´e limk′ →∞ Ak′ = D, protoˇze Bk′ = 0. Ak′ = Dk′ + Bk′ a lim ′ k →∞

Tedy det(λI − D) = lim det(λI − Ak′ ) = det(λI − A). ′ k →∞

Matice Ak′ a A jsou podobn´e, tedy det(λI − Ak′ ) = det(λI − A) pro vˇsechna k ′ . Takˇze D a A maj´ı stejn´e charakteristick´e polynomy, tedy i stejn´a vlastn´ı ˇc´ısla. D proto mus´ı b´ yt diagon´aln´ı, D = diag(λs(i) ) • Posloupnost {Dk }, kde Dk je vektor dimenze n2 , je ohraniˇcen´a, nebot’ ||Dk ||2 = (

n X

i,j=1

(k) |dij |2 )1/2

≤(

n X

i,j=1

(k)

|aij |2 )1/2 =

||Ak ||2 = ||A||2 Jsou tedy splnˇeny pˇredpoklady lemmatu 6.1.1 a posloupnost {Ak } konverguje. Pˇ r´ıklad 6.1.1. Klasickou Jacobiho metodou urˇcete vˇsechna vlastn´ı ˇc´ısla matice   8 −1 3 −1 1 6 2 0 , A= 3 2 9 1 −1 0 1 7 29

ˇ sen´ı 4. Maxim´aln´ı nediagon´ Reˇ aln´ı prvek (v absolutn´ı hodnotˇe) je 3 na pozici (1,3) ⇒ p = 1 q=3 a33 − a11 1 K= = 6= 0 ⇒ 2a13 6 t je koˇren (s menˇs´ı absolutn´ı hodnotou) polynomu t2 +

1 − 1 = 0, 3

t = 0.84712708838304, 1 t c= √ = 0.76301998247272 s = √ = 0.64637489613020 2 1+t 1 + t2 b13 = b31 = 0 b11 = a11 − t · a13 = 5.45861873485088, b33 = a33 + t · a13 = 11.54138126514912, b12 = c · a12 − s · a32 = −2.05576977473312 = b21 , b14 = c · a14 − s · a34 = −1.40939487860292 = b41 , b32 = c · a32 + s · a12 = 0.87966506881525 = b23 , b34 = c · a34 + s · a14 = 0.11664508634253 = b43 , b22 = a22

b44 = a44

b42 = b24 = a24 .

Pak dostaneme matici   5.45861873485088 −2.05576977473312 0 −1.40939487860292  −2.05576977473312 6 0.87966506881525 0 .   0 0.87966506881525 11.54138126514912 0.11664508634253  −1.40939487860292 0 0.11664508634253 7 Nyn´ı opˇet vybereme maxim´aln´ı prvek a stejn´ym zp˚ usobem postupujeme d´al .. . Po 7 kroc´ıch se dostamene k matici   3.79407218081762 0.07086171427580 −0.00393661412823 0.00516622055919  0.07086171427580 6.40219536739289 −0.08436498867668 −0.06428537120075 . B= −0.00393661412823 −0.08436498867668 11.76776520507119  o 0.00516622055919 −0.06428537120075 0 8.03596724671830 Zde uˇz je vidˇet , ˇze nediagon´aln´ı prvky konverguj´ı k nule.Po dalˇs´ıch sedmi kroc´ıch uˇz dostaneme diagon´ aln´ı matici   3.2957 0 0 0  0 6.5923 0 0  , B=  0 0 11.7043 0  0 0 0 8.4077 kde diagon´aln´ı prvky odpov´ıdaj´ı vlastn´ım ˇc´ısl˚ um zadan´e matice A. 30

Nyn´ı se budeme zab´ yvat konvergenc´ı vlastn´ıch vektor˚ u klasick´e Jacobiho metody, kterou dok´aˇzeme pomoc´ı n´asleduj´ıc´ı vˇety. Pˇripomeˇ nme, ˇze Ak+1 = STk Ak Sk = QTk AQk kde Qk = S1 . . . Sk . Vˇ eta 6.1.3. Pˇredpokl´adejme, ˇze vˇsechna vlastn´ı ˇc´ısla matice A jsou vz´ajemnˇe r˚ uzn´a. Potom posloupnost matic Qk , k = 1, 2 . . . , konstruovan´ych klasickou Jacobiho metodou konverguje k ortogon´aln´ı matici, jej´ıˇz sloupce tvoˇr´ı ortogon´aln´ı mnoˇzinu vlastn´ıch vektor˚ u matice A. D˚ ukaz. Opˇet pouˇzijeme lemma 6.1.1, ovˇeˇr´ıme jeho pˇredpoklady. • {Qk } m´a pouze koneˇcn´ y poˇcet hromadn´ ych bod˚ u, kter´e jsou nutnˇe ve tvaru [±ps(1) ± ps(2) ± . . . ± ps(n) ],

s ∈ Pn ,

kde p1 , . . . , pn jsou sloupce ortonorm´aln´ı matice Q, pro n´ıˇz QT AQ = diag(λi ). Necht’ {Qk′ } je podposloupnost posloupnosti {Qk }, Qk′ → Qk . Podle vˇety 6.1.2 existuj´ı s ∈ Pn tak, ˇze (QTk′ Ak′ Qk′ ) = QTk′ Ak′ Qk′ diag(λs(i) ) = lim Ak′ = lim ′ ′ k →∞

k →∞

coˇz bylo dok´az´ano. Vˇsechna vlastn´ı ˇc´ısla jsou r˚ uzn´a, tedy existuje pouze koneˇcnˇe mnoho hromadn´ ych bod˚ u. • Pro u ´hly urˇcuj´ıc´ı Sk m´ame (k)

tg 2αk =

2apq

(k)

(k)

aqq − app

,

|αk | ≤ π/4.

Podle vˇety 6.1.2 odtud plyne, ˇze existuje l tak, ˇze pro k ≥ l je (k) |a(k) qq − app | ≥

1 min |λi − λj | > 0. 2 i6=j

(k)

Protoˇze se dvojice (p,q) mˇen´ı s k, nem˚ uˇzeme dok´azat, ˇze posloupnosti aqq konverguj´ı. Ale lim a(k) pq = 0,

(k)

a app

k→∞

tedy lim αk = 0 a

k→∞

lim Sk = I

k→∞

Qk+1 − Qk = Qk (Sk − I) → 0.

A koneˇcnˇe posloupnost {Qk } je ohraniˇcen´a, protoˇze ||Qk || = 1. Pozn´ amka 6.1.3. Pˇri v´ ypoˇctu m˚ uˇzeme pr˚ ubˇeˇznˇe kontrolovat v´ ysledky t´ım, ˇze po kaˇzd´em kroku zjiˇst’ujeme, zda (k) a(k+1) + a(k+1) = a(k) pp qq pp + aqq . Nebo vypoˇc´ıt´ame matici SDST , kter´a by se mˇela rovnat matici A. 31

Pozn´ amka 6.1.4. Pˇresnost Jacobiho metody z´avis´ı na tom, jak pˇresnˇe se vypoˇc´ıtaj´ı odmocniny pro urˇcen´ı sin αk a cos αk . Pozn´ amka 6.1.5. Aˇckoliv se Jacobiho metoda pouˇz´ıv´a pˇrev´aˇznˇe pro symetrick´e matice, pracuje ˇcasto dobˇre i v pˇr´ıpadˇe nesymetrick´ ych matic. V tomto pˇr´ıpadˇe ovˇsem konverguje k troj´ uheln´ıkov´e matici a m´a-li v´ ychoz´ı matice komplexn´ı vlastn´ı ˇc´ısla, je nutn´e pouˇz´ıt m´ısto matic Sk vhodn´e unit´arn´ı matice.

6.2

Householderova matice zrcadlen´ı

Definice 6.2.1. Matice tvaru H(u) : = I −

2uuT 2uuT = I − uT u kuk2

se naz´ yv´a Householderova matice (nˇekdy t´eˇz element´arn´ı zrcadlen´ı nebo Householderova transformace). Vlastnosti: • oznaˇcen´ı matice zrcadlen´ı se pouˇz´ıv´a proto, ˇze aplikujeme-li matici H(u) pro nˇejak´e u na vektor x ∈ Rn , pak je vektor H(u)x soumˇern´ y s vektorem x podle nadroviny ortogon´aln´ı k vektoru v.

Obr´azek 6.1: Householderova transformace • matice I je speci´aln´ı pˇr´ıpad Householderovy transformace. Pro u = o je H(o) = I. • kHxk2 = kxk2 pro kaˇzd´e x ∈ Rn , tj. zrcadlen´ı tedy nemˇen´ı d´elku vektoru. • Hy = y pro kaˇzd´e y ∈ P = {v ∈ Rn | vT u = 0}. 32

• H m´a jednoduchou vlastn´ı hodnotu -1 a (n − 1)-n´asobnou vlastn´ı hodnotu 1. D˚ ukaz. Protoˇze y ∈ P = {v ∈ Rn | vT u = 0} m´a n − 1 line´arnˇe nez´avisl´ ych vektor˚ u y1 , . . . , yn−1 a Hyi = yi pro i = 1, 2, . . . , n − 1, pak 1 je (n − 1)-n´asobn´a vlastn´ı hodnota a H tak´e zrcadl´ı u na -u, tj. Hu = −u. Takˇze -1 je vlastn´ı hodnota matice H, kter´a mus´ı b´ yt jednoduch´a, nebot’ H m´a pouze n vlastn´ıch hodnot. • z vˇety o spektr´aln´ım rozkladu plyne det(H) = (−1)1 · · · 1 = −1, • Matice H je ortogon´aln´ı a symetrick´a. D˚ ukaz. Symetrie plyne z  uuT T 2uuT = H(u). =I− HT (u) = IT − 2 T u u kuk D´ale plat´ı 2

H (u) =



2uuT I− T u u



2uuT I− T u u



= I2 − 4

uuT uuT uuT + 4 = I, kuk2 kuk4

a proto je matice H(u) ortogon´aln´ı. Vˇ eta 6.2.1. Pro kaˇzd´e dva vektory y, z ∈ Rn takov´e, ˇze y 6= z a kyk2 = kzk2 , plat´ı y = H(y - z)z. Jin´ymi slovy, kaˇzd´e dva r˚ uzn´e vektory o stejn´e normˇe lze pˇrev´est jeden na druh´y Householderovou transformac´ı. D˚ ukaz. Plat´ı 

2(y − z)(y − z)T  yT z − kzk22 H(y − z)z = I − z = z − 2 (y − z) = ky − zk22 ky − zk22 kyk22 + kzk22 − 2yT z ky − zk22 =z+ (y − z) = z + (y − z) = y. ky − zk22 ky − zk22 D˚ usledek 6.2.1. Jsou-li y, z dva vektory o stejn´e normˇe, potom existuje ortogon´aln´ı matice Q takov´ a, ˇze y = Qz. D˚ ukaz. Pro y 6= z staˇc´ı vz´ıt Q = H(y − z), jinak Q = I.

33

Vˇ eta 6.2.2. Pro kaˇzd´e x ∈ Rn je ( H(x + sgn(x1 )kxk2 e1 ), pro x1 = 6 kxk2 , H= I, pro x1 = kxk2 , ortogon´aln´ı matice s vlastnost´ı Hx = kxk2 e1 .

Nebo-li, aplikujeme-li vhodnou matici H na vektor x, dostaneme vektor, kter´y m´a vˇsechny sloˇzky aˇz na prvn´ı nulov´e. D˚ ukaz. Je-li x1 = kxk2 , potom z x21 = x21 + · · · + x2n plyne, ˇze x2 = · · · = xn = 0. Tedy x = x1 e1 = kxk2 e1 = Ix = Hx. Je-li x1 6= kxk2 , potom x + sgn(x1 )kxk2 e1 6= 0, takˇze vektory y = sgn(x1 )kxk2 e1 a z = x jsou r˚ uzn´e a plat´ı pro nˇe kyk2 = kxk2 = kzk2 , a odtud je y = sgn(x1 )kxk2 e1 = H(y − z)z = H(−x − sgn(x1 )kxk2 e1 )x.

Pozn´ amka 6.2.1. Pro vektor urˇcuj´ıc´ı Householderovu matici lze volit bud’ +kxk2 e1 nebo −kxk2 e1 . Z d˚ uvodu minimalizace numerick´ ych chyb vol´ıme stejn´e znam´enko jako u prvn´ı sloˇzky vektoru x. Vˇ eta 6.2.3. Pro kaˇzd´e x takov´e, ˇze kxk2 = 1, je ( H(x + sgn(x1 )w1 ), H= I,

pro x 6= e1 , pro x = e1 .

ortogon´aln´ı matice, jej´ımˇz prvn´ım sloupcem je vektor x. D˚ ukaz. Pro x = e1 je zˇrejm´ y. Necht’ tedy x 6= e1 . Protoˇze kxk2 = 1 = ke1 k2 , je podle Vˇety 6.2.2 x = H(x + sgn(x1 )e1 ) = He1 = H•1 ,

coˇz je tvrzen´ım vˇety. D´ıky tˇemto vˇet´am tedy um´ıme naj´ıt vektor u tak, ˇze dan´ y nenulov´ y vektor x se transformuje na vektor, kter´ y m´a nenulovou pouze prvn´ı sloˇzku. Pˇ r´ıklad 6.2.1. Lze

H(u)

x = (−1, −2, 7)T −−−→ (α, 0, 0)T ? √ √ √ Protoˇze kxk2 = 3 6, poloˇz´ıme u = x − kxk2 e1 = (−1 − 3 6, −2, 7)T a kuk2 = 6(18 + 6 ). D´ ale √  √ √ √    55 + 6√ 6 2 + 6 6 −7 − 21 6 −1 − 3 6 √  (−1 − 3 6, −2, 7) =  2 + 6 6 , uuT =  −2 4 −14 √ 7 −7 − 21 6 −14 49 takˇze

√ √  √ −1 − 3√6 −2 − 6√ 6 7 + 21 6 1 √  −2 − 6 √6 50 + 3 6 H(u) = 14√  . 3(18 + 6 ) 14 5+3 6 7 + 21 6

Snadno lze ovˇeˇrit, ˇze



√ H(u)x = (3 6, 0, 0)T . 34

6.3

Givensova-Householderova metoda

Jedn´a se o metodu speci´alnˇe vhodnou k hled´an´ı nˇekter´ ych vlastn´ıch ˇc´ısel symetrick´ ych matic, napˇr. vˇsech vlastn´ıch ˇc´ısel obsaˇzen´ ych v pˇredem zadan´em intervalu. Umoˇzn ˇuje poˇc´ıtat vlastn´ı ˇc´ısla s r˚ uznou pˇresnost´ı. Na druh´e stanˇe n´am neposkytuje informace o vlastn´ıch vektorech. M´a dvˇe etapy: • Householderova metoda pro redukci symetrick´e matice na tˇr´ıdiagon´aln´ı tvar. • Givensova metoda (metoda bisekce) pro v´ ypoˇcet vlastn´ıch ˇc´ısel symetrick´e tˇr´ıdiagon´aln´ı matice.

6.3.1

Householderova metoda

Necht’ A je symetrick´a matice, postupnˇe se urˇcuje n − 2 ortogon´aln´ıch matic H1 , . . . , Hn−2 , tak, aby matice Ak = HTk−1 · Ak−1 · Hk−1 = (H1 . . . Hk−1 )T · A · (H1 . . . Hk−1 ),

k = 1, . . . , n − 2

byly ve tvaru

Tud´ıˇz matice

 • • • • •   • • •   • •   • Ak =       ak →  



• • |• |• |• |• |•

• • • • • •

• • • • • •

aTk • • • • • •

• • • • • •

      •  •  •  •  • •

An−1 = (H1 . . . Hn−2 )T · A · (H1 . . . Hn−2 ) je tˇr´ıdiagon´aln´ı a tak´e podobn´a matici A. Kaˇzd´a transformace Ak → Ak+1 = HTk · Ak · Hk se prov´ad´ı pomoc´ı matice

  Ik 0 Hk = fk , 0 H

fk = H(f e k )ak byla nenulov´a. fk byl zvolen tak, aby pouze prvn´ı sloˇzka H(v kde H vk ), kde v

35

Potom zˇrejmˇe  • • • • •   • • •   • • •   • •  HTk · Ak · Hk =  |•   |•   f T  Hk ak → |•   |• |•



fk aTk H • • • • • • • • • • • • • • • • • •

• • •

• • •

      •  , •   •  •  • •

fk m´ame dalˇs´ı ˇc´ast tˇr´ıdiagon´aln´ı matice. tj.po vhodn´e volbˇe v Matici Hk m˚ uˇzeme popsat tak´e jako Householderovu matici pˇr´ısluˇsnou vektoru fk ]T . vk = [0, . . . , 0, v

M´ame dvˇe moˇzn´e volby vektoru vk : vk =

(k) [0, . . . , 0, ak+1,k

±(

n X

i=k+1

(k)

(k)

(k)

|aik |2 )1/2 , ak+2 , . . . , an,k ]T ,

(k)

(k+1)

znam´enko se vol´ı stejn´e jako je znam´enko u ak+1,k . M´ame-li urˇcen vektor vk , prvky aij k + 1 ≤ i , j ≤ n matice Ak+1 = Postupnˇe urˇc´ıme vektory

(k+1) [aij ]

urˇc´ıme n´asledovnˇe:

wk = (vTk vk )−1/2 vk , (k)

jejichˇz sloˇzky oznaˇc´ıme wi

(k)

qk = 2(I − wk wTk )Ak wk ,

, qi . Potom matice Ak+1 m´a tvar Ak+1 = Ak − wk qTk − qk wTk

tj.

(k+1)

aij

(k)

(k) (k)

(k)

(k)

= aij − wi qj − qi wj

k + 1 ≤ i, j ≤ n. Pˇ r´ıklad 6.3.1. Householderovou transfornac´ı pˇreved’te matici   4 2 2 1 2 −3 1 1  A= 2 1 3 1 1 1 1 2 na tˇr´ıdiagon´aln´ı tvar.

36

,

ˇ sen´ı 5. Reˇ


T √ v0 = 0 2 + ( 22 + 22 + 12 ) 2 1 ,


T 1 w0 = (vT0 v0 )− 2 v0 = 0 0.912871 0.365148 0.182574 ,


T q0 = 2(I − w0 wT0 )Aw0 = 5.477224 −2.7386095 5.03904626 3.61496813 ,   4 −3 0 0 −3 2 −2.6 −1.8   A1 = A − w0 qT0 − q0 wT0 =   0 −2.6 −0.68 −1.24 , 0 −1.8 −1.24 0.68 p


−2.62 + (−1.8)2 −1.8 ,
T w1 = 0 0 −0.954514 −0.298168 ,
T q1 = 0 6.0365793 −0.3770516 1.207794 ,   4 −3 0 0 −3 2 3.162278 0  . A2 = A1 − w1 qT1 − q1 wT1 =   0 3.162278 −1.4 −0.2 0 0 −0.2 1.4 v1 = 0 0 −2.6 −

6.3.2

Givensova metoda

Metoda slouˇz´ı k urˇcen´ı vlastn´ıch ˇc´ısel  b1 c1   B=  

symetrick´e tˇr´ıdiagon´aln´ı matice  c1  b2 c2   ... ... ... .  cn−2 bn−1 cn−1  cn−1 bn

Pokud je nˇekter´e z ˇc´ısel ci nula, rozpad´a se matice B na dvˇe tˇr´ıdiagon´aln´ı matice stejn´eho typu. Tedy bez u ´jmy na obecnosti m˚ uˇzeme pˇredpokl´adat, ˇze ci 6= 0 , (i = 1, . . . , n − 1). Oznaˇcme   b1 c 1  c 1 b2 c 2      ... ... ... Bi =  ,    ci−1 bi−1 ci−1  ci−2 bi i = 1, . . . , n

Vˇ eta 6.3.1. Polynomy pi (λ), λ ∈ R, definovan´e pro i = 1, . . . , n rekurentnˇe p0 (λ) = 1 p1 (λ) = b1 − λ

pi (λ) = (bi − λ)pi−1 (λ) − c2i−1 pi−2 (λ),

maj´ı n´ asleduj´ıc´ı vlastnosti:

37

2≤i≤n

1. Polynom pi je charakteristick´y polynom matice Bi (pi (λ) = det(Bi − λI)). 2. lim pi (λ) = +∞,

λ→∞

i = 1, . . . , n

3. Jestliˇze pi (λ0 ) = 0, potom pi−1 (λ0 )pi+1 (λ0 ) < 0, i = 1, . . . , n − 1 4. Polynom pi m´a vz´ajemnˇe i r˚ uzn´ych koˇren˚ u, kter´e oddˇeluj´ı i + 1 koˇren˚ u polynomu pi+1 , i = 1, . . . , n. D˚ ukaz.

1. Plyne z rozvoje det(Bi − λI)

2. pi (λ) = (−1)i λi − . . . → ∞ pro λ → ∞ 3. Necht’ pi (λ0 ) = 0 pro nˇejak´e i, i = 1, . . . , n − 1, z definice pi plyne pi+1 (λ0 ) = −c2i · pi−1 (λ0 ). Protoˇze ci 6= 0, dostaneme bud’ pi−1 (λ0 ) · pi+1 (λ0 ) < 0 nebo pi−1 (λ0 ) = pi− (λ0 ) = pi+1 (λ0 ) coˇz by indukc´ı vedlo k tomu, ˇze pi (λ0 ) = pi−1 (λ0 ) = . . . = p1 (λ0 ) = p0 (λ0 ), coˇz je spor, protoˇze p0 (λ0 ) = 1. 4. Plyne z 2 a 3. Pozn´ amka 6.3.1. Posloupnost polynom˚ u splˇ nuj´ıc´ı 2-4 se naz´ yv´a Sturmova posloupnost (pouˇz´ıv´a se pˇri v´ ypoˇctu koˇrenu polynom˚ u). Pˇ r´ıklad 6.3.2. Pomoc´ı charakteristick´eho polynomu urˇcete vlastn´ı ˇc´ısla tˇr´ıdiagon´aln´ı matice A2 z pˇr´ıkladu 6.3.1.   4 −3 0 0 −3 2 3.162278 0   A2 =   0 3.162278 −1.4 −0.2 0 0 −0.2 1.4

38

ˇ sen´ı 6. Reˇ p0 (λ) = 1 p1 (λ) = 4 − λ p2 (λ) = (−2 − λ)(4 − λ) − 9 p3 (λ) = (−1.4 − λ)[(−2 − λ)(4 − λ) − 9] − 10(4 − λ) p4 (λ) = (1.4 − λ)[(−1.4 − λ)[(−2 − λ)(4 − λ) − 9] − 10(4 − λ)] − 0.04[(−2 − λ)(4 − λ) − 9] = λ4 − 2λ3 − 29λ2 + 58λ − 22

Koˇreny polynomu p4 (λ) jsou λ1 = −5.4355 λ2 = 5.4907

λ3 = 1.4289 λ4 = 0.5159

Vˇ eta 6.3.2. Bud’ i pˇrirozen´e ˇc´ıslo, 1 ≤ i ≤ n. Pro dan´e µ ∈ R poloˇzme ( sgnpi (µ) je-li pi (µ) 6= 0, sgnpi (µ) = sgnpi−1 (µ) je-li pi (µ) = 0. Potom N (i, µ), coˇz je poˇcet znam´enkov´ych zmˇen v posloupnosti po sobˇe jdouc´ıch prvk˚ u uspoˇr´adan´e mnoˇziny N (i, µ) = {+, sgnp1 (µ), . . . , sgnpi (µ)} se rovn´a poˇctu koˇren˚ u polynomu pi , kter´e jsou menˇs´ı neˇz µ. Tato vˇeta umoˇzn ˇuje aproximaci (s libovolnou pˇresnost´ı) vlastn´ıch ˇc´ısel matice B = Bn a dokonce pˇr´ım´ y v´ ypoˇcet vlastn´ıho ˇc´ısla na dan´e pozici. (n) Pˇredpokl´adejme napˇr´ıklad, ˇze chceme aproximaci i-t´eho vlastn´ıho ˇc´ısla λi = λi matice B ( jako pˇredt´ım pˇredpokl´ad´ame, ˇze λ1 , . . . , λn jsou vz´ajemnˇe r˚ uzn´a a uspoˇr´adan´a sestupnˇe). Krok 1: Urˇc´ıme interval ha0 , b0 i, v nˇemˇz leˇz´ı ˇz´adan´e vlastn´ı ˇc´ıslo, napˇr. −a0 = b0 = ||B||∞ . Krok 2: c0 = Potom bud’

a0 + b 0 , spoˇcteme N (n, c0 ). 2 N (n, c0 ) ≥ i a λi ∈< a0 , c0 )

nebo N (n, c0 ) < i a λi ∈< c0 , b0 > t´ım z´ısk´ame interval < a1 , b1 >, v nˇemˇz leˇz´ı koˇren λi . Postupnˇe z´ısk´ame posloupnost interval˚ u < ak , bk >, k ≥ 0 takov´ ych, ˇze λi ∈< ak , bk > a bk − ak = 2−k (b0 − a0 ), k ≥ 0.

39

6.4

QR-rozklad

Definice 6.4.1. Dvojici matic Q a R nazveme QR-rozkladem matice A, pokud plat´ı, ˇze A = QR, pˇriˇcemˇz Q je ortogon´aln´ı matice a R je horn´ı troj´ uheln´ıkov´a matice. Nyn´ı uvedeme vˇety o existenci QR-rozkladu a jeho jednoznaˇcnosti. Vˇ eta 6.4.1. K libovoln´e re´aln´e matici A ∈ Rm×n , kde m ≥ n, existuje ortogon´aln´ı matice Q ∈ Rm×m a horn´ı troj´ uheln´ıkov´a matice R ∈ Rm×n tak, ˇze plat´ı A = QR. Vˇ eta 6.4.2. Jsou-li sloupce matice A ∈ Rm×n , m ≥ n, line´ arnˇe nez´ avisl´e, potom v QRrozkladu jsou matice R a prvn´ıch n sloupc˚ u matice Q urˇceny aˇz na znam´enko jednoznaˇcnˇe. D˚ ukazy obou vˇet viz [2]

6.5 6.5.1

Konstrukce QR-rozkladu QR-rozklad pomoc´ı Gram-Schmidtova algoritmu

Vˇ eta 6.5.1 (Gram-Schmidt˚ uv QR-rozklad). K libovoln´e re´aln´e matici A ∈ Rm×n , kde m ≥ n, existuje ortogon´aln´ı matice Q ∈ Rm×m a horn´ı troj´ uheln´ıkov´a matice R ∈ Rm×n s nez´ aporn´ymi prvky na diagon´ ale tak, ˇze plat´ı A = QR. V pˇr´ıpadˇe line´ arnˇe nez´avisl´ych sloupc˚ u matice A jsou prvky na diagon´ ale kladn´e. Z´akladn´ı myˇslenka d˚ ukazu: M´ame-li matici A ∈ Rm×n , pak aplikac´ı zobecnˇen´eho GramSchmidtova ortogonalizaˇcn´ıho procesu na sloupce matice A (ty mohou b´ yt line´arnˇe z´avisl´e m i nez´avisl´e) a doplnˇen´ım tˇechto vektor˚ u na b´azi v R z´ısk´ame sloupce matice Q. Uvaˇzujme matici A = (a1 | . . . |an ) sloˇzenou ze sloupcov´ ych vektor˚ u. Pak u1 = a1 ,

e1 =

u2 = a2 − pe1 a2 ,

u1 , ku1 k

e2 =

u3 = a3 − pe1 a3 − pe2 a3 ,

u2 , ku2 k

e3 =

.. . uk = ak −

k−1 X

p ej a k ,

j=1

40

ek =

u3 , ku3 k

uk , kuk k

kde pu v =

u.

Po u ´pravˇe obdrˇz´ıme vzorce pro vektory ai a1 = e1 ku1 k, a2 = pe1 a2 + e2 ku2 k, a3 = pe1 a3 + pe2 a3 + e3 ku3 k, .. . ak =

k−1 X j=1

pej ak + ek kuk k.

Oznaˇcme Q = (e1 | . . . |en ). Nyn´ı m´ame   < e1 , a1 > < e1 , a2 > < e1 , a3 > · · · < e1 , an >  0 < e2 , a2 > < e2 , a3 > · · · < e2 , an >     0 0 < e3 , a3 > · · · < e3 , an >  R = QT A =  ,   .. .. .. . . .. ..   . . . 0 0 0 . . . < en , an > nebot’ QQT = I a < ej , aj >= kuj k, < ej , ak >= 0 pro j  12  6 Pˇ r´ıklad 6.5.1. Proved’me QR-rozklad matice A = −4

> k.

 −51 4 167 −68. 24 −41

ˇ sen´ı 7. Gram-Schmidtov´ym procesem dostaneme Reˇ   12 −69 −58 6 . U = (u1 | u2 | u3 ) =  6 158 −4 30 −165

Matici Q potom z´ısk´ame jako Q=



   6/7 −69/175 −58/175 u1 u2 u3 6/175  . =  3/7 158/175 ku1 k ku2 k ku3 k −2/7 6/35 −33/35

A = QQT A = QR, takˇze

  14 21 −14 R = QT A =  0 175 −70 . 0 0 35 Algoritmus Mˇejme matici A. Poloˇzme r11 = ka1 k,

41

q1 =

a1 , r11

pro k = 2, . . . , n spoˇc´ıtejme: rjk = < qj , ak > zk = ak − 2 rkn

k−1 X

pro j = 1, . . . , k − 1,

rjk qj ,

j=1

= < zk , zn > zk qk = . rkk

Metodu lze tak´e upravit tak, ˇze zamˇen´ıme poˇrad´ı operac´ı. Tedy poloˇzme A0 ≡ A. Pak pro k = 2, . . . , n, spoˇctˇeme rkk

(k−1)

, = ak 2 (k−1)

rki = qTk ai

(k−1)

a , qk = k rkk pro i = k + 1, . . . , n,

A(k) = A(k−1) − qk rkT . Z form´aln´ıho hlediska jde o zmˇenu poˇrad´ı operac´ı, ovˇsem z numerick´eho hlediska obdrˇz´ıme kvalitativnˇe r˚ uzn´e v´ ysledky.

6.5.2

QR-rozklad pomoc´ı Householderovy matice

Vˇ eta 6.5.2 (Householder˚ uv QR-rozklad). m×n Kaˇzdou matici A ∈ R lze pomoc´ı s = min{n, m − 1} Householderov´ych matic rozloˇzit na souˇcin QR, a to tak, ˇze plat´ı 
R1  m > n,  0 T Hs · · · H2 H1 A = Q A = (R1 , 0) m < n,   R m = n. D˚ ukaz. Konstrukce QR-rozkladu Mˇejme re´alnou matici A



a11  a21  A =  ..  .

··· ··· ...

 a1n a2n   ..  . . 

am1 · · · amn

Krok 1.: Zkonstruujme Householderovu matici H1 tak, aby H1 A mˇela v prvn´ım sloupci pouze sam´e 0 s v´ yjimkou pozice (1, 1), tj. aby   ⊞ ···  0 · · ·   H1 A =  .. ..  . . .  0 ··· 42

K tomu staˇc´ı z´ıskat vektor un (dle pˇredchoz´ıho) tak, ˇze pro un uTn H1 = I − 2 T un un plat´ı



   a11 ⊞  a21   0      H1  ..  =  ..  .  .  . am1

Oznaˇcme A(1) : = H1 A. A(1) je tvaru

A(1)

0

  a11 · · ·  0 · · ·   =  .. ..  .  . .  0 ···

Krok 2.: Zkonstruujme Householderovu matici H2 tak, ˇze H2 A(1) m´a ve druh´em sloupci 0 pod pozic´ı (2, 2) pˇri zachov´an´ı poˇzadavku prvn´ıho kroku, tj.   ⊞  ···  0 ⊞ · · ·     A(2) : = H2 A(1) =  0 0 .  .. .. ..  . . .  0 0 ···

Matici H2 z´ısk´ame tak, ˇze nejdˇr´ıve zkonstruujeme Householderovu matici o rozmˇeru (m − 1) × (n − 1) T b 2 : = In−1 − 2 un−1 un−1 H uTn−1 un−1 takovou, ˇze     ⊞ a22  a32   0     b2 H  ..  =  ..  ,  .  . 0 am2 a definujme

  1 0 ··· 0 0    H2 : =  .. . b .  H2 0

T´ım z´ısk´ame matici A(2) = H2 A(1) . Analogicky pokraˇcujeme d´ale.

Pro k ≤ s. Krok k-t´y: Obecnˇe vytv´aˇr´ıme Householderovu matici T

b k : = In−k+1 − 2 un−k+1 un−k+1 H uTn−k+1 un−k+1 43

o rozmˇeru (m − k + 1) × (n − k + 1) takovou, ˇze

  ⊞ akk 0 .    bk  H  ..  =  ..  . . amk 0 

Definujeme



  Ik−1 0 Hk : = bk , 0 H

ˇcili m˚ uˇzeme spoˇc´ıtat A(k) = Hk A(k−1) .

T´ımto zp˚ usobem po s kroc´ıch obdrˇz´ıme matici A(s) , kter´a bude v horn´ım troj´ uheln´ıkov´em tvaru a bude pr´avˇe matic´ı R. Protoˇze A(k) = Hk A(k−1) k = 2, . . . , s, m´ame Poloˇzme

R = A(s) = Hs A(s−1) = Hs Hs−1 A(s−2) = · · · = Hs Hs−1 · · · H2 H1 A. QT = Hs Hs−1 · · · H2 H1 .

M´ame hledanou ortogon´aln´ı matici (nebot’ kaˇzd´a z Hi je ortogon´aln´ı). Celkem R = QT A, tj. A = QR. (Zopakujme si, ˇze Q = HT1 HT2 · · · HTs = H1 H2 · · · Hs .) Pˇ r´ıklad 6.5.2. Uvaˇzme matici

ˇ sen´ı 8. Krok 1.: Konstrukce H1 . Reˇ

  0 1 1 A = 1 2 3 . 1 1 1     ⊞ 0    H1 1 = 0  . 0 1

Potom tedy dle Pˇr´ıkladu 6.2.1 spoˇcteme   √    1 0 2 √      u3 = 1 + 2 0 = 1 , 0 1 1 takˇze

   1 1 0 0 u3 uT3  H1 = I3 − 2 T = 0 1 0 −  √12 u3 u3 √1 0 0 1 2 44

√1 2 1 2 1 2



√1 2 1  2  1 2



 0 − √12 − √12  1 − 12  = − √12 . 2 1 1 1 − √2 − 2 2

Urˇceme A(1) Krok 2.: Zkonstruujeme

 √ √ √  − 2 − 3 2√2 2 √2   = H1 A =  0 − 1−2√2 − 2−2√2  . 0 − 1+2 2 − 2+2 2

    −0, 2071 ⊞ = , −1, 2071 0       −1, 4318 1 −0, 2071 , = − 1, 2247 u2 = −1, 2071 0 −1, 2071   −0, 1691 −0, 9856 b2 = , H −0, 9856 0, 1691 b2 = H

tzn.

 1 0 0 H2 = 0 −0, 1691 −0, 9856 , 0 −0, 9856 0, 1691 

a spoˇc´ıt´ ame A(2) = H2 A(1) Pro Q nyn´ı plat´ı

Celkem tedy



 −1, 4142 −2, 1213 −2, 8284 0 1, 2247 1, 6330 = R = H2 H1 A =  0 0 −0, 5774 

 0 0, 8165 0, 5774 0, 4082 −0, 5774 . Q = H2 H1 = −0, 7071 −0, 7071 −0, 4082 0, 5774

  0 1 1 A = 1 2 3 = 1 1 1    0 0, 8165 0, 5774 −1, 4142 −2, 1213 −2, 8284 −0, 7071 0, 4082 −0, 5774  0 1, 2247 1, 6330 = QR. −0, 7071 −0, 4082 0, 5774 0 0 −0, 5774

45

6.5.3

QR-rozklad pomoc´ı Givensovy matice

Definice 6.5.1. Matice tvaru  1 ··· 0 .  .. . . . .. .  0 · · · c . .. G(i, j, c, s) : =  .  .. 0 · · · −s  . ..  .. . 2

2

0 ···

0

··· ··· ... ··· ···

 0 ··· 0 .. ..  . .  s · · · 0 .. ..  i iT j jT i iT j jT . .  = I+(c−1)(e e +e e )+s(e e −e e ), c · · · 0  .. . . ..  . . . 0 ··· 1

kde c + s = 1, se naz´ yv´a Givensova matice, kter´a n´am mezi jin´ ymi popisuje Givensovu transformaci.

Obr´azek 6.2: Geometrick´ y v´ yznam Givensovy rotace Givensovu matici znaˇc´ıme G(i, j, α). Opˇet chceme setrojit matice Q1 , Q2 , . . ., Qs tentokr´at vˇsak pomoc´ı Givensov´ ych matic tak, aby A(1) = Q1 A mˇela nuly pod prvkem (1, 1) v prvn´ım sloupci, matice A(2) = Q2 A(1) mˇela nuly pod prvkem (2, 2) ve druh´em sloupci, atd. Kaˇzdou z matic Qi lze sestrojit jako souˇcin Givensov´ ych matic – ten je moˇzn´e sestrojit takto: Q1 : = G(1, m, α)G(1, m − 1, α) · · · G(1, 3, α)G(1, 2, α) Q2 : = G(2, m, α)G(2, m − 1, α) · · · G(2, 3, α) .. . Bud’ s = min{m − 1, n}. Pak R = A(s) = Qs A(s−1) = · · · = Qs Qs−1 · · · Q2 Q1 A = QT A. Nyn´ı m´ame A = QR, kde QT = Qs · · · Q2 Q1 . To lze zformulovat do n´asleduj´ıc´ı vˇety. 46

Vˇ eta 6.5.3 (Givens˚ uv QR-rozklad). Bud’ A matice m × n a necht’ s = min{m − 1, n}. Existuje s ortogon´aln´ıch matic Q1 , . . ., Qs definovan´ych jako Qi : = G(i, m, α)G(i, m − 1, α) · · · G(i, i + 1, α). Pro Q = QT1 QT2 · · · QTs

plat´ı

A = QR, kde R je matice m × n s nulami pod hlavn´ı diagon´ alou.

Zn´azornˇeme si sch´ematicky Givensovu metodu redukce matice A ∈ R3×3 na horn´ı troj´ uheln´ıkov´ y tvar (symbol • znaˇc´ı prvky, kter´e se transformac´ı nezmˇenily, a ± znaˇc´ı prvky, kter´e se zmˇenily):     • • • ± ± ± G(1, 3, α) G(1, 2, α) A = • • • −−−−−→  0 ± ± −−−−−→ • • • • • •     ± ± ± • • • G(1, 3, α) G(2, 3, α) −−−−−→  0 • •  −−−−−→ 0 ± ± = R. 0 ± ± 0 0 ±

Pˇ r´ıklad 6.5.3. Necht’

  0 1 1 A = 1 2 3 . 1 1 1

ˇ sen´ı 9. Krok 1.: Najdˇeme c a s tak, aby Reˇ      ⊞ a11 c s . = 0 a21 −s c Nebot’ a11 = 0 a a21 = 1, mus´ı b´yt c = 0 a s = 1, tedy   0 1 0 G(1, 2, α) = −1 0 0 . 0 0 1

Pak dostaneme

    1 2 3 0 1 0 0 1 1 e = G(1, 2, α)A = −1 0 0 1 2 3 = 0 −1 −1 . A 1 1 1 0 0 1 1 1 1 

Nyn´ı najdˇeme c a s tak, aby

 Nebot’ e a11 = 1 a e a31 = 1, bude c =

c s −s c √1 2

    e a11 ⊞ . = e a31 0

as=

√1 , 2



tedy

√1 2

G(1, 3, α) =  0 − √12 47

0 1 0

 √1 2

0 .

√1 2

Celkem 

√1 2

 √ √3 2 2 3 2  −1 0  0 −1 −1 =  0 √1 1 1 1 0 − √12 2  √1 1 2

0 1 0

e = 0 A(1) = G(1, 3, α)A − √12 Krok 2.: Urˇceme c a s tak, aby 

c s −s c



(1) a22 (1) a32

!

√  2 2  −1 √ . − 2

  ⊞ = . 0

q (1) (1) Nebot’ a22 = −1 a a32 = − √12 , bude c = − 23 a s = − √13 , tedy G(1, 3, Θ)A(1)

6.5.4

 √  1 0 0 √3 q 2 2   1 2  0 − 3 − √3   0 −1 =  q  1 0 − √12 √ 0 − 3 − 23

√  √  √ 2 √32 2 2 2 2 q q   3 2  = R. −1 = 2 0    2 3 √ − 2 √1 0 0 3

Srovn´ an´ı algoritm˚ u

Pˇri v´ ypoˇctu QR-rozkladu pomoc´ı Householderovy matice je poˇcet proveden´ ych operac´ı roven ˇc´ıslu n n2 (3 − ). 3 K explicitn´ımu vyj´adˇren´ı matice Q je nav´ıc potˇreba 1 2(m2 n − mn2 + n3 ) 3 operac´ı, tedy celkem

1 2m2 n − mn2 + n3 . 3 Zat´ımco pro QR-rozklad pomoc´ı Givensovy matice je tento poˇcet dvojn´asobn´ y, tj. 2n2 (3 −

n ). 3 

c s −s c



Ovˇsem pokud v metodˇe s Givensovou matic´ı nahrad´ıme matici rotace a matice       c s 1 a a 1 odrazu maticemi a s ortogon´aln´ımi sloupci, pak se n´am podaˇr´ı s −c −a 1 1 −a sn´ıˇzit poˇcet operac´ı na u ´roveˇ n metody vyuˇz´ıvaj´ıc´ı Householderovy matice – jedn´a se o tzv. matice rychl´ e Givensovy transformace. Householderova matice m´a vˇsak tu nev´ yhodu, ˇze v matici, kterou ji n´asob´ıme, n´am zmˇen´ı vˇsechny prvky (zat´ımco Givensova matice jen i-t´ y a k-t´ y ˇr´adek), takˇze n´am napˇr´ıklad m˚ uˇze z ˇr´ıdk´e matice vytvoˇrit matici plnou. V modifikovan´e metodˇe s Gram-Schmidtov´ ym algoritmem je poˇcet operac´ı mn2 .

48

6.5.5

QR-rozklad a vlastn´ı ˇ c´ısla matice A – QR-algoritmus

Z´akladn´ı QR-algoritmus: Mˇejme matici A. Sestrojme jej´ı QR-rozklad, tj. A = A0 = Q0 R0 , urˇc´ıme A1 : = R0 Q0 . Nyn´ı sestrojme QR-rozklad matice A1 , tj. A1 = Q1 R1 , a spoˇctˇeme A2 = R1 Q1 . Takto pokraˇcujme analogicky d´ale. Jistˇe plat´ı Ak+1 = Rk Qk = QTk Ak Qk = QTk Rk−1 Qk−1 Qk = = QTk QTk−1 Ak−1 Qk−1 Qk = · · · = (Q0 Q1 · · · Qk )T A(Q0 Q1 · · · Qk ). Tedy matice A0 , A1 , . . . jsou kongruentn´ı (tj. A ≡ B ⇐⇒ A = PT BP). Nav´ıc d´ıky ortogon´alnosti matic Q0 , Q1 , . . . jsou matice A0 , A1 , . . . tak´e podobn´e (tj. A ∼ B (t´eˇz A ≈ B) ⇐⇒ A = P−1 BP). Tyto matice maj´ı d´ıky podobnosti stejn´a vlastn´ı ˇc´ısla jako matice A. Posloupnost tˇechto matic konverguje za urˇcit´ ych pˇredpoklad˚ u k horn´ı troj´ uheln´ıkov´e (resp. horn´ı blokovˇe troj´ uheln´ıkov´e) matici, kter´a m´a vlastn´ı ˇc´ısla na diagon´ale (resp. diagon´aln´ı bloky maj´ı vlastn´ı ˇc´ısla se stejnou absolutn´ı hodnotou) seˇrazena podle velikosti poˇc´ınaje nejvˇetˇs´ım vlastn´ım ˇc´ıslem. Poddiagon´aln´ı prvky (resp. poddiagon´aln´ı bloky) konverguj´ı k nule. Ovˇsem d˚ ukazy konvergence existuj´ı jen pro nˇekter´e speci´aln´ı typy matic. Napˇr´ıklad m´a-li matice A kladn´a vlastn´ı ˇc´ısla, pak Qk konverguje k jednotkov´e matici a posloupnost matic Ak k horn´ı tro´ uheln´ıkov´e matici, pˇriˇcemˇz diagon´aln´ı prvky t´eto matice jsou vlastn´ı ˇc´ısla matice A. Pˇ r´ıklad 6.5.4. Urˇceme vlastn´ı ˇc´ısla matice  2 A = 3 0

1 3 5 −3 11 9

 1 1 . 5 3

Lze snadno ovˇeˇrit, ˇze vlastn´ı ˇc´ısla matice A jsou λ1 = 1, λ2 = −2 a λ3 = 3. V´ysledky z´ıskan´e QR-algoritmem: (k)

(k)

(k)

k

a11

a22

a33

1

2,0

-1,6666667

1,6666667

5

3,1781374

-2,2260322

1,0478949

10 2,9486278

-1,9471270

0,9984996

15 3,0003596

-2,0064061

1,0000468

20 2,9991547

-1,9991527

0,9999984

25 3,0001104

-2,0001098

0,9999999

49

Ke zrychlen´ı konvergence lze vyuˇz´ıt tzv. posunut´ı a poˇc´ıtat nikoli rozklad matice Ak = Qk Rk , n´ ybrˇz matice e kR e k. Ak − σk I = Q

P˚ uvodn´ı spektrum matice A se t´ımto posune o σk (je v´ yhodn´e volit jej jako nˇejakou aproximaci vlastn´ıho ˇc´ısla; matice Ak − σk I m´a vlastn´ı ˇc´ısla λj − σk , jsou-li λj vlastn´ı ˇc´ısla matice Ak ). Zaznamen´av´ame-li velikosti posunut´ı σk , snadno ze znalosti spektra matice Ak najdeme spektrum matice A. Protoˇze jeˇstˇe (v metodˇe bez posunut´ı) Ak = (Q0 Q1 · · · Qk−1 )T A(Q0 Q1 · · · Qk−1 ), je A = (Q0 Q1 · · · Qk−1 )T Ak (Q0 Q1 · · · Qk−1 ). Vlastn´ı vektory y matice A dostaneme z vlastn´ıch vektor˚ u matice Ak podle vzorce y = Q0 Q1 · · · Qk−1 . e j ) z˚ T´ yˇz vzorec (s maticemi Q ust´av´a v platnosti i pro metodu s posunut´ımi, protoˇze pˇri posunut´ı se vlastn´ı vektory nemˇen´ı. Vol´ı-li se posunut´ı speci´alnˇe, dost´av´ame v nˇekter´ ych d˚ uleˇzit´ ych pˇr´ıpadech i kubickou konvergenci (tzn. zhruba ˇreˇceno, poˇcet platn´ ych m´ıst se v kaˇzd´em kroku pˇribliˇznˇe ztrojn´asob´ı).

Pozn´ amka 6.5.1. Nejv´ yhodnˇejˇs´ı se jev´ı upravit nejdˇr´ıve matici A do tzv. Hessenbergova tvaru (tj. aij = 0 pro j < i − 1, i, j = 1, . . . , n) pomoc´ı Gaussovy eliminace a pak na tuto upravenou matici pouˇz´ıt QR-rozklad. Konvergence je potom rychlejˇs´ı (obzvl´aˇstˇe pouˇzijeme-li metodu posunu, kde za σk vol´ıme tzv. Rayleigh˚ uv pod´ıl T

T

ek Hek /ek ek , kde H je pr´avˇe matice A v Hessenbergovˇe tvaru). Nav´ıc plat´ı, ˇze je-li matice A v Hessenbergovˇe tvaru, pak kaˇzd´a z matic Hk je tak´e v Hessenbergovˇe tvaru, a to i pˇri metodˇe posunut´ı.

50

Kapitola 7 Podm´ınˇ enost probl´ emu vlastn´ıch ˇ c´ısel D˚ uleˇzitou charakteristikou libovoln´eho probl´emu je jeho podm´ınˇenost, kter´a ud´av´a, jak v´ yznamnˇe se zmˇen´ı ˇreˇsen´ı probl´emu, pokud zmˇen´ıme vstupn´ı hodnoty. Podm´ınˇenost probl´emu vlastn´ıch ˇc´ısel m˚ uˇzeme popsat pomoc´ı tzv. glob´aln´ıho ˇc´ısla podm´ınˇenosti.

7.1

Glob´ aln´ı ˇ c´ıslo podm´ınˇ enosti

Vzhledem k zaokrouhlov´an´ı ˇreˇs´ıme ve skuteˇcnosti probl´em (A + E − λI)x = 0 nam´ısto (A − λI)x = 0.

˜ x ˜ , kter´ V d˚ usledku zaokrouhlov´an´ı pˇri v´ ypoˇctu dost´av´ame ˇreˇsen´ı probl´emu λ, y je pˇresn´ ym ˇreˇsen´ım probl´emu s poruchou ˜ x = 0, (A − EM − λI)˜ kde EM zahrnuje zaokrouhlovac´ı chyby bˇehem v´ ypoˇctu. Pozn´ amka 7.1.1. Anal´ yza stability vlastn´ıho probl´emu je velmi sloˇzit´a a d´a se uspokojivˇe prov´est jen pro jednoduch´e vlastn´ı ˇc´ıslo a nebo pro matici, kter´a je diagonalizovateln´a. Definice 7.1.1. Matice A je diagonalizovateln´a, kdyˇz existuje matice X takov´a, ˇze X−1 AX = D, kde D je diagon´aln´ı matice. Vˇ eta 7.1.1 (Bauer, Fike). Pokud je A diagonalizovateln´ a matice s vlastn´ımi ˇc´ısly λ1 , . . ., λn , potom vlastn´ı ˇc´ısla matice A + E leˇz´ı v kruhu Ki = {z; |z − λi | ≤ c(X).kEk}, kde c(X) = kXkkX−1 k je ˇc´ıslo podm´ınˇenosti matice vlastn´ıch vektor˚ u v maticov´e normˇe k.k. D˚ ukaz. Necht’ (A + E)x = λx. Potom bud’ λ = λi pro nˇejak´ y index i a potom λ ∈ Ki , a nebo λ 6= λi pro i = 1, . . . , n. Potom λI − A je regul´arn´ı matice. Matice (λI − A)−1 (λI − A − E) = I − (λI − A)−1 E 51

je singul´arn´ı. Proto podle vztahu ρ[(λI − A)−1 E] ≥ 1 plat´ı 1 ≤ k(λI − A)−1 Ek = k(λI − XDX−1 )−1 Ek = kX(λI − D)−1 X−1 Ek ≤ ≤ kXkkX−1 kkEkk(λI − D)−1 k = c(X)kEkk(λI − D)−1 k = = c(X)kEk max[ i

1 ]. |λ − λi |

Pˇritom jsme vyuˇzili skuteˇcnost, ˇze maticov´a norma diagon´aln´ı matice je dan´a jej´ım maxim´aln´ım diagon´aln´ım prvkem v absolutn´ı hodnotˇe. Odtud min |λ − λi | ≤ c(X)kEk. i

ˇ ıslo c(X) charakterizuje m´ıru odchylky poruˇsen´ Pozn´ amka 7.1.2. C´ ych vlastn´ıch ˇc´ısel v z´avislosti na velikosti poruchy kEk. D˚ usledek 7.1.1. Probl´em vlastn´ıch ˇc´ısel norm´ aln´ıch matic je dobˇre podm´ınˇen´y. D˚ ukaz. Protoˇze norm´aln´ı matice jsou unit´arnˇe podobn´e diagon´aln´ı matici a ve spetkr´aln´ı normˇe kUk2 = kU∗ k2 = kU−1 k2 = 1, potom c2 (U) = 1.

7.2

Odhad chyby vypoˇ c´ıtan´ eho vlastn´ıho ˇ c´ısla

Pˇresnost vypoˇc´ıtan´eho vlastn´ıho ˇc´ısla a vlastn´ıho vektoru ovˇeˇrujeme pomoc´ı rezidu´ı. ˜ a k Vˇ eta 7.2.1 (Odhad chyby vypoˇc´ıtan´eho vlastn´ıho ˇc´ısla). Necht’ urˇcen´e vlastn´ı ˇc´ıslo λ ˜ d´avaj´ı (pˇresn´y) rezidu´aln´ı vektor nˇemu pˇr´ısluˇsn´y vypoˇc´ıtan´y vlastn´ı vektor x ˜ x. r = A˜ x − λ˜ ˜ x ˜ jsou pˇresn´e hodnoty vlastn´ıho ˇc´ısla a vlastn´ıho vektoru matice s poruchou A + E, Potom λ, kde r˜ x E=− kxk22 a plat´ı odhad

˜ ≤ |λ − λ|

kyk2 kxk2 krk2 , |yT x|k˜ xk2

kde y je lev´y vlastn´ı vektor matice A pˇr´ısluˇsn´y vlastn´ımu ˇc´ıslu λ.

52

(7.1)

D˚ ukaz.



Kdyˇz

 r˜ x k˜ xk22 ˜x ˜ A− r = A˜ x − r = λ˜ x = A˜ x − k˜ xk22 k˜ xk22  2 kyk kxk kEk kEk 2 2 2 2 ˜ ≤ |λ − λ| . +O |yT x| kAk2 kEk2 =

dosad´ıme do nerovnosti, dost´av´ame

krk2 krk2 k˜ x∗ k2 = k˜ xk2 k˜ xk

 2  2 kEk krk kEk krk kyk kxk 2 2 2 2 2 2 ˜ ≤ = c(λ) . +O +O |λ − λ| |yT x| k˜ xk2 kAk2 k˜ xk2 kAk2 Pozn´ amka 7.2.1.

1. Pro symetrick´e matice c(λ) =

.

kyk2 kxk2 ˜ ≤ krk2 = 1 a |λ − λ| T |y x| k˜ xk2

2. V praxi pˇresn´e hodnoty x, y nezn´ame, proto se ve vztahu (7.1) nahrazuj´ı hodnotami ˜, y ˜ , tj. x k˜ yk2 krk2 yk2 k˜ xk2 krk2 ˜ / k˜ = |λ − λ| T xk2 ˜ | k˜ ˜| |˜ y x |˜ yT x

7.3

(7.2)

Relativn´ı chyba vypoˇ c´ıtan´ eho vlastn´ıho ˇ c´ısla

Pro jednoduch´e vlastn´ı ˇc´ıslo λ 6= 0 m˚ uˇzeme pomoc´ı (7.1) vyj´adˇrit relativn´ı chybu takto kBk2 |∆λ| yT Bx 1 kyk2 kxk2 kBk2 ≈ ε T = εc(λ) ≈ ≤ε T |λ| y x λ |y x| |λ| |λ| ≈ εc(λ)

kAk2 kAk2 ρ(A) = εc(λ) . |λ| ρ(A) |λ|

Vid´ıme, ˇze relativn´ı chyba vypoˇc´ıtan´eho vlastn´ıho ˇc´ısla z´avis´ı nejen na chybˇe ε a ˇc´ıslu podm´ınˇenosti c(λ), ale i na pod´ılu ρ(A)/|λ| = max1 |λi |/|λ|, (ˇc´ıslo ||A||/ρ(A) je konstanta, kter´a je stejn´a pro vˇsechny vlastn´ı ˇc´ısla matice A). Menˇs´ı vlastn´ı ˇc´ısla se poˇc´ıtaj´ı s vˇetˇs´ı relativn´ı chybou neˇz vˇetˇs´ı vlastn´ı ˇc´ısla.

53

Z´ avˇ er Numerick´e metody uveden´e v m´e pr´aci jsou u vyˇsˇs´ıch ˇr´ad˚ u matice poˇcetnˇe n´aroˇcn´e a proto se v souˇcasn´e dobˇe u ´lohy tohoto typu ˇreˇs´ı pomoc´ı softwarov´eho vybaven´ı poˇc´ıtaˇc˚ u. Proto maj´ı zvl´aˇstn´ı d˚ uleˇzitost podm´ınky konvergence jednotliv´ ych metod, kter´e mus´ı provˇeˇrit uˇzivatel pˇredem, protoˇze vˇetˇsina program˚ u tyto propl´emy netestuje. Pr´avˇe tak je d˚ uleˇzit´e prov´adˇet zpˇetnou kontrolu a ovˇeˇrovat pravdivost z´ıskan´ ych v´ ysledku.

54

Seznam literatury ´ V. N.: Computation Methods of Linear Algebra. [1] FADDEJEV, D. K., FADDEJEVOVA, Moskva : Fizmatgiz, 1963. [2] FIEDLER, M.: Speci´aln´ı matice a jejich pouˇzit´ı v numerick´e matematice. Praha : SNTL, 1981. [3] HIGHAM, NICHOLAS, J.: Accuracy and stability of numerical algoritms. Philadelphia: Society for Industrial and Applied Mathematics. [4] HORN, R. A., JOHNSON, Ch. R.: Matrix Analysis. Cambridge : Cambridge University Press, 1986. (rusk´ y pˇreklad, Moskva : Mir, 1989) ´ I.: Numerick´e metody. Brno :Masarykova univerzita, 1999. [5] HOROVA, ´ ENCYKLOPEDIE: Aplikovan´a matematika A aˇz Z. ˇ Praha :SNTL, 1978. [6] OBOROVE [7] RALSTON, A.: A First Course in Numerical Analysis. N. Y. : Mc Graw-Hill Book Company, 1965 (ˇcesk´ y pˇreklad Praha : Academia, 1973) ˇ [8] SIK, F.: Line´arn´ı algebra zamˇeˇren´a na numerickou anal´ yzu. Brno :Masarykova univerzita, 1998. ´ [9] VITASEK, E.: Numerick´e metody, Praha: SNTL, 1987 [10] WILKINSON, J. H.: The Algebraic Eigenvalue Problem. Oxford : Clarendon Press, 1965. (rusk´ y pˇreklad Moskva : Nauka, 1970) [11] Biswa Nath Datta: Numerical linear algebra and applications. Brooks and Cole Publishing Company, Pacic Grove, California, 1995.

55