Numerick´ e metody line´ arn´ı algebry
´ Uvod
1 1.1
´ Ulohy line´ arn´ı algebry
ˇ sen´ı soustav line´arn´ıch rovnic A~x = ~b 1. Reˇ ˇ sen´ı soustavy s regul´arn´ı ˇctvercovou matic´ı A ˇr´adu n×n pro jednu • Reˇ nebo v´ıce prav´ych stran • V´ypoˇcet inverzn´ı matice A−1 • V´ypoˇcet determinantu ˇ sen´ı soustav n rovnic o m nezn´am´ych a singul´arn´ıch n×n soustav • Reˇ v nˇejak´em definovan´em smyslu (1 ˇreˇsen´ı + b´aze nulprostoru, ˇreˇsen´ı s nejmenˇs´ı normou, ˇreˇsen´ı ve smyslu nejmenˇs´ıch ˇctverc˚ u) 2. Hled´an´ı vlastn´ıch ˇc´ısel a vlastn´ıch vektor˚ u ´ y probl´em vlastn´ıch ˇc´ısel • Upln´ ˇ asteˇcn´y probl´em vlastn´ıch ˇc´ısel • C´ Standardn´ı numerick´e knihovny - LINPACK (speci´aln´ı pro ˇreˇsen´ı soustav line´arn´ıch rovnic), EISPACK (speci´aln´ı pro vlastn´ı ˇc´ısla a vektory), NAG (obecn´a), IMSL (obecn´a)
1
1.2
Nˇ ekter´ e z´ akladn´ı pojmy a oznaˇ cen´ı
Vektorem zde budeme rozumˇet sloupcov´y vektor, horn´ı index T oznaˇcuje transponovan´y vektor (nebo matici), matici oznaˇcujeme tuˇcn´ym p´ısmem x1 a11 a12 . . . a1m x a a ... a 2m 2 21 22 ~x = .. = (x1, x2, . . . , xn)T , A = .. .. = (aij ) . .. .. . . . . . an1 an2 . . . anm xn
Pokud nebude uvedeno jinak, budeme uvaˇzovat re´aln´e matice a vektory. M´ırou velikosti vektor˚ u a matic je jejich vektorov´a norma. Vektorov´a norma mus´ı splnit 3 podm´ınky a) k~xk ≥ 0 a k~xk = 0 ⇔ ~x = ~0, b) kλ~xk = |λ|k~xk, c) k~x + ~y k ≤ k~xk + k~yk. Pˇr´ıklady vektorov´ych norem: • maximov´a: kxkI = max |xi | i=1,...,n
• souˇctov´a: kxkII =
n X
|xi|
i=1
• Eukleidovsk´a: kxkIII
v u n uX =t x2 i
i=1
Maticov´a norma Vektorov´a norma matice se naz´yv´a maticovou normou, pokud nav´ıc plat´ı pro vˇsechny matice A, B podm´ınka d) pro souˇcin matic kA · Bk ≤ kAk · kBk. Maticov´a norma je souhlasn´a s vektorovou normou, pokud ∀ A, ~x plat´ı kA~xk ≤ kAk · k~xk 2
Pˇr´ıklady maticov´ych norem: • ˇr´adkov´a: kAkI = max
i=1,...,n
n X j=1
• sloupcov´a: kAkII = max
j=1,...,n
• Eukleidovsk´a: kAkIII
|aij |
n X
|aij |
i=1
v uX n u n X =t a2ij i=1 j=1
Kaˇzd´a z uveden´ych maticov´ych norem je souhlasn´a se stejnˇe oznaˇcenou vektorovou normou. 1.3
Metody ˇ reˇ sen´ı soustav line´ arn´ıch rovnic
1. pˇr´ım´e (finitn´ı) 2. iteraˇcn´ı 3. gradientn´ı
3
2
Pˇ r´ım´ e metody ˇ reˇsen´ı soustav line´ arn´ıch rovnic
Pˇr´ım´e metody ˇreˇsen´ı spoˇc´ıvaj´ı v ´upravˇe matice na troj´uheln´ıkov´y (pˇr´ıpadnˇe diagon´aln´ı) tvar (pˇr´ım´y bˇeh) a potom ˇreˇsen´ım soustavy s horn´ı troj´uheln´ıkovou matic´ı U nebo doln´ı troj´uheln´ıkovou matic´ı L (zpˇetn´y bˇeh). Zpˇetn´y bˇeh je podstatnˇe rychlejˇs´ı neˇz pˇr´ım´y. 2.1
ˇ sen´ı soustav s troj´ Reˇ uheln´ıkovou matic´ı
Rovnice s horn´ı troj´uheln´ıkovou matic´ı U se ˇreˇs´ı postupn´ym prov´adˇen´ım vzorce ve smˇeru klesaj´ıc´ıho indexu k n X 1 ukj xj . xk = bk − ukk j=k+1
K v´ypoˇctu libovoln´eho xk je tˇreba nejv´yˇse n vnitˇrn´ıch cykl˚ u (1 n´asoben´ı + 1 2 2 sˇc´ıtan´ı), poˇcet operac´ı roste tedy ∼ n (pˇresnˇeji ≃ 0.5 n vnitˇrn´ıch cykl˚ u). 2.2
Gaussova a Gauss-Jordanova eliminace
ˇ s´ım soustavu rovnic A~x = ~b. V prvn´ım kroku necht’ prvek a11 6= 0 (ˇcehoˇz Reˇ lze vˇzdy dos´ahnout pˇrehozen´ım rovnic). Prvek a11, pouˇzit´y k u ´pravˇe rovnic 2, . . . , n nazveme hlavn´ım prvkem (pivot). (1)
Od i-t´e rovnice odeˇcteme 1. rovnici n´asobenou multiplik´atorem mi = −ai1 /a11. ´ Modifikovan´a soustava bude m´ıt v 1. sloupci pod diagon´alou sam´e 0. Uprava prov´adˇen´a souˇcasnˇe s pravou stranou odpov´ıd´a n´asoben´ı rovnice matic´ı 1 0 ... 0 −a /a 1 . . . 0 21 11 D1 = .. . . . .. . .. . . . −an1 /a11 0 . . . 1 Rovnice po prvn´ı ´upravˇe m´a tvar D1 A~x = D1~b. Oznaˇc´ıme A(1) ≡ D1 A ~b(1) ≡ D1~b.
4
a
Po k − 1 ´uprav´ach m´a matice A(k−1) tvar a11 a12 . . . a1,k−1 a1k (1) 0 a(1) . . . a(1) a2k 22 2,k−1 .. .. . . . .. .. . . . . (k−2) (k−2) 0 0 . . . ak−1,k−1 ak−1,k (k−1) A = (k−1) 0 0 ... 0 akk (k−1) 0 0 ... 0 ak+1,k .. .. .. .. .. . . . . . (k−1) 0 0 ... 0 ank
... ... .. . ... ... ... .. . ...
a1n (1) a2n .. . (k−2) ak−1,n (k−1) ak,n (k−1) ak+1,n .. . (k−1) an,n
.
Zde horn´ı index znaˇc´ı poˇcet ´uprav dan´eho prvku. (k−1) (k) Pokud akk 6= 0, lze ho zvolit za hlavn´ı prvek, spoˇc´ıtat multiplik´atory mi = (k−1) (k−1) −aik /akk pro i = k + 1, . . . , n a upravit pˇr´ısluˇsn´e rovnice. V k-t´em kroku ´upravy pouˇz´ıv´am jako hlavn´ı prvek prvek k − 1-kr´at upraven´y (odeˇc´ıt´an´ı!) hlavn´ı prvek → ztr´ata pˇresnosti ⇒ v´ybˇer hlavn´ıho prvku. Bez v´ybˇeru hlavn´ıho prvku jsou pˇr´ım´e metody nepouˇziteln´e pro obecn´e matice! Poˇcet operac´ı Na kaˇzdou nulu ≤ n vnitˇrn´ıch cykl˚ u, potˇrebuji vynulovat 12 n(n − 1) prvk˚ u. 3 3 Celkov´y poˇcet vnitˇrn´ıch cykl˚ u ∼ n (pˇresnˇeji ≃ 1/3 n ), sloˇzitost algoritmu je 3 ˇr´adu n . Gauss-Jordanova eliminace Upravuj´ı se vˇsechny prvky mimo diagon´alu. Matice se pˇrevede na jednotkovou I. Pˇr´ımo spoˇctu inverzn´ı matici A−1. Vyˇsˇs´ı poˇcet operac´ı ≃ n3 vnitˇrn´ıch cykl˚ u.
5
2.3
V´ ybˇ er hlavn´ıho prvku (pivoting)
V kaˇzd´em kroku pˇr´ım´eho bˇehu vyb´ır´ame hlavn´ı prvek. Z moˇzn´ych prvk˚ u vybereme nejvˇetˇs´ı v absolutn´ı hodnotˇe. Mal´e ˇc´ıslo vzniklo pravdˇepodobnˇe odeˇc´ıt´an´ım pˇri ´uprav´ach, je tedy pravdˇepodobnˇe zat´ıˇzeno velkou relativn´ı chybou a proto nen´ı jako pivot vhodn´e. ´ y v´ybˇer hlavn´ıho prvku: v cel´e dosud neupraven´e ˇc´asti matice hled´ame • Upln´ max |aij |. Pomal´y. ˇ asteˇcn´y v´ybˇer hlavn´ıho prvku: v dan´em sloupci (sloupcov´y) nebo ˇr´adku • C´ (ˇr´adkov´y). • Implicitn´ı v´ybˇer hlavn´ıho prvku: Rychlejˇs´ı, vylepˇsen´a strategie sloupcov´eho v´ybˇeru. Pˇri v´ybˇeru porovn´av´am velikosti prvk˚ u v dan´em sloupci normovan´e na maximum absolutn´ıch hodnot prvk˚ u v dan´em ˇr´adku p˚ uvodn´ı matice. V´ybˇer hlavn´ıho prvku → pouˇzitelnost pˇr´ım´ych metod pro vˇetˇsinu matic. Pro obecn´e velk´e matice (N > 50) nutn´a dvojit´a pˇresnost. I tak ˇcasto probl´emy u velk´ych ˇspatnˇe podm´ınˇen´ych matic! Probl´emy: 1. Singul´arn´ı matice 2. Singul´arn´ı matice vznikl´a ztr´atou pˇresnosti pˇri ´uprav´ach 3. Ztr´ata pˇresnosti
6
2.4
LU metoda
Kaˇzdou regul´arn´ı matici A lze rozloˇzit do tvaru A = L · U, kde L, U jsou lev´a doln´ı, resp. prav´a horn´ı troj´uheln´ıkov´e matice. Potom ˇreˇsen´ı najdu postupn´ym ˇreˇsen´ım 2 soustav s troj´uheln´ıkovou matic´ı A~x = ~b ⇒ (LU)~x = ~b ⇒ L (U~x) = ~b ⇒ L~y = ~b, U~x = ~y . | {z } ~y
LU dekompozice a11 a12 a13 a a a 21 22 23 a a a 31 32 33 .. .. .. . . . an1 an2 an3
... ... ... .. . ...
a1n a2n a3n .. . ann
=
1 0 0 l21 1 0 l31 l32 1 .. .. .. . . . ln1 ln2 ln3
... ... ... ... ...
0 0 0 .. . 1
u11 u12 u13 0 u22 u23 0 0 u33 .. .. .. . . . 0 0 0
... ... ... ... ...
N´asoben´ı matic pro i ≤ j :
aij = uij +
i−1 X
lik ukj
k=1 j−1
pro i > j :
aij = lij ujj +
X
lik ukj
k=1
Crout˚ uv algoritmus - postupn´y v´ypoˇcet napˇr. odleva po sloupc´ıch a ve sloupc´ıch odshora. Nejdˇr´ıve uij = aij −
i−1 X
lik ukj
i = 1, . . . , j
k=1
uˇz´ıv´a l z pˇredchoz´ıch sloupc˚ u a u z pˇredchoz´ıch ˇr´adk˚ u, a potom ! j−1 X lik ukj /ujj i = j + 1, . . . , n lij = aij − k=1
uˇz´ıv´a l z pˇredchoz´ıch sloupc˚ u a u z naddiagon´aln´ı ˇc´asti sloupce. Sloupcov´e hled´an´ı hlavn´ıho prvku (´upln´e nelze). Prvky aij pouˇziji jen jednou, v´ysledn´e prvky matic L a U se vejdou do jedn´e matice. 7
u1n u2n u3n .. . unn
Vlastnosti LU metody: • Pˇr´ım´a (finitn´ı) metoda, stejnˇe krok˚ u jako pˇr´ım´y bˇeh Gaussovy eliminace • Hlavn´ı v´yhoda: pˇri dekompozici nepracuji s pravou stranou rovnice, rychl´e v´ypoˇcty pro postupnˇe z´ısk´avan´e prav´e strany • Lze iterativnˇe zpˇresnit v´ysledek 2.5
Iterativn´ı zpˇ resnˇ en´ı ˇ reˇ sen´ı
Hled´ame ˇreˇsen´ı ~x line´arn´ı rovnice A~x = ~b. Z´ısk´ame nepˇresn´e ˇreˇsen´ı ~x˜ ~ ⇒ A(~x + δx) ~ = ~b + δb ~ ⇒ Aδx ~ = δb ~ = A~x˜ − ~b ~x˜ = ~x + δx ~ ~x = ~x˜ − δx Oznaˇc´ıme ~x0 nepˇresn´e ˇreˇsen´ı syst´emu v prvn´ım kroku A~x0 ≃ ~b a pak prov´ad´ıme iteraci ~ i = ~b − A~xi, A(δx)
~ i. ~xi+1 = ~xi + (δx)
8
2.6
Podm´ınˇ enost ˇ reˇ sen´ı soustavy line´ arn´ıch rovnic
Vzhledem k chyb´am vstupn´ıch ´udaj˚ u a zaokrouhlen´ı m´ısto A~x = ~b ˇreˇs´ıme ve skuteˇcnosti ´ulohu (A + ∆A) (~x + ∆~x) = ~b + ∆~b. • Nejdˇr´ıve pˇr´ıpad ∆A = 0: ∆~x = A−1∆~b ⇒ k∆~xk ≤ kA−1kk∆~bk, k~bk A~x = ~b ⇒ k~xk ≥ . kAk Pro relativn´ı chybu ˇreˇsen´ı tedy plat´ı ~ k∆~xk −1 k∆bk . ≤ kAkkA k k~xk k~bk ˇ ıslo Cp = kAkkA−1k se naz´yv´a podm´ınˇenost matice. C´ • Pokud je i ∆A 6= 0, pak k∆~bk k~bk . Cp k∆Ak kAk
k∆Ak kAk
k∆~xk ≤ Cp k~xk 1−
+
Pro Cp ≫ 1 je soustava ˇspatnˇe podm´ınˇen´a a mal´e chyby vstupn´ıch dat i mal´e zaokrouhlovac´ı chyby ve v´ypoˇctu se projev´ı velkou chybou ˇreˇsen´ı.
9
2.7
V´ ypoˇ cet inverzn´ı matice a determinantu
Gauss-Jordanova eliminace vypoˇcte pˇr´ımo inverzn´ı matici. U Gaussovy eliminace a LU dekompozice z´ısk´ame inverzn´ı matici ˇreˇsen´ım pro n prav´ych stran, tvoˇren´ych vektory standardn´ı b´aze. Vˇsechny tˇri metody jsou rovnocenn´e jak pˇresnost´ı, tak i poˇctem n3 operac´ı. Determinant nepoˇc´ıt´ame ze vzorc˚ u (r˚ ust zaokrouhlovac´ı chyby), ale vyuˇzijeme vˇety, ˇze determinant souˇcinu matic je roven souˇcinu determinant˚ u. Pro LU dekompozici det(A) = det(L) · det(U) =
n Y
ujj .
j=1
Gaussova eliminace odpov´ıd´a n´asoben´ı maticemi a pokud se s´am ˇr´adek nen´asob´ı a pouze se pˇriˇc´ıtaj´ı n´asobky jin´ych ˇr´adk˚ u, je determinant matice kaˇzd´e u ´pravy Di = 1 a determinant A se rovn´a souˇcinu diagon´aln´ıch prvk˚ u troj´uheln´ıkov´e matice. Pokud doch´az´ı k pˇrehazov´an´ı ˇr´adk˚ u pˇri v´ybˇeru hlavn´ıch prvk˚ u, dojde k n´asoben´ı determinantu −1, je nutno si zapamatovat poˇcet zmˇen znam´enka.
10
2.8
Speci´ aln´ı typy matic
ˇ ıdk´a matice m´a vˇetˇsinu prvk˚ • R´ u = 0. Pro ˇreˇsen´ı soustav s ˇr´ıdkou matic´ı se ˇcasto pouˇz´ıvaj´ı gradientn´ı metody, spoˇc´ıvaj´ıc´ı v minimalizaci vhodn´eho kvadratick´e formy, napˇr. kA~x −~bk2III . Pro ˇr´ıdkou matici je totiˇz poˇcet operac´ı pro v´ypoˇcet A~x ∼ n, a ne ∼ n2 jako pro plnou matici. • Matice A je p´asov´a, pokud aij = 0 pro |i − j| > p. Tridiagon´aln´ı matice pro p = 1, pˇetidiagon´aln´ı matice pro p = 2. Pro p´asov´e matice se pouˇz´ıvaj´ı pˇr´ım´e metody, pro ostatn´ı ˇr´ıdk´e matice jsou pˇr´ım´e metody vˇetˇsinou neefektivn´ı. • Soustavy a1 c2 0 . .. .. . . .. 0 0
s tridiagon´aln´ı matic´ı b1 0 0 . . . 0 0 0 a2 b2 0 . . . 0 0 0 c3 a3 b3 . . . 0 0 0 .. . . . . . . . . . .. .. .. . . . . .. .. . . . . . . . . . .. .. . . . . .. .. . . . . . . . . . .. .. . . . . 0 0 0 . . . cn−1 an−1 bn−1 0 0 0 ... 0 cn an
·
x1 x2 x3 .. . .. . .. . xn−1 xn
=
f1 f2 f3 .. . .. . .. . fn−1 fn
Tridiagon´aln´ı matici zap´ıˇseme do 3 vektor˚ u ~a, ~b, ~c. V praxi t´emˇeˇr vˇzdy tridiagon´aln´ı matice, u kter´ych v´ybˇer hlavn´ıho prvku nen´ı potˇrebn´y (silnˇe regul´arn´ı matice). ˇ sen´ı: Pˇredpokl´ad´ame zpˇetn´y bˇeh xk = µk xk+1 + ρk . Dosad´ıme Reˇ ci (µi−1 xi + ρi−1) + ai xi + bixi+1 = fi, po ´upravˇe xi =
fi − ci ρi−1 −bi xi+1 + , ci µi−1 + ai ci µi−1 + ai
v´ysledek µi =
−bi , ci µi−1 + ai
ρi =
fi − ci ρi−1 . ci µi−1 + ai
Startov´an´ı c1 = 0, bn = 0, {µ0, ρ0 , xn+1} libovoln´e. 11
• Blokovˇe tridiagon´aln´ı matice: Ai, Bi , Ci jsou mal´e matice ⇒ µi , ρi jsou mal´e matice 2.9
´ Ulohy se ˇ z´ adn´ ym ˇ reˇ sen´ım nebo ∞ ˇ reˇ sen´ımi
m rovnic o n nezn´am´ych, line´arnˇe z´avisl´e n × n syst´emy. Metoda SVD (singular value decomposition): pokud A~x = ~b m´a ∞ ˇreˇsen´ı, urˇc´ı ˇreˇsen´ı s nejmenˇs´ı Eukleidovskou normou a b´azi nulprostoru, pokud ˇreˇsen´ı neexistuje, najde ˇreˇsen´ı ve smyslu nejmenˇs´ıch ˇctverc˚ u - vektor ~x minimalizuj´ıc´ı kA~x − ~bkIII . SVD: A, U matice m × n, W, V, I matice n × n, I jednotkov´a, U, V ortogon´aln´ı (UT · U = VT · V = I) A = U · W · VT
⇒
W diagon´aln´ı,
~x = V · [diag(1/wj )] · UT · ~b
Pokud wj = 0 (wj ≃ 0) nahrad´ıme
1 wj
→ 0 (signalizuje singul´arnost matice).
12
2.10
Rychlost ˇ reˇ sen´ı soustav line´ arn´ıch rovnic (matice n×n) pˇ r´ım´ ymi metodami
• Gauss–Jordanova eliminace: Pˇri v´ypoˇctu se prov´ad´ı ∼ n3 vnitˇrn´ıch cykl˚ u (v kaˇzd´em cyklu jedno n´asoben´ı a jedno sˇc´ıt´an´ı) • Gaussova eliminace a LU dekompozice: Pˇr´ım´y bˇeh je podstatnˇe n´aroˇcnˇejˇs´ı na poˇcet operac´ı. U obou metod je k v´ypoˇctu pˇr´ım´eho bˇehu potˇreba ∼ 13 n3 vnitˇrn´ıch cykl˚ u. Pˇri zpˇetn´em bˇehu 1 ∼ 2 n(n − 1) v Gaussovˇe eliminaci. Pˇri LU dekompozici jsou zapotˇreb´ı 2 zpˇetnˇe bˇehy, ale pˇr´ım´y bˇeh je nepatrnˇe kratˇs´ı, protoˇze se neupravuje prav´a strana, operac´ı je u Gaussovy eliminace i LU dekompozice stejnˇe. Pro v´ypoˇcet inverzn´ı matice jsou pracnosti Gauss–Jordanovy eliminace, Gaussovy eliminace i LU dekompozice stejn´e. • Metoda ˇr´adu < 3: Byla dok´az´ana (Strassen) existence metody, kde poˇcet operac´ı ∼ N log2 7 (log2 7 ≃ 2.807) a tedy roste s dimenz´ı n matice pomaleji neˇz u klasick´ych metod, kde poˇcet operac´ı ∼ n3 , Tato metoda vyˇzaduje komplikovanou pr˚ ubˇeˇznou archivaci napoˇc´ıtan´ych hodnot. Pro mal´e matice je tato metoda podstatnˇe pomalejˇs´ı neˇz klasick´e metody a jej´ı v´yhody se projev´ı aˇz matice ˇr´adu n ≫ 1000.
13
3
Gradientn´ı metody
Line´arn´ı rovnici A~x = ~b ˇreˇs´ıme napˇr´ıklad minimalizac´ı funkce 1 f (~x) = |A~x − ~b|2 . 2 V kaˇzd´em kroku hled´ame λ takov´e, aby f (~x + λ~u) bylo minim´aln´ı. Tedy λ=
−~u ∇f , kde |A~u|2
∇f (~x) = AT (A~x − ~b).
Pro ˇr´ıdk´e matice se sloˇzitost n´asoben´ı vektoru matic´ı sniˇzuje z poˇctu operac´ı ∼ n2 na poˇcet operac´ı ∼ n. Pozn. Nejpouˇz´ıvanˇejˇs´ı gradientn´ı metodou je metoda sdruˇzen´ych (konjugovan´ych) gradient˚ u, kterou si uk´aˇzeme v kapitole vˇenovan´e hled´an´ı extr´em˚ u.
14
4
Iteraˇ cn´ı metody ˇ reˇsen´ı soustav line´ arn´ıch rovnic
4.1
Nˇ ekter´ e typy matic
• Matice A typu n × n je diagon´ alnˇ e dominantn´ı, pokud |aii | >
n X
|aij |
∀i = 1, 2, . . . , n.
j=1,j6=i
ˇ • Ctvercov´ a matice A je pozitivnˇ e definitn´ı, pokud ∀~x 6= ~0 je skal´arn´ı souˇcin (~x, A~x) > 0. Symetrick´a matice je pozitivnˇe definitn´ı pr´avˇe tehdy, kdyˇz m´a vˇsechna vlastn´ı ˇc´ısla kladn´a. 4.2
Iteraˇ cn´ı proces
ˇ sen´ı ~x line´arn´ı rovnice A~x = ~b odhadneme vektorem ~x(0) a dalˇs´ı pˇribl´ıˇzen´ı Reˇ k pˇresn´emu ˇreˇsen´ı vypoˇcteme pomoc´ı pˇredpisu ~x(k+1) = Bk ~x(k) + ~ck . Pro ˇreˇsen´ı ~x mus´ı platit ~x = Bk ~x + ~ck . Odtud vypl´yv´a ~x(k+1) − ~x = Bk (~x(k) − ~x) = Bk Bk−1(~x(k−1) − ~x) = . . . . Pro konvergenci iteraˇcn´ıho procesu je tedy nutn´e a staˇc´ı, aby lim Bk Bk−1 . . . B0 = 0.
k→∞
Iteraˇcn´ı metody dˇel´ıme na metody stacion´arn´ı, kter´e maj´ı matici B konstantn´ı (Bk = B), a na metody nestacion´arn´ı. 4.3
Pˇ r´ıklad nestacion´ arn´ı iteraˇ cn´ı metody
Pˇr´ıkladem nestacion´arn´ı iteraˇcn´ı metody je ˇr´ızen´a relaxace. Pro jednoduchost ji uk´aˇzeme pro matici A, kter´a m´a diagon´aln´ı ˇcleny jednotkov´e (aii = 1). Necht’ po k-t´e iteraci je ze sloˇzek rezidua |A~x − ~b| maxim´aln´ı i-t´a sloˇzka. Budeme tedy nulovat i-tou sloˇzku rezidua k + 1-n´ı iterac´ı (k+1)
xi
(k)
(k)
(k)
= bi − ai1 x1 − . . . − aii−1xi−1 − aii+1xi+1 − . . . − ain xn(k) . 15
Matice Bk a vektor ~ck maj´ı tvar 1 ... 1 Bk = −ai1 . . . −aii−1 0 −aii+1 . . . −ain 1 ... 1
,
~ck =
0 .. . 0 bi 0 .. . 0
.
Tato metoda se vˇsak nehod´ı pro vyuˇzit´ı na poˇc´ıtaˇci, protoˇze vyˇzaduje v kaˇzd´em kroku hled´an´ı rovnice, kde odchylka od ˇreˇsen´ı je maxim´aln´ı, a to je ˇcasovˇe n´aroˇcn´e. 4.4
Stacion´ arn´ı iteraˇ cn´ı metody
Pro vˇsechna vlastn´ı ˇc´ısla λ matice B (tedy ˇc´ısla, pro kter´ a existuje takov´y ~v , ˇze plat´ı B~v = λ~v ) mus´ı platit |λ| < 1. Vˇ eta Nutnou a postaˇcuj´ıc´ı podm´ınkou konvergence metody je podm´ınka, aby spektr´aln´ı polomˇer ̺(B) byl ̺(B) = max |λi | < 1. i=1,...,n
Vˇ eta Pokud v nˇekter´e maticov´e normˇe plat´ı kBk < 1, pak iterace konverguje. Odhad chyby. Chceme dos´ahnout pˇresnosti ε a tedy iterovat, dokud nebude k~x(k) − ~xk ≤ ε. V´yraz k~x(k) − ~xk lze odhadnout k~x(k) − ~xk = k~x(k) − A−1~bk = kA−1(A~x(k) − ~b)k ≤ kA−1k kA~x(k) − ~bk.
16
4.5
Prost´ a iterace
Line´arn´ı rovnici A~x = ~b pˇrevedeme na tvar ~x = (I − A)~x + ~b, kde I je jednotkov´a matice. Prost´a iterace je tedy d´ana iteraˇcn´ım vzorcem ~x(k+1) = (I − A)~x(k) + ~b. Oznaˇc´ıme B = I − A, pak lze pˇresnost k-t´e iterace odhadnout v´yrazem " # ~ kbk k~x(k) − ~xk ≤ kBkk k~x(0) k + . 1 − kBk Metoda prost´e iterace se pro syst´emy line´arn´ıch rovnic prakticky nepouˇz´ıv´a. 4.6
Jacobiho metoda
Pˇredpokladem metody je, ˇze matice A m´a nenulov´e diagon´aln´ı prvky aii 6= 0. Sloˇzky k + 1 Jacobiho iterace ˇreˇsen´ı jsou d´any aii+1 (k) ain (k) bi aii−1 (k) ai1 (k) x1 − . . . − xi−1 − xi+1 − . . . − x + . aii aii aii aii n aii Matici A lze zapsat ve tvaru (k+1)
xi
=−
A = D + L + R, kde D je diagon´aln´ı, L je doln´ı troj´uheln´ıkov´a a R horn´ı troj´uheln´ıkovou matici (L a R maj´ı nulovou diagon´alu). Jacobiho iteraci lze zapsat ve tvaru ~x(k+1) = −D−1(L + R)~x(k) + D−1~b. Vˇ eta Pokud je matice A diagon´alnˇe dominantn´ı, pak Jacobiho metoda konverguje. D˚ ukaz: Pro diagon´alnˇe dominantn´ı matici plat´ı n X |aij | <1 |aii |
j=1,j6=i
a tedy ˇr´adkov´a norma matice kD−1(L + R)k < 1. 17
4.7
Gauss–Seidelova metoda
Gauss–Seidelova metoda je podobn´a Jacobiho metodˇe, ale na rozd´ıl od n´ı (k+1) pouˇz´ıv´a pˇri v´ypoˇctu sloˇzek vektoru xi jiˇz dˇr´ıve vypoˇcten´e sloˇzky k + 1. iterace. Iterace je d´ana vztahem (k+1)
xi
=−
ai1 (k+1) ain (k) bi aii−1 (k+1) aii+1 (k) x1 xi−1 − xi+1 − . . . − x + . −...− aii aii aii aii n aii
Vztah lze zapsat vektorovˇe ~x(k+1) = −(D + L)−1R~x(k) + (D + L)−1~b. Vˇ eta Pro konvergenci Gauss–Seidelovy metody staˇc´ı, kdyˇz plat´ı libovoln´a z n´asleduj´ıc´ıch dvou podm´ınek: 1. A je diagon´alnˇe dominantn´ı matice 2. A je symetrick´a pozitivnˇe definitn´ı matice 4.8
Superrelaxaˇ cn´ı metoda
Gauss–Seidelova metoda konverguje pro pomˇernˇe ˇsirokou tˇr´ıdu matic, ale jej´ı (k) (k+1) (k) konvergence m˚ uˇze b´yt velmi pomal´a. Oznaˇc´ıme-li ∆xi = xi − xi rozd´ıl mezi iteracemi Gauss–Seidelovou metodou, pak je superrelaxaˇcn´ı metoda d´ana vztahem (k+1)
xi
(k)
(k)
= xi + ω∆xi ,
kde relaxaˇcn´ı faktor ω je z intervalu (0, 2), obvykle ω ∈ h1, 2). Relaxaˇcn´ı faktor slouˇz´ı k urychlen´ı konvergence metody a jeho optim´aln´ı hodnotu lze vypoˇc´ıtat ze vztahu ωopt =
1+
p
2 , 1 − ̺2 (B)
B = −(D + L)−1R.
B je iteraˇcn´ı matice Gauss–Seidelovy metody. Gauss–Seidelova metoda je tedy speci´aln´ım pˇr´ıpadem superrelaxaˇcn´ı metody s ω = 1.
18
5
Hled´ an´ı vlastn´ıch ˇ c´ısel a vektor˚ u
5.1
´ Uvod
Necht’ pro ˇc´ıslo λ existuje vektor ~x 6= ~0 takov´y, ˇze A~x = λ~x. Pak λ je vlastn´ı ˇc´ıslo a ~x je vlastn´ı vektor matice A. Dva typy ´uloh ´ y probl´em vlastn´ıch ˇc´ısel – hled´an´ı vˇsech vlastn´ıch ˇc´ısel matice a popˇr´ıpadˇe 1. Upln´ i pˇr´ısluˇsn´ych vlastn´ıch vektor˚ u ˇ asteˇcn´y probl´em vlastn´ıch ˇc´ısel – hled´an´ı 1 nebo nˇekolika vlastn´ıch ˇc´ısel 2. C´ (obvykle nejvˇetˇs´ıch) Charakteristick´ y polynom matice A: determinant det(A − λI). Pozn. Je-li A matice n × n, je charakteristick´y polynom n-t´eho stupnˇe, a tedy m´a n koˇren˚ u (mohou b´yt i v´ıcen´asobn´e). Ke kaˇzd´emu vlastn´ımu ˇc´ıslu ∃ alespoˇn 1 vlastn´ı vektor. Poˇcet l line´arnˇe nez´avisl´ych (LN) vlastn´ıch vektor˚ u je l ≤ k (k je n´asobnost vlastn´ıho ˇc´ısla). Pozn. Matice defektn´ı m´a < n line´arnˇe nez´avisl´ych vlastn´ıch vektor˚ u. Pˇr´ıklad 1 0 , det(A − λI) = (1 − λ)2 = 0, a tedy λ1,2 =1. A= 1 1 0 Vektor ~x = je jedin´ym vlastn´ım vektorem A. 1 Pozn. Re´aln´a matice m˚ uˇze m´ıt komplexnˇe sdruˇzen´a vlastn´ı ˇc´ısla a vlastn´ı vektory. Napˇr´ıklad 1 −1 A= , det(A − λI) = (1 − λ)2 + 1 = 0, a tedy λ1,2 = 1 ± i. 1 1 1 1 Vlastn´ı vektory jsou ~x1 = a ~x2 = . −i i Norm´ aln´ı matice A je takov´a, ˇze AT A = AAT . Norm´aln´ı matice ˇr´adu n m´a n LN vlastn´ıch vektor˚ u. 19
Pozn. Symetrick´a matice (A = AT ) m´a vˇsechna vlastn´ı ˇc´ısla re´aln´a. Pozn. Troj´uheln´ıkov´e matice maj´ı vˇsechna vlastn´ı ˇc´ısla na diagon´ale. Vˇ eta Podobn´e matice A a P−1AP maj´ı stejn´a vlastn´ı ˇc´ısla (stejn´e spektrum). D˚ ukaz: det(P−1AP − λI) = det[P−1(A − λI)P] = = det(P−1) det(A − λI) det(P) = det(A − λI) Pokud je vektor ~x vlastn´ım vektorem matice A, pak vektor P−1~x je vlastn´ım vektorem matice P−1AP. Vˇ eta Ke kaˇzd´e matici existuje j´ı podobn´a matice v tvaru λ J1 0 . . . 0 1 0 J 0 2 J = .. kde Ji = 0 , . . . . .. . . . . 0 0 . . . Js 0
Jordanovˇe norm´aln´ım 0 0 ... λ 0 1 λ ...
0 0 0 .. . 0 0 ... λ
.
Pozn. Pro kaˇzdou norm´aln´ı matici existuje podobn´a matice diagon´aln´ı. uv Pozn. Neexistuje finitn´ı postup jak prov´est transformaci matice na Jordan˚ norm´aln´ı tvar. Numerick´ e metody ˇ reˇ sen´ı u ´ pln´ eho probl´ emu vlastn´ıch ˇ c´ısel 1. Sekvenc´ı element´arn´ıch transformac´ı pˇrevedeme matici na pˇribliˇznˇe diagon´aln´ı tvar (pˇr´ıp. Jordan˚ uv norm´aln´ı tvar) nebo pˇribliˇznˇe speci´aln´ı typ (napˇr. tridiagon´aln´ı nebo Hessenbergovu matici). ˜ = 2. Rozloˇz´ıme matici A na souˇcin dvou matic A = FL · FR . Matice A FR FL je podobn´a matici A. −1 Odvozen´ı: FR FL = F−1 L FL FR FL = FL AFL .
20
5.2
Jacobiho transformace (metoda)
Jacobiho metoda nalezne vˇsechna vlastn´ı ˇc´ısla a vektory symetrick´ e matice. C´ılem je postupnˇe zmenˇsovat ˇc´ısla mimo diagon´alu. V kaˇzd´em kroku nalezneme v absolutn´ı hodnotˇe maxim´aln´ı ˇc´ıslo mimo diagon´alu |apq | = max |aij |. Pootoˇc´ıme osy p, q tak, aby se submatice [app apq ; aqp aqq ] i,j
transformovala na diagon´aln´ı. Jeden (k-t´y) krok iterace Jacobiho metodou m´a tvar 1 0 ... cos ϕ . . . − sin ϕ .. .. (k) T (k−1) A = Tp,q A Tp,q , kde Tp,q = . 1 . sin ϕ . . . cos ϕ ... 0 1 Po k-t´e iteraci je prvek matice
(k) (k−1) (k−1) (k−1) apq = (cos2 ϕ − sin2 ϕ) apq − cos ϕ sin ϕ (app − aqq ). (k)
Aby prvek apq = 0, ´uhel ϕ mus´ı splˇnovat vztah tg 2ϕ =
2apq . app − aqq
D˚ ukaz konvergence: Pˇri posloupnosti Jacobiho transformac´ı nejsou nuly mimo diagon´alu trval´e. D˚ ukaz konvergence je zaloˇzen na konvergenci sumy mimodiagon´aln´ıch prvk˚ u k nule. n P a2ij . Pak Oznaˇcme t(A) = i,j=1;i6=j
(k)
t(A ) = t(A
(k−1)
)−
2a2pq
∨
a2pq
t(A(k−1)) ≥ . n(n − 1)
Tedy (k)
t(A ) ≤
k 2 2 (k−1) t(A). t(A )≤ 1− 1− n(n − 1) n(n − 1) 21
.
Posloupnost je tedy shora odhadnuta geometrickou posloupnost´ı s kvocientem < 1.
5.3
LU rozklad pro u ´ pln´ y probl´ em vlastn´ıch ˇ c´ısel
Tato metoda konverguje velmi pomalu, na v´ypoˇcet je tˇreba velk´y poˇcet operac´ı. Matice A je nult´ym krokem iterace, t.j. A0 ≡ A. V k-t´em kroku rozloˇz´ıme matici Ak = Lk Uk . Vytvoˇr´ıme matici Ak+1 = Uk Lk , ta je podobn´a matici Ak . Pokud posloupnost Bk = L0 L1 . . . Lk → k regul´arn´ı matici, potom matice Ak → k horn´ı troj´uheln´ıkov´e matici, kter´a m´a na diagon´ale vlastn´ı ˇc´ısla. Existuj´ı ale speci´aln´ı rozklady matic vhodn´e pro hled´an´ı vlastn´ıch ˇc´ısel a vektor˚ u. Pro tyto rozklady konverguje obdobn´y postup rychle.
22
5.4
ˇ asteˇ C´ cn´ y probl´ em vlastn´ıch ˇ c´ısel
Hled´ame vlastn´ı ˇc´ıslo nejvˇetˇs´ı v absolutn´ı hodnotˇe. Zvol´ıme libovoln´y vektor ~x(0) . D´ale poˇc´ıt´ame iterace ~x(k+1) =
1 A~x(k), ̺k
kde
̺k = ~eT x(k) 1 A~
(pˇr´ıp. ̺k = kA~x(k)k).
Pak plat´ı lim ̺k = λ1
k→∞
∨
lim ~x(k) = ~x1.
k→∞
Pokud chceme dalˇs´ı vlastn´ı ˇc´ıslo, redukujeme matici na ˇr´ad (n−1). Je-li vektor ~x1 = (u1, . . . , un)T , plat´ı u1 0 . . . 0 u 1 0 λ1 ~qT 2 −1 P = .. , P AP = . . . ... ~0 B , . un 0 . . . 1 a tedy hled´ame maxim´aln´ı vlastn´ı ˇc´ıslo matice B. Pozn. Nev´yhodou je postupn´a ztr´ata pˇresnosti. Pozn. Pro v´ypoˇcet nejmenˇs´ıho vlastn´ıho ˇc´ısla lze nal´ezt nejvˇetˇs´ı vlastn´ı ˇc´ıslo matice A−1. Pro hled´an´ı vlastn´ıho ˇc´ısla v urˇcit´e oblasti lze prov´est posun, nebot’ (A + µI)~x = (λ + µ)~x.
23