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 1 • 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.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
~x =
x1 x2 .. . xn
= (x1, x2, . . . , xn)T
A=
a11 a12 . . . a1m a21 a22 . . . a2m .. .. .. .. . . . . an1 an2 . . . anm
= (aij )
Pokud nebude uvedeno jinak, budeme uvaˇzovat re´aln´e matice a vektory. 1
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 a c) k~x + ~yk ≤ 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
x2i
i=1
Maticov´a norma Vektorov´a norma matice se naz´yv´a maticovou normou, pokud nav´ıc plat´ı pro ∀ matice A, B podm´ınka 4) pro souˇcin matic kA · Bk ≤ kAk · kBk Maticov´a norma je souhlasn´a s vektorovou normou, pokud pro ∀ A, ~x plat´ı kA~xk ≤ kAk · k~xk Pˇr´ıklady maticov´ych norem • ˇr´adkov´a - kAkI = max
n X
i=1,...,n j=1
• sloupcov´a - kAkII = max
|aij | n X
j=1,...,n i=1
• Eukleidovsk´a - kAkIII =
|aij |
v u n n uX X u a2ij t i=1 j=1
Kaˇzd´a z uveden´ych maticov´ych norem je souhlasn´a se stejnˇe oznaˇcenou vektorovou normou.
2
1.3
Metody ˇ reˇ sen´ı soustav line´ arn´ıch rovnic
1. pˇr´ım´e (finitn´ı) 2. iteraˇcn´ı 3. gradientn´ı
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 sˇc´ıtan´ı), poˇcet operac´ı roste tedy ∼ n2 (pˇresnˇeji ≃ 0.5 · n2 vnitˇrn´ıch cykl˚ u).
3
2.2
Gaussova a Gauss-Jordanova eliminace
ˇ s´ım soustavu rovnic A~x = ~b. V prvn´ım kroku necht’ prvek a11 6= 0 (lze vˇzdy Reˇ dos´ahnout pˇrehozen´ım rovnic). Prvek a11 , pouˇzit´y k ´upravˇ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´ı
D1 =
1 −a21/a11 .. . −an1 /a11
0 1 .. . 0
0 ... 0 ... .. .. . . 0 ...
0 0 .. . 1
Rovnice po prvn´ı ´upravˇe m´a tvar D1 A~x = D1~b, oznaˇc´ıme D1 A ≡ A(1) a D1~b ≡ ~b(1) . Po k − 1 ´uprav´ach m´a matice A(k−1) tvar
A(k−1) =
a11 a12 (1) 0 a22 .. .. . . 0 0 0 0 0 0 .. .. . . 0 0
. . . a1,k−1 (1) . . . a2,k−1 .. .. . . (k−2) . . . ak−1,k−1 ... 0 ... 0 .. .. . . ... 0
a1k (1) a2k .. . (k−2) ak−1,k (k−1) akk (k−1) ak+1,k .. . (k−1) 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) (k−1) 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 - pˇr´ım´e metody nepouˇziteln´e pro obecn´e matice!!
4
Poˇcet operac´ı Na kaˇzdou 0 →≤ n vnitˇrn´ıch cykl˚ u, potˇrebuji n(n − 1) prvk˚ u → 0. Celkov´y 3 3 poˇcet vnitˇrn´ıch cykl˚ u ∼ n (pˇresnˇeji ≃ 1/3 n ), sloˇzitost algoritmu je ˇr´adu n3 . 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.
2.3
V´ ybˇ er hlavn´ıho prvku (pivoting)
V ∀ kroku pˇr´ım´eho bˇehu → v´ybˇer hlavn´ıho prvku. ´ y v´ybˇer hlavn´ıho prvku - v cel´e dosud neupraven´e ˇc´asti matice, po• Upln´ mal´y, strategie - max |aij | ˇ asteˇcn´y v´ybˇer hlavn´ıho prvku - v dan´em sloupci (sloupcov´y) nebo ˇr´adku • C´ (ˇr´adkov´y), rychlejˇs´ı, vylepˇsen´a strategie napˇr. 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 u ´prav´ach 3. Ztr´ata pˇresnosti
5
2.4
LU metoda
∀ matice 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 a21 a31 ... an1
a12 a22 a32 ... an2
a13 a23 a33 ... an3
... ... ... ... ...
a1n a2n a3n ... ann
=
1 0 0 l21 1 0 l31 l32 1 ... ... ... ln1 ln2 ln3
... 0 ... 0 ... 0 ... ... ... 1
u11 u12 0 u22 0 0 ... ... 0 0
u13 u23 u33 ... 0
... ... ... ... ...
N´asoben´ı matic i ≤ j
aij = uij +
i−1 X
lik ukj
k=1 j−1 X
i > j
aij = lij ujj +
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
lij = aij −
j−1 X
k=1
lik ukj /ujj
i = j + 1, . . . , n
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 1 ×, v´ysledn´e prvky matic L a U se vejdou do 1 matice. 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 6
u1n u2n u3n ... unn
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 ~xi+1 = ~xi + (δx) 2.6
~ i = ~b − A~xi A(δx)
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 ´ulohu (A + ∆A) (~x + ∆~x) = ~b + ∆~b Nejdˇr´ıve pˇr´ıpad ∆A = 0 ∆~x = A−1∆~b ⇒ k∆~xk ≤ kA−1k · k∆~bk A~x = ~b ⇒ k~xk ≥ k~bk/kAk Pro relativn´ı chybu ˇreˇsen´ı tedy plat´ı k∆~xk k∆~bk −1 ≤ kAk · kA k · k~xk k~bk ˇ ıslo Cp = kAk · kA−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´ı.
7
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 metody 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 ´upravy 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. 2.8
Speci´ aln´ı typy matic
ˇ ıdk´a matice m´a vˇetˇsinu prvk˚ u = 0. R´ Pro ˇreˇsen´ı soustav s ˇr´ıdkou matic´ı se ˇcasto pouˇz´ıvaj´ı gradientn´ı metody, spoˇc´ıvaj´ıc´ı v minimalizaci 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. Tridiadon´aln´ı matice pro p = 1, pˇetidiagon´aln´ı matice pro p = 2. Soustavy s tridiagon´aln´ı matic´ı
a1 c2 0 .. . 0 0
b1 a2 c3 .. . 0 0
0 b2 a3 .. . 0 0
0 0 b3 .. . 0 0
... 0 0 0 ... 0 0 0 ... 0 0 0 .. .. .. .. . . . . . . . cn−1 an−1 bn−1 ... 0 cn an
·
x1 x2 x3 .. . xn−1 xn
=
f1 f2 f3 .. . fn−1 fn
Tridiagon´aln´ı 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). 8
ˇ 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. Blokovˇe tridiagon´aln´ı matice - Ai, Bi, Ci - mal´e matice ⇒ µi , ρi - 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´ı, 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, W diagon´aln´ı, I jednotkov´a, U, V ortogon´aln´ı (UT · U = VT · V = I) A = U · W · VT
⇒
~x = V · [diag(1/wj )] · UT · ~b
Pokud wj = 0 (wj ≃ 0) nahrad´ıme 1/wj → 0 (signalizuje singul´arnost matice).
9
2.10
Rychlost ˇ reˇ sen´ı soustav line´ arn´ıch rovnic (matice N x N ) pˇ r´ım´ ymi metodami
• Gauss–Jordanova eliminace Pˇri v´ypoˇctu se prov´ad´ı ∼ N 3 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 N 3 vnitˇrn´ıch cykl˚ u. Pˇri zpˇetn´em 1 bˇehu ∼ 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´ı ∼ N 3 , 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.
3
Gradientn´ı metody
Line´arn´ı rovnici A~x = ~b ˇreˇs´ıme napˇr´ıklad minimalizac´ı funkce f (~x) = 21 |A~x − ~b|2 . V kaˇzd´em kroku λ takov´e, aby f (~x + λ~u) bylo minim´aln´ı. Tedy λ=
−~u∇f |A~u|2
kde
∇f (~x) = AT (A~x − ~b)
Pro ˇr´ıdk´e matice se sloˇzitost n´asoben´ı vektoru matic´ı sniˇzuje z poˇctu operac´ı ∼ N 2 na poˇcet operac´ı ∼ N . Pozn. Existuje ˇrada modern´ıch ˇcasto pouˇz´ıvan´ych gradientn´ıch metod.
10
4 4.1
Iteraˇ cn´ı metody ˇ reˇsen´ı soustav line´ arn´ıch rovnic 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
Symetrick´a matice A typu n × n je pozitivnˇ e definitn´ı, pokud pro ∀~x 6= ~0 je skal´arn´ı souˇcin (~x, A~x) > 0. 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) + ~c(k) 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)
11
Matice Bk a vektor ~ck maj´ı tvar
Bk =
1
−ai1
... 1 . . . −aii−1 0 −aii+1 . . . 1 ...
−ain ,
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
12
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
k~bk (k) k (0) k~x − ~xk ≤ kBk k~x 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 (k+1)
xi
=−
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 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 j=1,j6=i |aii | a tedy maximov´a norma matice kD−1(L + R)k < 1. 13
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. Matice A je pozitivnˇe definitn´ı (symetrick´a s kladn´ymi vlastn´ımi ˇc´ısly). 4.8
Superrelaxaˇ cn´ı metoda
Gauss–Seidelovy 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 =
2 1 + 1 − ̺2 (B) q
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.
14
5
Hled´ an´ı vlastn´ıch ˇ c´ısel a vektor˚ u
5.1
´ Uvod
Necht’ pro ˇc´ıslo λ ∃~x 6= ~0 takov´y, ˇze A~x = λ~x. Pak λ je vlastn´ı ˇc´ıslo a vektor ~x je vlastn´ı vektor matice A. 2 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
A=
1 0 1 1
Vektor ~x =
0 1
,
det(A − λI) = (1 − λ)2 = 0,
tedy λ1,2 = 1.
je jedin´ym vlastn´ım vektorem A.
uˇze m´ıt komplexnˇe sdruˇzen´a vlastn´ı ˇc´ısla a vlastn´ı vekPozn. Re´aln´a matice m˚ tory. Napˇr´ıklad
A=
1 −1 1 1
,
det(A − λI) = (1 − λ)2 + 1 = 0,
Vlastn´ı vektory jsou ~x1 =
1 −i
a ~x2 =
tedy λ1,2 = 1 ± i.
1 . 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.
15
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). Odvozen´ı 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´a matici ∃ j´ı podobn´a matice v Jordanovˇe norm´aln´ım tvaru
J=
J1 0 . . . 0 0 J2 0 .. . . . ... . 0 0 . . . Js
,
kde
Ji =
λ 1 0 .. . 0
0 0 ... λ 0 1 λ ...
0 0 0 .. . 0 0 ... λ
.
Pozn. Pro ∀ norm´aln´ı matici ∃ 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. Sekvence element´arn´ıch transformac´ı −→ 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 Hessenbergova matice). ˜ = 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 .
16
5.2
Jacobiho transformace (metoda)
Jacobiho metoda nalezne ∀ 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 matice [app apq ; aqp aqq ] i,j
transformovala na diagon´aln´ı. Jeden (k-t´y) krok iterace Jacobiho metodou m´a tvar 1 A
(k)
=
(k−1) TT Tp,q , p,q A
kde
Tp,q =
0
... c . . . −s .. .. . 1 . s ... c
...
0 1 V t´eto matici plat´ı c = cos ϕ a s = sin ϕ. Po k-t´e iterace je prvek matice
.
(k) (k−1) (k−1) (k−1) apq = (c2 − s2 )apq − cs(app − aqq ) (k) Aby se prvek apq = 0, ´uhel ϕ mus´ı splˇnovat vztah 2apq tg2ϕ = app − aqq
D˚ ukaz konvergence Pˇri posloupnosti Jacobiho transformac´ı nejsou 0 mimo diagon´alu trval´e. D˚ ukaz konvergence je zaloˇzen na konvergenci sumy mimodiagon´aln´ıch prvk˚ u k 0. 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
2 2 t(A(k−1) ) ≤ 1 − t(A) t(A ) ≤ 1 − n(n − 1) n(n − 1) Posloupnost je tedy shora odhadnuta geometrickou posloupnosti s kvocientem < 1. (k)
17
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 0-t´ym krokem iterace, tj. 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 = L0L1 . . . Lk → k regul´arn´ı matici ⇒ 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. 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) . A d´ale poˇc´ıt´ame iterace 1 ~x(k+1) = A~x(k) , kde ̺k = ~eT x(k) ( pˇr´ıp. 1 A~ ̺k Pak plat´ı lim ̺k = λ1
k→∞
∨
̺k = kA~x(k) 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´ı
P=
u1 0 . . . 0 u2 1 0 .. . . . ... . un 0 . . . 1
,
P−1AP =
λ1 ~qT ~0 B
.
d´ale 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. 18
