Masarykova Univerzita v Brně Přírodovědecká fakulta
Přímé metody výpočtu charakteristických čísel matic Bakalářská práce
Vedoucí bakalářské práce RNDr. Ladislav Adamec, CSc.
Brno 2007
Roman Melichar
Prohlašuji, že jsem tuto bakalářskou práci vypracoval samostatně pod odborným vedením RNDr. Ladislava Adamce, CSc. a s použitím uvedené literatury. V Brně 20. 5. 2007
................................. Roman Melichar
Poděkování Na tomto místě bych rád poděkoval vedoucímu bakalářské práce RNDr. Ladislavu Adamcovi, CSc. za cenné rady, připomínky a čas, který mi věnoval. Dále bych chtěl poděkovat svým prarodičům za veškerou podporu, které se mi v průběhu studia dostalo a Janě Modlitbové za pečlivé přečtení textu a za morální oporu.
Obsah 1 Úvod
4
2 Vlastní čísla a vlastní vektory matic 2.1 Matice speciálního typu . . . . . . . . . . . . 2.2 Charakteristický polynom, vlastní čísla matice 2.3 Podobné matice . . . . . . . . . . . . . . . . . 2.4 Vlastní čísla maticového polynomu . . . . . . 2.5 Lineární operátory . . . . . . . . . . . . . . . 2.6 Vlastní vektory a vlastní čísla operátoru . . . 2.7 Konjugovaný operátor . . . . . . . . . . . . . 3 Základní přímé metody 3.1 Metoda Leverrierova 3.2 Metoda interpolační 3.3 Metoda Krylovova . 3.4 Metoda Danilevského
. . . . . . .
. . . . . . .
. . . . . . .
řešení problému vlastních . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
čísel . . . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . .
5 . 5 . 6 . 8 . 9 . 9 . 13 . 15
. . . .
18 18 20 23 30
. . . .
4 Stabilita problému vlastních čísel
36
5 Příloha - kódy algoritmů 5.1 Kód algoritmu Leverrierovy metody . . . . . . . 5.2 Kód algoritmu pro výpoče koeficientů cmi . . . . 5.3 Kód algoritmu interpolační metody . . . . . . . 5.4 Kód algoritmu Krylovovy metody . . . . . . . . 5.5 Kód algoritmu Danilevslkého metody . . . . . . 5.6 Kód algoritmu pro měření přesnosti jednotlivých
39 39 40 40 41 43 45
. . . . . . . . . . . . . . . . . . . . metod
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
Závěr
46
Tabulka použitých symbolů
48
Literatura
49
Rejstřík
50
3
Kapitola 1 Úvod Jednou z nejstarších matematicky zpracovaných úloh o vlastních hodnotách je úloha vyšetřovaná již roku 1744 Leonhardem Eulerem pojednávající o kritickém zatížení štíhlého prutu namáhaného tlakem, při němž by mohlo nastat vybočení prutu ohybem. S rozvojem klasické ”matematické fyziky” se během 19. století objevily četné úlohy o vlastních hodnotách při problémech kmitání. Minulé století pak v různých fyzikálních a technických oborech přineslo takovou hojnost úloh, které vedly právě k výpočtu vlastních hodnot, že se tato úloha stala jednou z nejdůležitějších v lineární algebře. Proto bylo věnováno poměrně velké úsilí na vypracování metod jejich výpočtu, založených na nejrůznějších myšlenkách. Prvními z takto vypracovaných metod byly metody přímé. Protože však nebyly co se týče přesnosti vhodné pro náročnější úlohy, byly později rozpracovány i nepřímé metody řešení výpočtu vlastních čísel a vektorů, jejichž použití už bylo vhodné i pro rozsáhlejší úlohy. Úkolem této práce je vyložit několik základních přímých metod výpočtu, navrhnout algoritmy implementující tyto metody a vzájemně je srovnat z hlediska vyžadovaného počtu operací a jejich přesnosti. Tyto algoritmy budou zapsány v programovacím jazyce matlab. Omezíme se zde však pouze na výpočet koeficientů charakteristického polynomu. Výpočet jeho kořenů (vlastních čísel) bude implementován jednou ze základních funkcí matlabu. Poznamenejme ještě, že všechny zde probírané metody kromě Leverrierovy (z roku 1840) vznikly ve třicátých letech tohoto století. Tato práce je členěna do čtyř kapitol. Po úvodu následuje v druhé kapitole seznámení s problematikou vlastních čísel a vektorů. Ve třetí kapitole jsou pak vyloženy čtyři základní přímé metody řešení problému vlastních čísel a ve čtvrté kapitole je krátce pojednáno o stabilitě problému vlastních čísel a vlastních vektorů. Pátou kapitolu tvoří kódy algoritmů pro jednotlivé metody. V závěru je uvedeno srovnání přesností dosahovaných jednotlivými metodami ve srovnání s nepřímou metodou implementovanou v matlabu. Všechny uváděné teoretické poznatky jsou pro jejich názornost doplněny o praktické výpočty.
4
Kapitola 2 Vlastní čísla a vlastní vektory matic V této kapitole zavedeme základní pojmy týkající se vlastních čísel a vlastních vektorů a vyšetříme jejich vlastnosti. Poté o nich pojednáme v souvislosti s lineárními operátory. Poslední část této kapitoly věnujeme některým vlastnostem vlastních čísel a vlastních vektorů konjugovaného operátoru. Cílem této kapitoly je nejen zavedení těchto pojmů, ale také vyšetření vlastností, které tvoří podstatu metod jejich určování.
2.1
Matice speciálního typu
Souhrn čísel uspořádaných do obdélníkového schématu nazýváme maticí. Zpravidla ji zapisujeme ve tvaru a11 a12 . . . a1n a21 a22 . . . a2n .. .. .. . . . am1 am2 . . . amn nebo též
A = (aij ) ,
i = 1, 2, . . . , m;
j = 1, 2, . . . , n ,
kde čísla m, n udávají počet řádků a sloupců. O této matici pak říkáme, že je typu m × n. Čísla aij (většinou reálné nebo komplexní) se nazývají prvky matice. Pokud je m = n, nazýváme pak tuto matici čtvercovou maticí řádu n a o jejich prvcích aij , pro které je i = j, mluvíme jako o diagonálních prvcích. Tyto pak tvoří tzv. hlavní diagonálu matice. Čtvercovou matici tvaru 1 0 ... 0 0 1 . . . 0 .. .. .. , . . . 0 0 ... 1 5
jejíž diagonální prvky jsou rovny jedné a jsou ostatní nulové, budeme nazývat jednotkovou maticí a budeme ji označovat E. Další velmi významnou čtvercovou maticí je matice nulová, jejíž všechny prvky jsou rovny nule. Tu budeme označovat symbolem 0. Nechť A je čtvercová matice typu n. Zaměníme-li v této matici řádky za sloupce, získáme tak matici transponovanou, kterou značíme A′ . Pokud bychom v matici A zaměnili všechny prvky za prvky k nim komplexně sdružené, dostali bychom tak matici komplexně sdruženou. Taková matice se značí symbolem A. Matici transponovanou a komplexně sdruženou k matici A nazvěme maticí konjugovanou a označme ji A∗ . Mezi čtvercovými maticemi existuje několik velmi významných skupin se specifickými vlastnostmi. Pro naše účely bude vhodné definovat pouze tři skupiny. První z nich jsou tzv. normální matice. Těmi budeme nazývat čtvercové reálné matice A, pro které platí A′ A = AA′ . Dále čtvercovou komplexní matici A nazveme hermitovskou nebo též hermitovsky symetrickou , platí-li aij = a ¯ji, tj. v maticovém zápisu A = A∗ . A konečně čtvercovou komplexní matici A nazveme unitární, je-li A∗ A = E.
2.2
Charakteristický polynom, vlastní čísla matice
Nechť A = (aij ) je čtvercová matice typu n a t je reálná proměnná. Deteminant tvaru det (A − tE), tedy a11 − t a . . . a 12 1n a21 a22 − t . . . a2n .. .. .. , . . . an1 an2 . . . ann − t
je polynomem n-tého stupně v t a nazývá se charakteristický polynom matice A. Položíme-li det (A − tE) = 0, získáme tak charakteristickou rovnici matice A. Poznamenejme, že přímý výpočet koeficientů pi charakteristického polynomu (−1)n [tn − p1 tn−1 − p2 tn−2 − · · · − pn ] je z hlediska počtu provedených operací velmi zdlouhavý, neboť jak známo platí p1 = a11 + a22 + · · · + ann , .. . pn = (−1)n−1 det A , 6
kde pk pro 2 ≤ k ≤ n − 1 je až na znaménko (−1)k−1 součtem všech minorů k-tého řádu matice A, jejichž hlavní diagonála je částí hlavní diagonály matice A. Definice 1. Nechť A je čtvercová matice řádu n. Kořeny λ1 , . . . , λn charakteristické rovnice matice A se nazývají vlastními čísly nebo též charakteristickými čísly matice A. Ze známých vztahů mezi koeficienty algebraické rovnice a jejími kořeny vyplývá, že λ1 + λ2 + · · · + λn = p1 = a11 + a22 + · · · + ann , λ1 λ2 . . . λn = (−1)n−1 pn = det A . Vyslovíme nyní zajímavou větu, která ukazuje, že k dané matici existuje nenulový polynom, který se anuluje po dosazení této matice. Věta 1 (Cayleyova-Hamiltonova). Nechť A je čtvercová matice řádu n. Je-li ϕ(t) její charakteristický polynom, potom platí ϕ(A) = 0. Důkaz. Vyjděme z označení ve výše uvedené větě a nechť dále B = adj (A − tE). Zřejmě každý algebraický doplněk matice A−tE je polynomem stupně nejvýše n−1, a proto můžeme psát B = Bn−1 + Bn−2 t + · · · + B0 tn−1 , kde Bn−1 , . . . , B0 jsou matice již neobsahující t. Vzhledem k základní vlastnosti adjungované matice, která říká, že pro každou čtvercovou matici M platí adj M.M = det M.E, můžeme psát (Bn−1 + Bn−2 t + · · · + B0 tn−1 ).(A − tE) = det (A − tE).E = = (−1)n (tn − p1 tn−1 − · · · − pn )E . což je Bn−1 A−Bn−1 t+Bn−2 At−Bn−2 t2 +· · ·+B0 Atn−1 −B0 tn = (−1)n (tn −p1 tn−1 −· · ·−pn )E. To je ekvivalentní soustavě rovností Bn−1 A = (−1)n+1 pn E , Bn−2 A − Bn−1 = (−1)n+1 pn−1 E , .. . B0 A − B1 = (−1)n+1 p1 E , −B0 = (−1)n E, Odtud násobením zprava každou z rovností maticí Ai−1 , kde i značí i-tou rovnost, a sečtením získáme 0 = (−1)n [−pn E − pn−1 A − pn−2 A2 − · · · + An ] = ϕ(A) , což jsme měli dokázat. 7
Zřejmě anuluje-li se polynom ψ(t) po dosazení nějaké čtvercové matice, pak se takto anuluje i každý jeho násobek. Nenulový polynom f (t) nejmenšího stupně s touto vlastností nazýváme minimálním polynomem matice A. Charakteristický polynom ϕ(t) libovolné čtvercové matice A je dělitelný minimálním polynomem ψ(t) této matice. Podle věty o dělení polynomů totiž existují polynomy q(t) a r(t), kde r(t) je buď nulový, nebo má menší stupeň než polynom ψ(t), tak, že ϕ(t) = ψ(t)q(t) + r(t). Po dosazení matice A do této rovnosti dostáváme r(A) = ϕ(A) − ψ(A)q(A) = 0. Protože polynom ψ(t) byl nenulovým polynomem nejmenšího stupně anulujícím se po dosazení matice A, musí být r(t) nutně polynom nulový. Charakteristický polynom ϕ(t) je tedy skutečně dělitelný minimálním polynomem ψ(t). Poznámka 2.2.1. Úplně stejným způsobem se dokáže, že libovolný polynom ω(t), pro který platí ω(A) = 0, je dělitelný minimálním polynomem matice A.
2.3
Podobné matice
Nechť A a B jsou dvě libovolné čtvercové matice stejného řádu. Existuje-li regulární matice C téhož řádu taková, že B = C −1 AC, řekneme, že matice B je podobná matici A. Zřejmě platí C −1 A1 C + C −1 A2 C + · · · + C −1 Am C = C −1 (A1 + A2 + · · · + Am )C , C −1 A1 C.C −1 A2 C . . . C −1 Am C = C −1 (A1 A2 . . . Am )C , a odtud i (C −1 AC)m = C −1 Am C . A proto pro libovolný polynom f (t) platí f (C −1 AC) = C −1 f (A)C . Odtud snadno nahlédneme, že minimální polynomy podobných matic jsou až na nenulovou multiplikativní konstantu totožné. Je-li totiž f (t) minimální polynom matice A, pak je zřejmě anulujícím polynomem matice B. Pokud by existoval minimální polynom g(t) matice C −1 AC nižšího řádu než f (t), pak z rovnosti g(C −1 AC) = 0 = C −1 g(A)C ihned plyne, že g(A) anuluje i matici A, neboť C a C −1 jsou nenulové, což vede ke sporu. Platí dokonce, že podobné matice mají stejné charakteristické polynomy. Platí totiž det (B − tE) = det (C −1 AC − tE) = det (C −1 AC − tC −1 EC) = = det (C −1 (A − tE)C) = det (C −1 ) det (A − tE) det C = = (det C)−1 det (A − tE) det C = det (A − tE) . 8
2.4
Vlastní čísla maticového polynomu
Nyní ukážeme zajímavou vlastnost maticových polynomů, na které je založena Leverrierova metoda hledání vlastních čísel. Věta 2. Nechť A je čtvercová matice řádu n a λ1 , . . . λn její vlastní čísla, z nichž některá mohou splývat. Nechť dále f (t) je libovolný polynom. Potom vlastními čísly matice f (A) jsou f (λ1 ), . . . , f (λn ). Důkaz. Pro nulový polynom f (t) věta zřejmě platí, neboť nulová matice libovolného řádu má jediné charakteristické číslo nula. Mějme tedy f (t) nenulový polynom s kořeny µ1 , . . . , µm . Potom můžeme psát f (t) = a(t − µ1 ) . . . (t − µm ) , kde a je vedoucí koeficient. Nechť ϕ(t) = (λ1 −t) . . . (λn −t) je charakteristický polynom matice A. Dosazením matice A do f (t) dostaneme f (A) = a(A − µ1E) . . . (A − µm E) a odtud už pro determinant det f (A) získáváme det f (A) =an det (A − µ1 E) . . . det (A − µm E) = an ϕ(µ1 ) . . . ϕ(µm ) = n m n m n Y Y Y Y Y n f (λi ) . (a (λi − µj )) = (λi − µj ) = =a i=1 j=1
i=1
j=1
i=1
Rovnost det f (A) = f (λ1 ) . . . (f (λn ) je identita vzhledem ke koeficientům polynomu f (t). Užijeme-li této identity pro polynom f (t) − u, dostaneme det (f (A) − uE) = (f (λ1 ) − u) . . . (f (λn ) − u). To však právě znamená, že vlastní čísla matice f (A) jsou čísla f (λ1 ), . . . , f (λn ). Poznámka 2.4.1. Speciálně pro matici Ak (k ∈ N) dostáváme vlastní čísla λk1 , . . . , λkn .
2.5
Lineární operátory
Nechť R je lineární prostor. Je-li na množině R dána funkce s oborem hodnot R, nazýváme tuto funkci operátorem nad R. Operátor, který splňuje následující podmínky: 1. A(αu) = αAu pro každý vektor u ∈ dom(A) a libovolné číslo α (reálné nebo komplexní podle druhu prostoru), 2. A(u1 + u2 ) = Au1 + Au2 pro každé dva vektory u1 , u2 ∈ dom(A), nazýváme lineárním operátorem. Definujme nyní několik základních operací s lineárními operátory nad týmž prostorem. První z nich je součin AB lineárních operátorů A a B (v tomto pořadí), pro nějž platí AB = A ◦ B. 9
Tento součin je zřejmě opět lineárním operátorem, neboť pro libovolné u1 ,u2 ∈ dom(B) platí AB(u1 + u2 ) = A(B(u1 + u2 )) = A(Bu1 + Bu2 ) = = A(Bu1 ) + A(Bu2 ) = ABu1 + ABu2 . Součinem lineárního operátoru A s číslem α nazveme operátor αA, pro který platí αAu = α(Au) pro každé u ∈ dom(A). Dále součtem lineárních operátorů A a B nazveme operátor A + B, kde (A + B)u = Au + Bu pro každé u z dom(A). Součet i číselný násobek lineárních operátorů jsou zřejmě opět lineárními operátory. Mezi lineárními operátory existují dva význačné, které ve smyslu výše uvedených operací plní funkci jednotkového a nulového prvku. Je to identický operátor E, který každému vektoru prostoru přiřazuje tentýž vektor a nulový operátor 0, který každému vektoru přiřadí nulový vektor. Pro libovolný operátor A pak zřejmě platí EA = AE = A , A+0 = 0+A= 0 . Přirozeným způsobem můžeme definovat i n-tou mocninu An lineárního operátoru A jakožto n-násobný součin stejného lineárního operátoru. Konečně definujme polynom operátoru jako operátor f (A) = a0 An + a1 An−1 + · · · + an−1 A + an E , kde f (t) je polynom f (t) = a0 tn + a1 tn−1 + · · · + an−1 t + an . V následujícím textu se budeme zabývat jen lineárními operátory a proto budeme slovo lineární vynechávat. Nyní ukážeme, že ke každému operátoru A : R → R při pevně zvolené bázi v R lze jednoznačně přiřadit matici a naopak. Nechť tedy u1 , u2 , . . . , un je libovolná báze v R a Au1 = a11 u1 + a21 u2 + · · · + an1 un , Au2 = a12 u1 + a22 u2 + · · · + an2 un , ................................. Aun = a1n u1 + a2n u2 + · · · + ann un . Pak pro libovolný vektor x = x1 u1 + x2 u2 + · · · + xn un z R platí Ax = x1 Au1 + · · · + xn Aun 10
a odtud y = Ax = = x1 (a11 u1 + · · · + an1 un ) + · · · + xn (a1n u1 + · · · + ann un ) = = (x1 a11 + · · · + xn a1n )u1 + · · · + (x1 an1 + · · · + xn ann )un . To lze maticově zapsat ve tvaru y = Ax, kde y a x jsou sloupce souřadnic vektorů y a x a A je matice složená ze sloupců souřadnic vektorů Au1 , . . . , Aun v bázi u1 , . . . , un , tj. a11 a12 . . . a1n a21 a22 . . . a2n A = .. .. .. . . . . an1 an2 . . . ann
Vzhledem k jednoznačnosti souřadnic libovolného vektoru v libovolné bázi je matice A určena jednoznačně. Naopak definujme přiřazení, které vektoru x se souřadnicemi x v nějaké bázi přiřazuje vektor y, jehož souřadnice y v té samé bázi jsou dány vztahem y = Ax. Toto přiřazení je zřejmě lineárním operátorem a je určeno jednoznačně. Toto jednoznačné přiřazení mezi maticemi a operátory je zachováno i při operacích s operátory. Matice součtu operátorů je totiž rovna součtu matic obou sčítaných operátorů a matice součinu operátorů je rovna součinu matic jednotlivých faktorů (v tom samém pořadí). Vidíme tedy, že existuje izomorfní zobrazení mezi množinou operátorů n-rozměrného prostoru a množinou čtvercových matic n-tého řádu. Toto zobrazení přiřazuje operátoru matici při pevně zvolené bázi. Je-li tedy báze prostoru pevně zvolena, můžeme ztotožnit operátor s příslušnou maticí, stejně jako lze ztotožnit vektor s příslušným sloupcem jeho souřadnic. Ve smyslu tohoto přiřazení si pak odpovídá vektor vzniklý provedením operátoru na vektoru a sloupec vzniklý vynásobením příslušné matice příslušným sloupcem. Nyní vyšetříme vztah souřadnic vektorů v různých bázích. Nechť R je lineární prostor a x ∈ R je vektor. Označme (x)α , kde α = (u1 , . . . , un ) je báze prostoru R, souřadnice vektoru x v bázi α. Nechť β = (v1 , . . . , vn ) je druhou bází tohoto prostoru. Potom matici A řádu n nazveme maticí transformace souřadnic v bázích α a β, jestliže platí (x)β = A(x)α . Podívejme se, jak vypadá tato matice transformace souřadnic. Vyjádříme-li vek-
11
tory báze α jako lineární kombinace vektorů báze β u1 = a11 v1 + · · · + an1 vn , u2 = a12 v1 + · · · + an2 vn , ............... un = a1n v1 + · · · + ann vn , snadno pak vyjádříme souřadnice vektoru x = x1 u1 + · · · + xn un v bázi β: x = x1 (a11 v1 + · · · + an1 vn )+ + x2 (a12 v1 + · · · + an2 vn )+ + ..................... + + xn (a1n v1 + · · · + ann vn ) = = (x1 a11 + x2 a12 + · · · + xn a1n )v1 + = (x1 a21 + x2 a22 + · · · + xn a2n )v2 + + .............................. + = (x1 an1 + x2 an2 + · · · + xn a1n )vn . Odtud už je zřejmé, že
Platí tedy
a11 a12 . . . a1n x1 a21 a22 . . . a2n x2 (x)β = .. .. .. . .. . . . . . an1 an2 . . . ann xn A = ( (u1 )β (u2 )β . . . (un )β ) .
Poznámka 2.5.1. Pro matici A operátoru A v bázi α zřejmě platí B = ( (Au1 )α (Au2 )α . . . (Aun )α ). Tyhle poznatky nám už dovolí lehce odvodit, jak se změní matice operátoru při změně báze prostoru. Mějme dvě matice A a B operátoru A v souřadnicích α a β, a C matici transformace souřadnic v bázích β a α. Ukážeme, že platí B = C −1 AC , tedy že matice operátoru v různých bázích jsou navzájem podobné a tato podobnost je realizována příslušnou maticí transformace. Zřejmě pro libovolný vektor x ∈ R platí (y)α = A(x)α , (x)α = C(x)β , (y)α = C(y)β . 12
Odtud vyplývá, že C(y)β = A(x)α = AC(x)β , a tedy (y)β = C −1 AC(x)β . Avšak (y)β = B(x)β a proto musí platit B = C −1 AC .
2.6
Vlastní vektory a vlastní čísla operátoru
Definice 2. Vlastním číslem (nebo také vlastní hodnotou) operátoru A nazveme každé číslo λ pro které existuje nenulový vektor x ∈ dom(A) takový, že platí Ax = λx. Každý takový nenulový vektor x nazveme vlastním vektorem operátoru A, příslušným vlastnímu číslu λ. Podívejme se, jak se najdou vlastní vektory a vlastní čísla nějakého operátoru A. K tomuto účelu nechť A je matice tohoto operátoru v nějaké bázi a nechť x1 , . . . , xn jsou souřadnice některého jeho vlastního vektoru x v této bázi. Souřadnice vektoru Ax budou zřejmě ve tvaru n X
a1k xk , . . . ,
ank xk .
k=1
k=1
Proto můžeme psát
n X
a11 x1 +a12 x2 + . . . a21 x1 +a22 x2 + . . .
+a1n xn = λx1 , +a2n xn = λx2 , .. .
an1 x1 +an2 x2 + . . . +ann xn = λxn nebo též (a11 − λ)x1 + a12 x2 +···+ a21 x1 + (a22 − λ)x2 + · · · + an1 x1 + an2 x2
a1n xn = 0 , a2n xn = 0 , .. .
+ · · · + (ann − λ)xn = 0 .
Zřejmě tato homogenní soustava lineárních rovnic o neznámých x1 , . . . , xn má nenulové řešení právě tehdy, když det (A − Eλ) = 0 a to je právě tehdy, když λ je vlastním číslem matice A. Pokud bychom zvolili matici operátoru A v jiné bázi, dostali bychom tytéž vlastní čísla, neboť, jak již víme, matice operátoru v různých bázích jsou si navzájem podobné a jejich charakteristické polynomy jsou totožné. Vidíme tedy, že vlastní čísla operátoru splývají s kořeny charakteristického polynomu matice, která tomuto operátoru v libovolné bázi odpovídá. To nás opravňuje vyslovit následující definici. 13
Definice 3. Charakteristický polynom matice operátoru v libovolné bázi se nazývá charakteristickým polynomem operátoru. Mějme A matici operátoru A a ϕ(t) jeho charakteristický polynom. Podle CayleyHamiltonovy věty platí, že ϕ(A) = 0. Polynomu ϕ(A) operátoru A tedy odpovídá nulová matice ϕ(A) a proto ϕ(A) = 0. Vidíme tedy, že můžeme Cayley-Hamiltonovu větu zobecnit i na operátory. Stejným způsobem můžeme operátoru přiřadit i jeho minimální polynom, který bude odpovídat nenulovému polynomu nejmenšího stupně, který po dosazení daného operátoru dá nulový operátor. Takový polynom je zřejmě i minimálním polynomem matice tohoto operátoru v libovolné bázi a proto je charakteristický polynom operátoru dělitelný minimálním polynomem tohoto operátoru. Každý kořen minimálního polynomu je tedy i kořenem charakteristického polynomu. Avšak platí to i obráceně, že každý kořen charakteristického polynomu je kořenem minimálního polynomu. Tedy že množiny kořenů charakteristického polynomu a minimálního polynomu jsou totožné. Abychom tohle ukázali, předpokládejme, že λ je vlastní číslo nějakého operátoru A a x je jeho vlastní vektor příslušný tomuto vlastnímu číslu. Je-li ψ(t) minimálním polynomem daného operátoru, pak podle Bézoutovy věty o zbytku platí pro nějaký polynom p(t), že ψ(t) = p(t)(t − λ) + ψ(λ). Dosazením A a násobením x pak dostáváme, že ψ(A)x = p(A)(A − λE)x + ψ(λ)x. Avšak (A − λE)x = 0 a ψ(A)x = 0x = 0. Proto je ψ(λ)x = 0. Vzhledem k tomu, že x je nenulový, musí být ψ(λ) = 0. Definice 4. Vlastním vektorem matice A nazveme každý nenulový vektor x, který vyhovuje rovnosti Ax = λx pro nějaké vlastní číslo λ matice A. Z předešlého plyne, že vlastní vektory matice A operátoru A v nějaké bázi α jsou právě sloupci souřadnic vlastních vektorů tohoto operátoru v bázi α. Platí totiž A(x)α = (Ax)α = (λx)α = λ(x)α . Na závěr této části vyšetřeme některé další užitečné vlastnosti vlastních vektorů. Předně podotkněme, že vektor, jehož souřadnice jsou komplexně sdruženy s příslušnými souřadnicemi vlastního vektoru reálné matice, je také vlastním vektorem této matice a odpovídá komplexně sdruženému vlastnímu číslu. To ukážeme snadno. ¯ Je-li totiž A reálná čtvercová matice a λ nějaké její vlastní číslo, pak je zřejmě i λ jejím vlastním číslem a platí ¯x . A¯ x = Ax = λx = λ¯ Dále ukážeme zajímavou vlastnost, že odpovídá-li několik vlastních vektorů navzájem různým vlastním číslům, jsou tyto vektory lineárně nezávislé. Abychom to dokázali, označme A libovolný operátor, λ1 , . . . , λs jeho navzájem různá vlastní čísla a x1 , . . . , xs některé jim odpovídající vlastní vektory. Předpokládejme, že tyto vektory jsou lineárně závislé. Bez újmy na obecnosti můžeme považovat vektory 14
x1 , . . . , xj , kde j < s, za lineárně nezávislé a zbylé vektory za jejich lineární kombinace. Pro xs platí j X ci xi . xs = i=1
Pak
Axs =
j X
ci Axi =
ci λi xi ,
i=1
i=1
ale zároveň také platí
j X
Axs = λs xs =
j X
λs ci xi .
i=1
Dohromady nám to tedy dává j X
ci λi xi −
j X
λs ci xi = 0
i=1
i=1
j X
(λs − λi )ci xi = 0 .
i=1
Z lineární nezávislosti vektorů x1 , . . . , xj vyplývá, že všechny koeficienty (λs − λi )ci jsou rovny nule a protože podle předpokladu λs 6= λi pro i = 1, 2, . . . , j, jsou všechna ci rovna nule. To ale znamená, že xs = 0, což je spor s definicí vlastního vektoru. Takže vektory x1 , . . . , xs jsou skutečně lineárně nezávislé. Má-li tedy operátor nad n-dimenzionálním prostorem n různých vlastních čísel, tvoří jeho vlastní kořeny bázi tohoto prostoru.
2.7
Konjugovaný operátor
Nechť A je operátor v n-rozměrném unitárním prostoru R a nechť A∗ je operátor nad stejným prostorem, pro který je splněna rovnost (Ax, y) = (x, A∗ y) pro všechna x, y ∈ R. Pak tento operátor nazveme konjugovaným operátorem k operátoru A. Ukážeme o tomto operátoru, nejen že vždy existuje, ale také že je určen jednoznačně. K tomuto účelu je ale nejprve nutné ukázat, že skalární součin má v souřadnicích libovolné ortonormální báze α = (u1 , u2 , . . . , un ) vyjádření (x, y) = x1 y¯1 + x2 y¯2 + · · · + xn y¯n , kde (x)α = (x1 , x2 , . . . , xn )′ a (y)α = (y1 , y2, . . . , yn )′ pro libovolné dva vektory x, y. To vidíme ihned z rovnosti n n n n X X X X (x, y) = ( xi ui , yj uj ) = xi y¯j (ui , uj ) = xi y¯i i=1
j=1
i,j=1
15
i=1
Nechť tedy A je libovolný operátor a11 a21 A = .. . an1
nad n-rozměrným prostorem R a a12 . . . a1n a22 . . . a2n .. .. . . an2 . . . ann
je jemu odpovídající matice v některé ortonormální bázi α = (u1 , . . . , un ). Nechť dále x, y ∈ R jsou dva libovolné vektory se souřadnicemi (x1 , x2 , . . . , xn )′ , (y1 , y2, . . . , yn )′ v bázi α. Pak a11 x1 + · · · + a1n xn a21 x1 + · · · + a2n xn (AX)α = .. . an1 x1 + · · · + ann xn a proto můžeme psát
(AX, Y) =a11 x1 y¯1 + · · · + a1n xn y¯1 + +a21 x1 y¯2 + · · · + a2n xn y¯2 + + ..................... + +an1 x1 y¯n + · · · + ann xn y¯n = =x1 (a11 y¯1 + a21 y¯2 + · · · + an1 y¯n ) + + ........................... + +xn (a1n y¯1 + a2n y¯2 + · · · + ann y¯n ) = =(X, A∗Y). Odtud vidíme, že za konjugovaný operátor lze zvolit operátor, který má v některé ortonormální bázi matici konjugovanou s maticí původního operátoru v téže bázi. Zbývá nám dokázat už jen jednoznačnost konjugovaného operátoru. Nechť A∗1 a A∗2 jsou dva operátory konjugované k operátoru A. Potom pro ně platí, že (Ax, y) = (x, A∗1 y) = (x, A∗2y). Odtud můžeme psát (x, (A∗1 − A∗2 )y) = 0. To ale znamená, že vektor (A∗1 − A∗2 )y je ortogonální ke každému vektoru x prostoru R, a tedy musí být nutně (A∗1 − A∗2 )y = 0 pro každý vektor y ∈ R. Tahle rovnost může být splněna jen v případě, že A∗1 − A∗2 = 0, tj. A∗1 = A∗2 . V následujícím textu vyšetříme vzájemné vztahy mezi vlastními čísly a vlastními vektory navzájem konjugovaných operátorů A a A∗ . Označíme-li A matici operátoru A v libovolné bázi, pak zřejmě platí det(A − λE) = det((A − λE)′ ) = det(A′ − λE) = ¯ = = det (A′ − λE) = det (A¯′ − λE) ¯ = det (A∗ − λE) . Odtud plyne, že číslo λ je kořenem charakteristického polynomu operátoru A právě ¯ kořenem charakteristického polynomu operátoru A∗ . když je λ 16
Označme dále xi vlastní vektor operátoru A nad n-rozměrným prostorem odpovídající vlastnímu číslu λi a yi vlastní vektor operátoru A∗ odpovídající vlastnímu ¯ j . Pak platí číslu λ (Axi , yj ) = (λi xi , yj ) = λi (xi , yj ) ¯ j yj ) = λj (xi , yj ) (Axi , yj ) = (xi , A∗yj ) = (xi , λ a tedy (λi − λj )(xi , yj ) = 0 . Pokud λi 6= λj , tj. λi − λj 6= 0, dostáváme (xi , yj ) = 0. Nechť nyní operátor A má n navzájem různých vlastních čísel. Pak pro λi = λj pokud by platilo (xi , yj ) = 0, tj. (xi , yi ) = 0, byl by vektor yi ortogonální ke všem vlastním vektorům x1 , . . . , xn operátoru A, a tedy i ke všem vektorům prostoru R. To by znamenalo, že yi je nulový vektor, což je spor s předpokladem, že yi je vlastní vektor. Tedy (xi , yj ) 6= 0. Vidíme tedy, že (xi , yj ) = 0 právě tehdy, když λi = λj . Ukažme ještě, že při libovolně vybraných vlastních vektorech y1 , . . . , yn , lze vektory x1 , . . . , xn normovat tak, že pro λi 6= λj je (xi yi ) = 1. Stačí místo nich vzít totiž vektory 1 1 x1 , . . . , xn , α1 αn kde αi = (xi , yi ) 6= 0, neboť
1 xi , yi αi
=
1 (xi , yi ) = 1 . αi
17
Kapitola 3 Základní přímé metody řešení problému vlastních čísel Úplným problémem vlastních čísel budeme rozumět úlohu určit všechna vlastní čísla dané čtvercové matice. V druhé kapitole jsme viděli, že koeficienty charakteristického polynomu pi jsou až na znaménko rovny součtům všech hlavních minorů i-tého řádu matice A. Přímý výpočet těchto koeficientů tedy vyžaduje při větším řádu matice obrovský počet operací. Proto bylo velmi přirozené vypracovat speciální numerické metody, které zjednodušují výpočet uvedené úlohy. V této kapitole se s některými z nich seznámíme. Poznamenejme ještě, že všechny zde probírané metody kromě metody Leverrierovy (z roku 1846) vznikly v třicátých letech minulého století nebo později. Při výkladu o numerických metodách budeme zpravidla předpokládat, že prvky matic jsou reálná čísla.
3.1
Metoda Leverrierova
Výhodou této metody je její jednoduchost a také to, že nemohou nastat výjimečné případy při výpočtu, jak tomu bude například u metody Krylovovy či Danilevského. Avšak cenou za tyto přednosti je právě větší náročnost při výpočtu. Samotná myšlenka spočívá v použití Newtonových vzorců. Je-li totiž ϕ(t) = (−1)n [tn − p1 tn−1 − p2 tn−2 − · · · − pn ] charakteristický polynom matice A s kořeny λ1 , λ2 , . . . , λn , potom označíme-li n X
λki = sk ,
i=1
platí Newtonovy vzorce:
p1 = s1 , kpk = sk − p1 sk−1 − · · · − pk−1s1 18
(k = 2, . . . , n) .
Odtud už lehce vyjádříme koeficienty charakteristického polynomu p1 = s1 , 1 p2 = (s2 − p1 s1 ) , 2 .. .
(3.1)
1 (sn − p1 sn−1 − · · · − pn−1 s1 ) . n Zbývá nám tedy ještě určit koeficienty sk . Zřejmě platí pn =
s1 = λ1 + λ2 + · · · + λn = Tr A . Jak již víme (viz. poznámka 2.4.1) jsou vlastní čísla matice Ak rovna právě číslům λk1 , λk2 , . . . , λkn a proto můžeme psát, že sk = Tr Ak . Vidíme tedy, že pro určení koeficientů charakteristického polynomu stačí postupně počítat mocniny dané matice až do n-tého řádu, jejich stopy a poté už můžeme tyto koeficienty vyjádřit jako řešení rekurentní soustavy rovnic (3.1). Přitom při počítání n-té mocniny stačí určit pouze diagonální prvky. Určování těchto mocnin je samo o sobě jednoduchou operací, avšak může být velmi pracné z hlediska počtu provedených operací. Celkově Leverrierova metoda vyžaduje 21 (n − 1)(2n3 − 2n2 + n + 2) operací násobení a dělení, což značně převyšuje počet stejného druhu operací při výpočtu metodou Danilevského a Krylovovou, a dokonce je i vyšší než u metody interpolační. Podotkněme ještě, že tahle metoda nijak neulehčuje výpočet vlastních vektorů, jak je tomu třeba u metody Krylovovy nebo Danilevského. Uveďme si pro ilustraci příklad (z [2]) výpočtu koeficientů charakteristického polynomu matice 1 2 3 4 2 1 2 3 A= 3 2 1 2 . 4 3 2 1 Výpočet bude probíhat následovně:
s1 = 4 , 30 22 A2 = 18 20 208 178 A3 = 192 242
22 18 16 18
18 16 18 22
20 18 , s2 = 96 , 22 30
178 148 154 192
192 154 148 178
242 192 , s3 = 712 , 178 208
19
2108 1388 A4 = 1388 p1 p2 p3 p4
3.2
, s4 = 6992 , 2108
=4, = (96 − (4).(4)′ )/2 = 40 , = (712 − (4 40).(96 4)′ )/3 = 56 , = (6992 − (4 40 56).(712 96 4)′ )/4 = 20 .
Metoda interpolační
Interpolační metody lze užít i na obecnější úlohy než jen na výpočet charakteristického polynomu. Umožňuje nám převést determinant tvaru f11 (t) . . . f1n (t) .. , (3.2) F (t) = ... . fn1 (t) . . . fnn (t) kde fik (t) jsou dané polynomy v t na polynomiální tvar. Protože polynom k-tého řádu lze s použitím různých interpolačních vzorců jednoznačně určit jeho hodnotami v k + 1 bodech, spočítáme nejprve k + 1 hodnot číselných determinantů F (ti ), kde k je předpokládaný stupeň charakteristického polynomu, pro nějaká obecně libovolně zvolená, čísla t0 , t1 , . . . , tk a poté sestrojíme tento polynom. K výpočtu je nejvýhodnější (viz. [1]) použít Newtonova interpolačního vzorce pro ekvidistantní argumenty ti (i = 0, 1, . . . , k). Uvedeme (podle [1]) Newtonův interpolační vzorec pro ti = i, i = 0, 1, . . . , k F (t) =
k X ∆i F (0) i=0
i!
t(t − 1) . . . (t − i + 1),
(3.3)
kde ∆i F (l) je i-tá diference vypočítaných hodnot polynomu F (t) definovaná rekurzivně takto: ∆0 F (l) = F (l), ∆i F (l) = ∆i−1 F (l + 1) − ∆i−1 F (l). Označme nyní i X 1 cmi tm pro 1 ≤ i ≤ k, t(t − 1) . . . (t − i + 1) = i! m=1
20
kde
0 X
cmi tm = 1 .
m=1
Potom můžeme (3.3) přepsat na tvar, jež je nazýván interpolačním vzorcem A. A. Markova: F (t) =
k X
∆i F (0)(
i X
cmi tm ) =
m=1
i=0 0
=∆ F (0) +
k X
i X
cmi ∆i F (0)tm =
i=1 m=1 =F (0) + (c11 ∆1 F (0)t)+ +(c12 ∆2 F (0)t + c22 ∆2 F (0)t2 )+ +(c13 ∆3 F (0)t + c23 ∆3 F (0)t2 + c33 ∆3 F (0)t3 )+
... +(c1k ∆k F (0)t + c2k ∆k F (0)t2 + · · · + ckk ∆k F (0)tk ) = k k X X =F (0) + ( cmi ∆i F (0))tm .
(3.4)
m=1 i=m
Vyšetřeme nyní koeficienty cmi . Napřed se pokusíme vyjádřit koeficienty polynomu t(t − 1)(t − 2) . . . (t − i + 1) a cmi dostaneme jejich podělením číslem i!. Označme (i)
(i)
(i)
t(t − 1)(t − 2) . . . (t − (i − 1)) = e0 ti − e1 ti−1 + · · · + (−1)i−1 ei−1 t1 . (i)
Zřejmě koeficient em bude součtem součinů všech m-tic různých prvků z {1, . . . , i−1} a proto můžeme psát (i)
e0 = 1, X (i) e1 = j, 1≤j≤i−1
(i) e2
=
X
j1 . j2 ,
1≤j1 <j2 ≤i−1
.. . (i)
ei−2 =
X
j1 . . . ji−2 ,
1≤j1 <···<ji−2 ≤i−1 (i)
ei−1 =
X
j1 . . . ji−1 = (i − 1)! .
1≤j1 <···<ji−1 ≤i−1
21
Pro výpočet těchto vztahů se dá použít rekurentního vzorce i X X X s j1 . . . jk = j1 . . . jk−1 , 1≤j1 <···<jk ≤i
1≤j1 <···<jk−1 ≤s−1
s=k
tedy
(i) ek
=
i−1 X
(j)
j . ek−1 .
j=k
Pro 0 < k < i − 1 můžeme výše uvedený vzorec ještě upravit (i) ek
=
i−1 X
j.
j=k
=
(i−1) ek
(j) ek−1
i−2 X (j) (i−1) j . ek−1 ) + (i − 1)ek−1 = =( j=k
(i−1)
+ (i − 1)ek−1 .
Při určování koeficientů interpolačního vzorce A. A. Markova se však zpravidla používá již předpočítaná tabulka koeficientů. Podívejme se nyní na počet provedených operací při výpočtu touto metodou. Nejprve je třeba vypočítat n+1 determinantů n-tého řádu, což samo o sobě vyžaduje 1 (n + 1)(n − 1)(n2 + n + 3) operací násobení a dělení. Poté, bereme-li koeficienty 3 Markovova interpolačního vzorce z tabulky, bude třeba ještě k výpočtu koeficientů provést 21 n(n + 1) operací násobení. Celkem tedy máme 1 1 (n + 1)(n − 1)(n2 + n + 3) + n(n + 1) , 3 2 operací násobení a dělení, což je sice nižší než u Leverrierovy metody, avšak pořád značně vyšší než u Danilevského nebo Krylovovy metody. Stejně jako v předešlém případě, nezjednodušuje nijak tato metoda výpočet vlastních vektorů. Jako příklad opět uveďme matici 1 2 3 4 2 1 2 3 A= 3 2 1 2 . 4 3 2 1 Potom tabulka determinantů a diferencí bude vypadat následovně: t 0 1 2 3 4 F (t) −20 −119 −308 −575 −884 ∆1 F (t) −99 −189 −267 −309 ∆2 F (t) −90 −78 −42 ∆3 F (t) 12 36 ∆4 F (t) 24 22
Dále doplňme koeficienty cmi Markovova interpolačního vzorce: i ∆i F (0) c4i c3i c2i c1i 1 −99 1, 00000 2 −90 0, 50000 −0, 50000 3 12 0, 16667 −0, 50000 0, 33333 4 24 0, 04167 −0, 25000 0, 45833 −0, 25000 Odtud už lehce spočítáme koeficienty p1 , p2 , p3 : p1 p2 p3 p4
3.3
= −c33 ∆3 F (0) + c34 ∆4 F (0) = 4 , = −c22 ∆2 F (0) + c23 ∆3 F (0) + c24 ∆4 F (0) = 40 , = −c11 ∆1 F (0)+ = c12 ∆2 F (0)+ = c13 ∆3 F (0)+ = c14 ∆4 F (0) = 56 , = −F (0) = 20 .
Metoda Krylovova
Nechť A je čtvercová matice n-tého řádu a11 − t a12 a21 a 22 − t ϕ(t) = .. .. . . an1 an2
a =0 . . . ann − t ... ...
a1n a2n .. .
(3.5)
její charakteristická rovnice. A. N. Krylov navrhl metodu založenou na myšlence, že se nejprve výše uvedená charakteristická rovnice transformuje na rovnici obecně s ní ekvivalentní tvaru b11 − t b12 . . . b1n b21 − t2 b22 . . . b2n D(t) = (3.6) .. .. .. = 0 , . . . bn1 − tn bn2 . . . bnn jejíž rozvoj podle mocniny t je pomocí Laplaceova rozvoje již značně jednodušší a odtud už bude možno snadno určit koeficienty charakteristického polynomu jako číselné násobky. Podívejme se tedy nejprve na podstatu Krylovovy transformace. Zřejmě rovnost (3.5) je nutnou a postačující podmínkou k tomu, aby soustava homogenních rovnic tx1 = a11 x1 + a12 x2 + · · · + a1n xn , tx2 = a21 x1 + a22 x2 + · · · + a2n xn , ................................. txn = an1 x1 + an2 x2 + · · · + ann xn , měla nenulové řešení x1 , x2 , . . . , xn . Násobme první rovnici číslem t t2 x1 = a11 tx1 + a12 tx2 + · · · + a1n txn 23
(3.7)
a dosaďme za tx1 , . . . , txn jejich vyjádření z (3.7) pomocí x1 , . . . , xn . Dostáváme tak t2 x1 =a11 (a11 x1 + a12 x2 + · · · + a1n xn )+ +a12 (a21 x1 + a22 x2 + · · · + a2n xn )+ .. . +a1n (an1 x1 + an2 x2 + · · · + ann xn ) = =(a11 a11 + a12 a21 + · · · + a1n an1 )x1 + +(a11 a12 + a12 a22 + · · · + a1n an2 )x2 + .. . +(a11 a1n + a12 a2n + · · · + a1n ann )xn a označíme-li b2k =
n X
a1s ask
(3.8)
s=1
můžeme psát t2 x1 = b21 x1 + b22 x2 + · · · + b2n xn .
(3.9)
Poté násobíme rovnici (3.9) znovu číslem t a opět dosadíme za tx1 , . . . , txn jejich vyjádřením pomocí x1 , . . . , xn . Tím dostaneme t3 x1 = b31 x1 + b32 x2 + · · · + b3n xn , kde b3k =
n X
b2s ask .
s=1
Je tedy zřejmé, že zopakujeme-li výše uvedený postup ještě (n − 3)-krát, dostaneme následující soustavu tx1 = b11 x1 + b12 x2 + · · · + b1n xn , t x2 = b21 x1 + b22 x2 + · · · + b2n xn , .. . tn xn = bn1 x1 + bn2 x2 + · · · + bnn xn , 2
(3.10)
jejíž koeficienty bik jsou určeny rekurentními vzorci b1k = a1k , n X bik = bi−1,s ask
(i = 2, . . . , n; k = 1, . . . , n).
(3.11)
s=1
Determinant soustavy (3.10) má zřejmě tvar (3.6). Přitom soustava (3.10) má nenulové řešení pro všechna t, pro která má nenulové řešení soustava (3.7), tj. pro 24
taková t, která vyhovují rovnici ϕ(t) hodnotách t, které jsou kořeny rovnice Ukážeme, že platí 1 0 b11 b12 D(t) = .. .. . . bn−1,1 bn−1,2
= 0. To znamená, že D(t) = 0 při všech ϕ(t) = 0. ϕ(t) = Nϕ(t), . . . bn−1,n
... ...
0 b1n .. .
(3.12)
tj., že pro N 6= 0 se D(t) liší od hledaného charakteristického polynomu jen číselným násobkem. Předpokládejme nejprve, že všechny kořeny polynomu ϕ(t) jsou navzájem různé. Protože tyto kořeny jsou zároveň kořeny D(t), je D(t) dělitelný ϕ(t). Avšak oba polynomy mají týž stupeň. To znamená, že polynom D(t) je násobkem polynomu ϕ(t), a přitom podíl N je konstanta (nezávislá na t). Snadno nahlédneme, že Laplaceovým rozvojem determinantu v (3.6) podle prvního sloupce dostaneme koeficient u tn ve tvaru b12 . . . b 1n b22 . . . b2n (−1)n .. .. . . . bn−1,2 . . . bn−1,n Koeficient polynomu ϕ(t) u tn je (−1)n . Odtud tedy platnost (3.12). Má-li polynom ϕ(t) vícenásobné kořeny, platí rovněž D(t) = Nϕ(t),
(3.13)
což lze ověřit přímo vynásobením determinantů, které se v této rovnosti vyskytují a užitím vztahů (3.11). Pro N = 0 z rovnosti (3.13) plyne, že D(t) se identicky rovná nule. Je tedy zřejmé, že v tomto případě uvedená transformace neumožňuje určení koeficientů charakteristického polynomu ϕ(t). I pro tento případ navrhl A. N. Krylov speciální metodu. Věnujme se ale nejprve vyšetření vlastností koeficientů bik určujících D(t). Za tímto účelem si definujeme vektory bi o souřadnicích bi1 , bi2 , . . . , bin . Z rovností (3.11) vyplývá, že bi−1,1 a11 a21 . . . an1 bi1 bi2 a12 a22 . . . an2 bi−1,2 .. = .. .. . .. , .. . . . . . bi−1,n a1n a2n . . . ann bin tedy
bi = A′ bi−1 . 25
(3.14)
A odtud již zřejmě bi = A′i b0 ,
(3.15)
neboť platí b1 = A′ b0 , kde b0 = (1, 0, . . . , 0)′ . Odtud už vidíme, že soustavu (3.7) je možné transformovat také tak, že vyjdeme např. od druhého sloupce této soustavy. Stejnými úvahami jako v předešlém textu bychom dospěli k determinantu D(t), ve kterém by se t vyskytovalo ve druhém sloupci a příslušné koeficienty bik by se určily podle vzorců (3.15), kde b0 = (0, 1, 0, . . . , 0)′ . Krylovovu metodu však lze, jak dále ukážeme, přirozeným způsobem zobecnit, budeme-li volit místo vektoru B0 speciálního typu libovolný vektor b0 = (b01 , b02 , . . . , b0n )′ . Je-li x1 , x2 , . . . , xn řešení soustavy (3.7), označme u = b01 x1 + b02 x2 + · · · + b0n xn . Násobením a dosazováním (stejně jako v předešlém případě) dostaneme následující soustavu rovnic u =b01 x1 + b02 x2 + · · · + b0n xn , ut =b11 x1 + b12 x2 + · · · + b1n xn , ut2 =b21 x1 + b22 x2 + · · · + b2n xn , .. . n ut =bn1 x1 + bn2 x2 + · · · + bnn xn ,
(3.16)
kde bi = (bi1 , bi2 , . . . , bin )′ = A′i b0 . Tato soustava n+1 lineárních homogenních rovnic o n + 1 neznámých u, x1 , . . . , xn má nenulové řešení právě když determinant 1 b01 . . . b0n t b11 . . . b1n D(t) = .. .. = 0 . .. . . n . t bn1 . . . bnn I zde bychom analogicky předešlým úvahám došli opět ke vztahu D(t) = Nϕ(t) , kde N je tentokrát b01 b . . . b 02 0n b11 b12 . . . b1n N = .. .. . .. . . . bn−1,1 bn−1,2 . . . bn−1,n
(3.17)
Pro N = 0, stejně jako v předešlém případě, uvedená transformace neumožňuje určení koeficientů polynomu ϕ(t). Nechť tedy nejprve N 6= 0. Označme Ni algebraické doplňky prvků tn−i v determinantu D(t). Pak z rovnosti D(t) = Nϕ(t) vyplývá, že koeficienty pi charakteristického polynomu jsou určeny poměrem (−1)n−1 Ni /N. 26
Podstatou práce A. N. Krylova je určení koeficientů charakteristického polynomu právě z těchto vztahů. Avšak uvedené úvahy umožňují najít hledané koeficienty bez výpočtu minorů, což podstatně snižuje náročnost výpočtu. Podmínka N 6= 0 je ekvivalentní lineární nezávislosti řádků determinantu (3.17), což jsou právě vektory b0 , b1 , . . . , bn−1 . Protože jich je n, tvoří bázi prostoru a proto můžeme zbylý vektor bn vyjádřit pomocí jejich lineární kombinace ve tvaru bn = q1 bn−1 + · · · + qn b0 .
(3.18)
Jak dále ukážeme, jsou koeficienty qi rovny koeficientům pi charakteristického polynomu ϕ(t) = (−1)n [tn − p1 tn−1 − · · · − pn ], což nám umožní určit tyto koeficienty jako řešení soustavy lineárních rovnic ekvivalentní této vektorové rovnosti. Odečtěme od posledního řádku determinantu D(t) lineární kombinaci předcházejících řádků s koeficienty q1 , q2 , . . . , qn . Tím dostaneme s využitím rovnosti (3.18), že 1 b01 . . . b0n t b11 . . . b1n .. .. .. = D(t) = . . . n−1 t bn−1,1 . . . bn−1,n n t − q1 tn−1 − · · · − qn 0 ... 0 =(−1)n [tn − q1 tn−1 − · · · − qn ] N ,
odkud už zřejmě ϕ(t) =
D(t) = (−1)n [tn − q1 tn−1 − · · · − qn ] . N
Koeficienty pi a qi jsou si tedy opravdu rovny. Podívejme se ještě, jak lze k témuž výsledku dojít pomocí Cayleyovy-Hamiltonovy věty. Aplikujeme-li ji na matici A′ , dostaneme A′n = p1 A′n−1 + · · · + pn E , a násobením b0 pak získáváme právě tvar A′n b0 = p1 A′n−1 b0 + · · · + pn b0 , což je bn = p1 bn−1 + · · · + pn b0 .
(3.19)
Aplikujeme-li však tuto větu na matici A, můžeme zřejmě místo soustavy (3.19) použít k určení koeficientů pi též soustavy cn = p1 cn−1 + · · · + pn c0 , kde ck = Ak c0 . 27
(3.20)
Určení koeficientů pi Krylovovou metodou s použitím soustavy (3.19) nebo (3.20) vyžaduje provést 32 n2 (n + 1) operací násobení a dělení, což je podstatně méně než u obou předcházejících metod. Pokud bychom určovaly tyto koeficienty A. N. Krylovovou metodou, tak jak ji původně navrhoval ve své práci, tj. prostřednictvím výpočtů minorů, vyžadovalo by to provést 13 (n4 +4n3 +2n2 −n−3) operací násobení a dělení. Pro ilustraci uveďme příklad opět pro matici 1 2 3 4 2 1 2 3 A= 3 2 1 2 4 3 2 1 a počáteční vektor c0 = (1, 0, 0, 0)′. Potom vektory c0 , c1 , . . . , c4 budou nabývat hodnot: c0 c1 c2 c3 c4 1 1 30 208 2108 22 178 1704 0 2 0 3 18 192 1656 0 4 20 242 1992 1 1 30 208 2108 0 1 11 89 852 0 0 −15 −75 −900 0 0 −24 −114 −1416 1 1 30 208 2108 0 1 11 89 852 1 5 60 0 0 0 0 0 6 24 1 1 30 208 2108 0 1 11 89 852 1 5 60 0 0 0 0 0 1 6
28
c0 c1 1 1 0 1 0 0 0 0 1 1 0 1 0 0 0 0 1 1 0 1 0 0 0 0 1 0 0 1 0 0 0 0
c2 c3 c4 30 208 2108 11 89 852 1 5 60 0 1 4 30 0 1276 11 0 496 1 0 40 0 1 4 0 0 76 0 0 56 1 0 40 0 1 4 0 0 20 0 0 56 1 0 40 0 1 4
Odsud vidíme, že p1 p2 p3 p4
=4 = 40 = 56 = 20 .
Vyšetřeme nyní situaci, kdy N = 0. V tomto případě soustava nehomogenních lineárních rovnic určená rovností (3.19) již neumožňuje určit koeficienty pi charakteristického polynomu, neboť determinant této soustavy je roven nule. Pro tento případ navrhl A. N. Krylov postup spočívající v určení koeficientů nenulového polynomu nejmenšího stupně ϑ(t) takového, že ϑ(A)c0 = 0, tj. koeficienty polynomu nejmenšího stupně anulujícího vektor c0 . Tento polynom bude roven minimálnímu polynomu matice A a tedy jeho kořeny, jak již bylo dokázáno v druhé kapitole, budou totožné se všemi kořeny charakteristického polynomu. Při nevhodné volbě vektoru c0 však můžeme dostat namísto minimálního polynomu některého jeho dělitele, čímž ztratíme část kořenů charakteristické rovnice. Je zřejmé, že dostaneme-li nějakou volbou počátečního vektoru c0 dělitele minimálního polynomu, potom musíme volit jiný počáteční vektor, který bude lineárně nezávislý s původním vektorem c0 . Není-li minimální polynom ψ(t) matice A totožný s jejím charakteristickým polynomem (až na nenulový faktor), je stupně menšího než n, a protože platí ψ(A)c0 , musí být vektory c0 , Ac0 , . . . , An−1 c0 lineárně závislé. Proto platí N = 0 při libovolné volbě vektoru c0 . Krylovova metoda tedy umožňuje pro N 6= 0 určit koeficienty charakteristického polynomu a pro N = 0 pouze některého jeho dělitele (obecně minimálního polynomu). Vztah N = 0 se projeví, řešíme-li soustavu (3.20) Gaussovou metodou tím, že se v části rovnic eliminují všechny koeficienty zároveň, takže tyto rovnice dávají 29
identitu 0 = 0. Tyto rovnice (označme jejich počet n−m) je nutno vynechat. Ve zbylé soustavě je třeba vynechat n−m posledních sloupců, počínaje sloupcem absolutních členů (tj. sloupcem souřadnic vektoru cn ). Poslední ze zbylých sloupců, složený ze souřadnic vektoru cm , se pak zvolí za absolutní člen nové soustavy. Řešení soustavy pak určí koeficienty lineární závislosti vektoru cm na vektorech c0 , c1 , . . . , cm−1 , tj. koeficienty polynomu nejmenšího stupně anulujícího vektor c0 . Pro názornost si to ukažme na příkladu s maticí 5 30 −48 B = 3 14 −24 . 3 15 −25 c0 1 0 0 1 0 0
c1 5 3 3 5 3 0
c2 c3 −29 125 −15 63 −15 63 −29 125 −15 63 0 0
Vyřešíme-li soustavu s vynecháním posledního sloupce, získáme koeficienty dělitele charakteristického polynomu: p3 + 5p2 = −29 , 3p2 = −15 , Odsud q2 = −5, q1 = −4 a tedy t2 + 5t + 4 .
3.4
Metoda Danilevského
Tato metoda spočívá v tom, že danou matici považujeme za matici operátoru v kanonické bázi (e1 , e2 , . . . , en ). Z Cayleyovy-Hamiltonovy věty plyne, že An e1 = p1 An−1 e1 + p2 An−2 e1 + · · · + pn e1 , kde p1 , p2 , . . . , pn jsou hledané koeficienty charakteristického polynomu. Jsou-li vektory e1 , Ae1 , . . . , An−1 e1 lineárně nezávislé, můžeme je považovat za bázi prostoru, kterou označíme α. Vyjádřeme nyní matici P uvažovaného operátoru v této bázi. Z první kapitoly již víme, že P = ( (Ae1 )α (AAe1 )α (AA2 e1 )α . . . (AAn−1 e1 )α ), kde (AAk e1 )α jsou sloupce souřadnic vektorů AAk e1 v bázi α. Pak je již zřejmé, že
30
matice P má tzv. Frobeniův tvar 0 1 0 P = .. . 0 0
. . . 0 pn . . . 0 pn−1 . . . 0 pn−2 .. . . . .. . . . 0 . . . 1 p2 0 . . . 0 p1
0 0 1 .. .
Nyní popíšeme přechod od báze e1 , e2 , . . . , en k bázi α. Postupně v každém z n − 1 kroků budeme přecházet od báze e1 , Ae1 , . . . , Ak−1 e1 , ek+1 , . . . , en k bázi e1 , Ae1 , . . . , Ak e1 , ek+2 , . . . , en , avšak za předpokladu, že všechny tyto soustavy vektorů jsou skutečně lineárně nezávislé. Případy degenerace, kdy nějaká ze soustav bude lineárně závislá, budeme vyšetřovat později. Matici, kterou dostaneme po (k−1)-ním kroku budeme značit A(k) . Tedy A = A(1) P = A(n) . Z výše řečeného je zřejmé, že sloupce matice A(k) jsou souřadnicemi vektorů Ae1 , A2 e1 , . . . ,Ak e1 , Aek+1 , . . . , Aen v bázi e1 , Ae1 , . . . , Ak−1 e1 , ek+1 , . . . , en . Prvních k −1 sloupců matice A(k) bude tedy totožných s týmiž sloupci Frobeniovy matice P . Samotný přechod se bude realizovat podle následujícího předpisu A(k+1) = Sk−1 A(k) Sk , kde Sk je matice transformace e1 ,Ae1 ,. . . ,Ak−1 e1 ,ek+1 ,. . . ,en . Je 1 0 0 1 Sk = ... ... 0 0 0 0
souřadnic v bázích e1 ,Ae1 ,. . . ,Ak e1 ,ek+2 ,. . . ,en a zřejmé, že bude mít tento tvar . . . s1,k+1 . . . 0 0 . . . s2,k+1 . . . 0 0 .. .. .. , . . . . . . sn−1,k+1 . . . 1 0 . . . sn,k+1 . . . 0 1
kde s1,k+1, s2,k+1 , . . . , sn,k+1 jsou souřadnice vektoru Ak e1 v bázi e1 , Ae1 , . . . , Ak−1 e1 , (k) ek+1 , . . . , en , tedy prvky aik v k-tém sloupci matice A(k) . Inverzní matice Sk−1 bude mít tvar s1,k+1 a1,k 1 0 . . . − sk+1,k+1 1 0 . . . − ak+1,k ... 0 0 ... 0 0 0 1 . . . − s2,k+1 . . . 0 0 0 1 . . . − a2,k . . . 0 0 sk+1,k+1 ak+1,k . . . . .. .. .. .. .. .. . . .. .. . . . . . . . . s a 0 0 . . . − k,k+1 . . . 0 0 0 0 . . . − k,k . . . 0 0 sk+1,k+1 ak+1,k Sk−1 = = . 1 1 . . . 0 0 0 0 . . . ak+1,k . . . 0 0 0 0 . . . sk+1,k+1 0 0 . . . − ssk+2,k+1 . . . 0 0 0 0 . . . − aak+2,k . . . 0 0 k+1,k+1 k+1,k .. .. .. .. .. .. .. .. .. .. . . . . . . . . . . sn,k+1 an,k 0 0 . . . − sk+1,k+1 . . . 0 1 0 0 . . . − ak+1,k . . . 0 1 31
Při výpočtu matice A(k+1) budeme nejprve počítat pomocnou matici B (k) B (k) = Sk−1 A(k) = 1 0 .. . 0 = 0 0 . .. 0
(k)
0 ... −
(k)
a1k (k)
ak+1,k (k)
a
1 . . . − (k)2k ak+1,k .. .. . . 0 ... −
(k) akk (k) ak+1,k
1 ak+1,k (k) ak+2,k
0 ...
0 . . . − (k) ak+1,k .. .. . . 0 ... −
(k)
(k) ak+1,1
(a11 − a1,k
(k)
ank (k)
ak+1,k
(k)
. . . 0 . . . 0 .. (k) (k) . a11 a12 . . . (k) (k) . . . 0 a21 a22 . . . . .. . . . . . 0 .. (k) (k) an1 an2 . . . . . . 0 .. . ... 1 (k)
(k) ak+1,2
), (a12 − a1,k
(k) a1n (k) a2n = .. . (k)
ann
(k)
(k)
(k) ak+1,n
), . . . , (a1n − a1,k
)
(k) (k) (k) ak+1,k ak+1,k ak+1,k (k) (k) (k) (k) (k) ak+1,2 (k) (k) ak+1,n (k) ak+1,1 (k) (a21 − a2,k a(k) ), (a22 − a2,k a(k) ), . . . , (a2n − a2,k a(k) ) k+1,k k+1,k k+1,k . . . .. .. .. . = (k) (k) (k) ak+1,1 ak+1,2 ak+1,n , , . . . , (k) (k) (k) ak+1,k ak+1,k ak+1,k . . . .. .. .. (k) (k) (k) a a a (k) (k) k+1,1 (k) (k) k+1,2 (k) (k) (an1 − an,k (k) ), (an2 − an,k (k) ), . . . , (ann − an,k k+1,n ) (k)
ak+1,k
ak+1,k
ak+1,k
Z výše řečeného, že prvních k − 1 sloupců matice A(k) je totožných s týmiž sloupci Frobeniovy matice P , je zřejmé, že se při těchto operacích prvních k − 1 sloupců nezmění. V k-tém sloupci dostaneme na (k + 1)-ním místě jedničku a na všech ostatních místech nuly. Takto vypočítanou matici B (k) budeme nyní zprava násobit maticí Sk 0 1 0 . = .. 0 .. . 0
0 0 1 .. . 0 .. .
(k)
. . . 0 b1,k+1 (k) . . . 0 b2,k+1 (k) . . . 0 b3,k+1 .. . . .. . . . (k) . . . 1 bk+1,k+1 .. .. . .
0 ... 0
(k)
bn,k+1
... ... ... ... ...
B (k) Sk = (k) b1,n (k) b2,n 1 (k) b3,n 0 .. .. . . . (k) bk+1,n 0 .. 0 . (k) bn,n 32
0 ... 1 ... .. .
s1,k+1 s2,k+1 .. .
0 ... 0 ...
sn−1,k+1 sn,k+1
... 0 0 ... 0 0 .. .. . .
. ... 1 0 ... 0 1
Tímto násobením se změní pouze (k + 1)-ní sloupec. Jeho prvky budou zřejmě čísla n X
(k) (k) b1j ajk ,
j=1
n X
(k) (k) b2j ajk ,
...,
j=1
n X
(k) (k)
bnj ajk .
j=1
Celkem tedy přechod od matice A(k) k matici A(k+1) bude realizován následující soustavou rekurentních vzorců (k)
(k) bk+1,j
=
ak+1,j (k)
ak+1,k
(k)
(k)
(k+1) aij
(k) bij n X
, (k) (k)
bi,j = ai,j − aik bk+1,j
(k+1)
=
ai,k+1 =
(i 6= k + 1), (j 6= k + 1),
(k) (k)
bij ajk .
j=1
Tato metoda je, z hlediska počtu operací násobení a dělení potřebných k výpočtu koeficientů charakteristického polynomu, nejlepší ze všech v této kapitole uvedených metod. Při výpočtu matice A(k) je třeba provést n − k + (n − k)(n − 1) + (n − k)n = 2n(n − k) těchto kroků. Celkem tedy máme n3 − n2 operací násobení a dělení. Jako příklad uveďme opět matici 1 2 3 4 2 1 2 3 A= 3 2 1 2 . 4 3 2 1 Výpočet v k-tém kroku budeme zapisovat do následující tabulky: (k)
bk+1,j A(k) (k+1) B (k) ai,k+1 A(k+1)
33
Tabulka celého výpočtu: 1, 000 1, 000 2, 000 3, 000 4, 000 2, 000 1, 000 2, 000 3, 000 0, 500 3, 000 2, 000 1, 000 2, 000 1, 000 4, 000 3, 000 2, 000 1, 000 1, 500 0, 000 1.500 2.000 2.500 19, 000 11, 000 1, 000 0.500 1.000 1.500 0, 000 0.500 −2.000 −2.500 −15, 000 0, 000 1.000 −2.000 −5.000 −24, 000 0, 000 19, 000 2, 000 2, 500 0, 000 1, 000 11, 000 1, 000 1, 500 1, 000 0, 000 −15, 000 −2, 000 −2, 500 1, 133 0, 000 −24, 000 −2, 000 −5, 000 1, 167 0, 000 0, 000 0, 533 −0.667 24, 000 1, 000 0, 000 −0.467 −0.333 34, 000 0, 000 1, 000 0, 133 0, 167 5, 000 6, 000 0, 000 0, 000 1, 200 −1, 000 0, 000 0, 000 24, 000 −0.667 0, 000 1, 000 0, 000 34, 000 −0, 333 0, 000 0, 000 1, 000 5, 000 0, 167 1, 000 0, 000 0, 000 6, 000 −1, 000 −0, 167 0, 000 0, 000 0, 000 3, 333 20, 000 1, 000 0, 000 0, 000 5, 333 56, 000 0, 000 1, 000 0, 000 1, 000 40, 000 0, 000 0, 000 1, 000 −0, 167 4, 000 0, 000 0, 000 0, 000 20, 000 1, 000 0, 000 0, 000 56, 000 0, 000 1, 000 0, 000 40, 000 0, 000 0, 000 1, 000 4, 000 Koeficienty pi charakteristického polynomu matice A budou tedy p1 p2 p3 p4
=4, = 56 , = 40 , = 20 .
Ukažme nyní, jak postupovat v případech degenerace. Může se stát, že při přechodu od báze e1 ,Ae1 ,. . . ,Ak−1 e1 ,ek+1 , . . . ,en k bázi e1 ,Ae1 ,. . . ,Ak e1 ,ek+2 ,. . . ,en , budou tyto vektory lineárně závislé. Protože e1 ,Ae1 ,. . . ,Ak−1 e1 ,ek+2 , . . . ,en jsou lineárně nezávislé, znamená to, že Ak e1 = q1 e1 + q2 Ae1 + · · · + qk Ak−1 e1 + qk+2 ek+2 + · · · + qn en . (k)
Zřejmě čísla qi jsou totožná s prvky matice A(k) v k-tém sloupci. A proto je ak+1,k = 0. V tomto případě už nemůžeme pokračovat obvyklou Danilevského metodou. Mohou nastat tyto dva případy: 34
(k)
1. Prvek ajk je různý od nuly pro nějaké j > k + 1. Potom v matici A(k) vyměníme (k + 1)-ní a j-tý řádek a zároveň (k + 1)-ní a j-tý sloupec. Tato záměna je ekvivalentní přechodu od báze e1 , Ae1 , . . . , Ak−1 e1 , ek+1, . . . , ej , . . . , en k bázi e1 , Ae1 , . . . , Ak−1e1 , ej , . . . , ek+1 , . . . , en . Snadno nahlédneme, že jsou-li vektory e1 , Ae1 , . . . , An−1 e1 lineárně nezávislé, pak nutně existuje takový index j. Po této transformaci můžeme pokračovat ve výpočtu. Hledáme-li pak vlastní vektory, musíme uvedenou transformaci brát ovšem v úvahu. Ve vhodném okamžiku musíme vyměnit (k + 1)-ní a j-tou souřadnici příslušných vektorů. (k) 2. Jsou-li všechna ajk pro j ≥ k+1 nulová, znamená to, že vektory e1 , Ae1 , . . . , Ak e1 jsou lineárně závislé, takže matice A(k) má tvar (k) 0 0 . . . 0 a1k (k) ! 1 0 . . . 0 a2k (k) . . A1 C .. .. .. .. . . . C (k) = 0 A2 (k) 0 0 . . . 1 akk (k) 0 A2 Charakteristický polynom matice A je tedy roven součinu charakteristických poly(k) (k) (k) nomů matice A1 a A2 . Přitom matice A1 už je ve Frobeniově kanonickém tvaru, takže tuto matici můžeme ihned rozvinout v polynom. K dokončení výpočtu charak(k) teristického polynomu nám tedy zbývá upravit metodou Danilevského matici A2 . To znamená, že se v tomto případě proces výpočtu charakteristického polynomu v každé etapě zjednoduší.
35
Kapitola 4 Stabilita problému vlastních čísel V této kapitole odpovíme na otázku, jak se změní vlastní čísla a vlastní vektory matice, změní-li se její prvky v mezích přípustné chyby a ukážeme, že v některých případech nemusí být problém vlastních čísel stabilní. Nechť A je čtvercová matice řádu n a A + dA je k ní blízká matice. Budeme předpokládat, že všechna vlastní čísla matice jsou navzájem různá a že prvky matice dA (a obdobně dx a dλ) jsou velmi malé. Máme-li Axi = λi xi pro i = 1, 2, . . . , n, dostáváme přibližně (dA) xi + A dxi = λi dxi + (dλi) xi . Uvažujme nyní matici A∗ a její vlastní vektory v1 , . . . , vn příslušné vlastním číslům λ¯1 , . . . , λ¯n . Pak skalárním násobením dostáváme ((dA) xi , vj ) + (A dxi , vj ) = (λi dxi , vj ) + ((dλi)xi , vj ) ((dA) xi , vj ) + (A dxi , vj ) = λi (dxi , vj ) + dλi(xi , vj ) , odkud pro i = j dostaneme ((dA) xi , vi ) + (A dxi , vi ) = (dxi , λ¯i vi ) + dλi (xi , vi ) a protože (xi , vi ) 6= 0 a A∗ vi = λ¯i vi , platí dλi =
((dA) xi , vi ) . (xi , vi )
(4.1)
Pro i 6= j je pak zřejmě (xi , vi ) = 0 a vzhledem k rovnosti (A dxi , vj ) = λj (dxi , vj ) dostáváme ((dA) xi , vj ) + λj (dxi , vj ) = λi (dxi , vj ), (λi − λj )(dxi , vj ) = ((dA) xi , vj ), ((dA) xi , vj ) . (dxi , vj ) = λi − λj Nechť dxi =
n X j=1
36
αij xj .
(4.2)
Potom zřejmě (dxi , vj ) = αij (xj , vj ) a tedy pro i 6= j αij =
((dA)xi , vj ) . (xj , vj )(λi − λj )
(4.3)
Koeficient αii je vzhledem k nejednoznačnosti vlastního vektoru xi neurčen. Bez újmy na obecnosti můžeme předpokládat, že αii = 0. Normováním ze vzorce (4.1) dostáváme následující vztah
dA .|xi |.|vi |
= ci dA , |dλi| ≤ |(xi , vi )|
kde
ci =
|xi |.|vi| . |(xi , vi )|
Z Cauchyho-Buňakovského nerovnosti je zřejmě ci ≥ 1. V případě reálných vlastních vektorů máme 1 ci = , | cos ϕi | kde ϕi je úhel mezi vektory xi a vi . Definice 5. Za výše uvedeného označení nazveme číslo ci koeficientem deformace matice A odpovídající vlastnímu číslu λi .
Koeficient deformace ci tedy vyjadřuje, že při dané normě dA může být změna λi tím větší, čím větší je tento koeficient deformace.
U hermitovských, specielně symetrických matic, platí ci = 1, tedy |dλi | ≤ dA , protože xi = vi a proto je v tomto případě úloha určit vlastní čísla vždy stabilní. Pro libovolnou matici platí, že je tato úloha nestabilní jen při velkém koeficientu deformace. Přejděme nyní ke stabilitě výpočtu vlastních vektorů. Z (4.2) a (4.3) plyne, že
dXi =
X
j∈{1,...,n}\{i}
((dA)Xi, Vj ) Xj (Xj , Vj )(λi − λj )
a odtud
|dXi| ≤ dA .|Xi |
= dA .|Xi |
X
j∈{1,...,n}\{i}
X
j∈{1,...,n}\{i}
|Vj |.|Xj | = |(Xj , Vj )|.|λi − λj | cj . |λi − λj |
Vidíme tedy, že úloha určit vlastní vektory matice může být nestabilní jen v případě, že alespoň jeden koeficient deformace ci je velký nebo jsou-li některá vlastní čísla příliš blízká. 37
Jako příklad (z [1]) uveďme nejprve 5 7 A= 6 5
a matici k ní blízkou
symetrickou matici 7 6 5 10 8 7 8 10 9 7 9 10
5, 1 7 6 5 7 10 8 7 , A + dA = 6 8 10 9 5 7 9 10
kde dA = 0, 1. Vlastní čísla první matice jsou po zaokrouhlení na čtyři desetiná místa 30, 2887 3, 8581 0, 8431 0, 0102 a vlastní čísla druhé matice jsou (opět po zaokrouhlení) 30, 3032 3, 8740 0, 8441 0, 0787. Nejvíce se změna projevilo u posledního vlastního čísla a to o 0, 0685. Vidíme tedy, že se opravdu změna vlastních čísel neprojevila více než o kolik se změnila norma matice A. Uveďme ještě jeden příklad tentokrát však na nesymetrickou matici. Stačí, když vezmeme matici 5 7 6 7 7 10 8 7 6 8 10 9 , 5 7 9 10
která se od předešlé liší pouze prvkem a15 . Budeme-li opět uvažovat změnu prvku a11 o 0, 1, budou vlastní čísla první matice po zaokrouhlení 30, 6789 3, 3601 0, 6454 0, 3157. Vlastní čísla druhé matice budou 30, 6943 3, 3658 0, 5555 0, 48448, což v odpovídá změně vlastních čísel až o 0, 1687. To opět souhlasí s našimi výsledky.
38
Kapitola 5 Příloha - kódy algoritmů Jak již bylo řečeno v úvodu, jednotlivé algoritmy pro výpočet koeficentů charakteristického polynomu jsou psány v programovacím jazyce matlab. U každého algoritmu je popsán jeho vstup a výstup a jsou uvedeny i pomocné funkce, jejichž specifikaci však zde už neuvádím. Ještě doplňme, že výpočet kořenů charakteristického polynomu je možno v matlabu provést pomocí jeho standardní funkce roots pro výpočet kořenů polynomu. Vstupním parametrem této funkce je právě vektor koeficientů charakteristického polynomu, jež je výstupem funkcí jednotlivých metod.
5.1
Kód algoritmu Leverrierovy metody
Vstup: čtvercová matice A řádu n. Výstup: vektor koeficientů p = (p1 , p2 , . . . , pn+1) takový, že p1 tn + p2 tn−1 + · · · + pn+1 je charakteristický polynom matice A.
function p = leverrier(A) [n,n]=size(A); B=A; s(1)=trace(B); for i=2:n-1, B=B*A; s(i)=trace(B); end s(n)=0; for i=1:n, for j=1:n, s(n)=s(n)+B(i,j)*A(j,i); end end 39
p(1)=s(1); for i=2:n, p(i)=(s(i)-p(1:i-1)*s(i-1:-1:1)’)/i; end p=[1 -p]*(-1)^n;
5.2
Kód algoritmu pro výpoče koeficientů cmi
Vstup: číslo n ∈ N. Výstup: tabulka (matice) C koeficintů cmi polynomu i X 1 t(t − 1) . . . (t − (i − 1)) = cmi tm , i! m=1
kde C(m, i) = cmi .
function C = markov(n) e=zeros(n+1,n+1); for i=1:n e(1,i+1)=1; for j=1:i-2 e(j+1,i+1)=e(j+1,i)+(i-1)*e(j,i); end if i>1 e(i,i+1)=(i-1)*e(i-1,i); end end c=zeros(n,n); for i=1:n for m=1:i C(m,i)=e(i-m+1,i+1)/factorial(i)*((-1)^(i-m)); end end
5.3
Kód algoritmu interpolační metody
Vstup: čtvercová matice A řádů n a čtvercová matice C koeficientů cmi o rozměru alespoň n × n. Výstup: vektor p = (p1 , . . . , pn+1 ) takový, že p1 tn + p2 tn−1 + · · · + pn+1 je charakteristický polynom matice A.
40
function p = interpolate(A,C) [n,n]=size(A); for i=0:n F(1,i+1)=det(A-i*eye(n,n)); end for i=2:n+1 for j=1:n-i+2 F(i,j)=F(i-1,j+1)-F(i-1,j); end end p(n+1)=F(1,1); for k=0:n-1 p(k+1)=C(n-k,[n-k:n])*F([n-k+1:n+1],1); end
5.4
Kód algoritmu Krylovovy metody
Vstup: čtvercová matice A řádu n a 1 ≤ r ≤ n udávající vektor kanonické báze er , který bude při výpočtu Krylovovy metody použit jako počáteční vektor. Výstup: vektor p = (p1 , . . . , pn+1 ) takový, že p1 tn + p2 tn−1 + · · · + pn+1 je charakteristický polynom matice A.
function p = krylov(A,r) [n,n]=size(A); c=vektoryC(A,r,n); x=zeros(1,n); for i=1:n x(i)=n-i+2; end p=zeros(1,n+1); [c,x,p]=primyChod(c,x,p); function c = vektoryC(A,r,n) c=zeros(n,n+1); c(1,r)=1; for i=2:n+1 for j=1:n for k=1:n c(j,i)=c(j,i)+A(j,k)*c(k,i-1); end end end function [c,x,p] = primyChod(c,x,p) 41
[n,m]=size(c); i=1; j=1; while i
while (j<=n) & (c(i,j)==0) j=j+1; end function c=vymenRadky(c,i,j) c=c1; pom=c(i,:); c(i,:)=c(j,:); c(j,:)=pom; function [c,x] = vymenSloupce(c,x,i,j) p1=c(:,i); c(:,i)=c(:,j); c(:,j)=p1; p2=x(i); x(i)=x(j); x(j)=p2;
5.5
Kód algoritmu Danilevslkého metody
Vstup: čtvercová matice A řádu n. Výstup: vektor p = (p1 , p2 , . . . , pn+1 ) takový, že p1 tn + p2 tn−1 + · · · + pn+1 je charakteristický polynom matice A.
function p = danilevskij(A) [n,n]=size(A); B=zeros(n,n); A1=zeros(n,n); bool=0; for k=1:n-1 if A(k+1,k)==0 j=nenulovyIndex(A,k,n); if j~=0 A=vymen(A,k,j); else p1=[1 (-1)*A([k:-1:1],k)]; p2=danielovskij(A([k+1:n],[k+1:n])); bool=1; break end end B=etapa1(A,k,n); A1=etapa2(A,B,k,n); A=A1; end if bool==1 p=vynasob(p1,p2); else p=[1 (-1)*A([n:-1:1],n)’]; 43
end function s = nenulovyIndex(A,k,n) s=0; for i=k+2:n if A(i,k)~=0 s=i; break; end end function A = vymen(A,k,j) pom=A(k+1,:); A(k+1,:)=A(j,:); A(j,:)=pom pom=A(:,k+1); A(:,k+1)=A(:,j); A(:,j)=pom; function p = vynasob(a,b) n=length(a); s=length(b); z=zeros(1,n+2*(s-1)); z([s:s+(n-1)])=a; p=zeros(1,s+n-1); for i=n+(s-1):-1:1 p(i)=z([i+(s-1):-1:i])*b’; end function B = etapa1(A,k,n) B=zeros(n,n); for j=1:n B(k+1,j)=A(k+1,j)/A(k+1,k); end for i=1:k for j=1:n B(i,j)=A(i,j)-A(i,k)*B(k+1,j); end end for i=k+2:n for j=1:n B(i,j)=A(i,j)-A(i,k)*B(k+1,j); end end function A1 = etapa2(A,B,k,n) for i=1:n for j=1:k A1(i,j)=B(i,j); end A1(i,k+1)=B(i,:)*A(:,k); for j=k+2:n A1(i,j)=B(i,j); end end
44
5.6
Kód algoritmu pro měření přesnosti jednotlivých metod
Vstup: číslo n ∈ N. Výstup: matice n×2, která má v prvním sloupci čísla 1, 2, . . . , n a v druhém sloupci maximální rozdíly vlastních čísel náhodných matic řádu 1, 2, . . . , n nalezené jednou z přímých metod (v tomto případě to bude metoda interpolační) a nepřímou metodou implementovanou standardní funkcí matlabu s názvem eig. Maximum se pro každý řád matice bere přes 20 různých náhodně vygenerovaných matic.
function b = compare(n) b=zeros(n,2); c=mar(n); for i=2:n x=zeros(1,5); for j=1:20 A=rand(i,i); funkceEig=sort(eig(A)); U=sort(roots(dan(A))); x(j)=max(abs(funkceEig-U)); end b(i,1)=i; b(i,2)=max(x); end
45
Závěr Závěrem porovnáme přesnosti jednotlivých metod s metodou nepřímou implementovanou standardní funkcí matlabu s názvem eig. Uvedeme tabulku, v jejímž prvním sloupci bude uveden řád matice a ve zbývajících sloupcích budou maximální rozdíli vlastích čísel matice daného řádu nalezených jednotlivými přímými metodami a nepřímou metodou. Maximum se bere přes dvacet náhodně vygenerovaných matic daného řádu a všechny vlastní čísla.
46
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
0 4, 4408e − 16 3, 5527e − 15 5, 7731e − 15 6, 4487e − 14 1, 4191e − 12 1, 1269e − 11 9, 0483e − 11 1, 6061e − 09 3, 4830e − 08 6, 8753e − 07 3, 9735e − 06 0, 0004 0, 0019 0, 0289 1, 3185 2, 1249 2, 4094 2, 4996 2, 8750 3, 1426 3, 1199 3, 4795 3, 9509 4, 1146 4, 5415 4, 6843 4, 8981 5, 3433 5, 7731
0 3, 3306e − 16 3, 9968e − 15 6, 4948e − 14 6, 2183e − 13 3, 7099e − 11 1, 2804e − 09 3, 1480e − 08 2, 1971e − 06 8, 4217e − 05 0, 0009 0, 1059 2, 0170 2, 7416 3, 2017 3, 7374 4, 9093 5, 8007 7, 2971 9, 4175 22, 1730 25, 7873 32, 8355 38, 8408 50, 6444 65, 9749 97, 1297 172, 0240 258, 2600 1867, 9096
0 2, 2204e − 16 2, 6645e − 15 8, 8762e − 14 9, 3458e − 13 2, 2421e − 11 4, 7829e − 11 7, 4972e − 10 1, 2331e − 08 1, 4315e − 07 3, 3273e − 06 9, 4285e − 05 0, 0024 0, 0657 1, 6584 40, 9188 24, 4693 151, 1945 97, 2332 156, 1146 54, 0166 17386, 0257 966, 7910 1169, 8838 2502, 5208 99, 5765 151, 5563 66, 6517 244, 7848 160, 1445
0 4, 4408e − 16 1, 3766e − 14 7, 3674e − 13 1, 6600e − 12 9, 1511e − 12 8, 2770e − 11 3, 4638e − 14 1, 8148e − 11 1, 4255e − 10 2, 9248e − 08 4, 7015e − 13 4, 7171e − 11 1, 0314e − 12 1, 2036e − 10 1, 0504e − 11 5, 0276e − 13 5, 4339e − 08 1, 1485e − 11 7, 4094e − 11 1, 0591e − 12 5, 3760e − 10 1, 1937e − 10 2, 4075e − 10 4, 4453e − 11 3, 5194e − 10 7, 7121e − 10 1, 0057e − 10 9, 0962e − 09 8, 8473e − 10
Odsud vidíme, že pro matice řádu větší než 10 jsou první tři metody prakticky nepoužitelné. Metoda Danilevského však vykazuje mnohem lepší výsledky a je prakticky, z hlediska přesnosti, použitelná i pro matice vyžšího řádu.
47
Tabulka použitých symolů
a ¯ A′ A¯ A∗ E det A adj A A−1 f ◦g dom(f ) (x)α (., .) Tr A ∆i f (t) |a| |x| A
komplexně sdružené číslo transponovaná matice komplexně sdružená matice konjugovaná matice A jednotková matice determinant matice A adjugovaná matice matice inverzní k matici A složení funkcí f a g definiční obor funkce f souřadnice vektoru x v bázi α skalární součin stopa matice A i-tá diference funkce f (t) absolutní hodnota čísla a norma vektoru x norma matice A
48
Literatura [1] Faddejev, Dmitrij Konstantinovič - Faddejeva, V. N. Numerické metody lineární algebry. 2. vyd. Praha : Státní nakladatelství technické literatury, 1964. [2] Demidovič, Boris Pavlovič - Maron, Isaak A. Základy numerické matematiky. 1. vyd. Praha : Státní nakladatelství technické literatury, 1966. [3] Lothar Collatz, Problémy charakteristických hodnot s technickými aplikacemi. 2. vyd. Praha : Státní nakladatelství technické literatury, 1965.
49
Rejstřík operátor, 9 identický, 10 konjugovaný, 15 nulový, 10 operátory lineární, 9
charakteristická rovnice, 6 charakteristické číslo, 7 charakteristický polynom, 6 charakteristický polynom operátoru, 14 hermitovsky symetrická matice, 6 hermitovská matice, 6 hlavní diagonála, 5
podobná matice, 8 polynom charakteristický, 6 operátoru, 14 minimální, 8 operátoru, 10
indentický operátor, 10 jednotková matice, 6 koeficient deformace, 37 koeficient deformace, 37 komplexně sdružená matice, 6 konjugovaná matice, 6 konjugovaný operátor, 15
rovnice charakteristická, 6 součet lineárních operátorů, 10 součin lineárních operátorů, 9
lineární operátor, 9 matice, 5 hermitovsky symetrická, 6 hermitovská, 6 jednotková, 6 komplexně sdružená, 6 konjugovaná, 6 normální, 6 nulová, 6 podobná, 8 transformace, 11 transponovaná, 6 unitární, 6 čtvercová, 5 minimální polynom, 8
transponovaná matice, 6 unitární matice, 6 vektor vlastní, 13 vlastní hodnota, 13 vlastní vektor, 13 vlastní číslo, 7, 13 čtvercová matice, 5 čísla charakteristická, 7 vlastní, 7, 13
normální matice, 6 nulová matice, 6 nulový operátor, 10 50
