Tomáš
Oberhuber
Iterativní
metody
Tomáš Oberhuber
Faculty of Nuclear Sciences and Physical Engineering
Czech Technical University in Prague
1 / 62
Tomáš
Oberhuber
Iterativní
metody
Výhody iterativních metod
Jaké jsou hlavní duvody
˚
studia iterativních metod pro úlohu
~x = ~b?
A
• GEM má složitost n3 , pro velké matice je témeˇ
ˇr
nepoužitelná
∞
• iterativní metody generují posloupnost ~
x (k ) k =1 , která
konverguje k ˇrešení soustavy
• napoˇcítání nového ~
x (k ) je založeno na násobení matice
a vektoru nebo velice podobné operaci, která má
složitost n2
• pokud dostaneme dostateˇcneˇ dobrou aproximaci po
méneˇ než n krocích, muže
˚ být iterativní metoda
rychlejší
• výhoda iterativních metod je také v tom, že nemení
ˇ
samotnou maticí
implementaci
A, díky tomu jsou snažší na
2 / 62
Tomáš
Oberhuber
Výhody iterativních metod
Iterativní
metody
• pˇri šikovném uložení ˇrídkých matic muže
˚ mít násobení
matice a vektoru ješteˇ menší nároˇcnost, než n2 , ale pˇri
aplikaci GEM na ˇrídkou matici muže
˚ vzniknout mnoho
nových nenulových prvku,
˚ což hraje ješteˇ více v
ˇ GEM
neprospech
• v nekterých
ˇ
situacích ani nemáme matici
A
explicitneˇ
ˇ pak je GEM nemožné použít
uloženou v pameti,
3 / 62
Tomáš
Oberhuber
Iterativní
metody
Výhody iterativních metod
• pokud ˇrešíme nelineární úlohu F ~
x , ~d
= 0, cˇ asto
provádíme linearizaci, jejímž výsledkem je, že ˇrešíme
posloupnost soustav
A(l)~x = ~b(l)
• oznaˇcíme-li ˇrešení této soustavy~
x (l) ,pak posloupnost
(l) ∞
~x
konverguje k ˇrešení F ~x , ~d = 0
l=1
• jednotlivé matice
• pokud je
A(l) a pravé strany ~b(l) se liší jen málo
ˇ
A(l) rozumneˇ podmínená,
pak ˇrešení ~x (l−1) je
blízko ˇrešení ~x (l) a tudíž je ~x (l−1) velmi dobrým
poˇcáteˇcním odhadem pro iterativní metody ˇrešící
soustavu (l)~x = ~b(l)
• pokud u GEM zmeníme
ˇ
jen trochu matici soustavy
musíme provést celou eliminaci znovu
A
A,
4 / 62
Tomáš
Oberhuber
Nevýhody iterativních metod
Iterativní
metody
• vetšina
ˇ
z nich nedá pˇresné ˇrešení po koneˇcném poˇctu
kroku˚
• víme ale, že GEM prakticky také ne
• musíme umet
ˇ rozhodnout, kdy zastavit napoˇcítávání
dalších iterací ~x (k )
• poˇcet nutných iterací silneˇ závisí na konkrétní matici a
tím pádem i doba výpoˇctu, u GEM doba výpoˇctu na
ˇ r nezávisí
samotné matici témeˇ
• proto se GEM lépe hodí do nekterých
ˇ
kritických
aplikací, kde je nutné zajistit dostateˇcneˇ rychlou odezvu
• pro malé matice je GEM cˇ asto rychlejší
5 / 62
Tomáš
Oberhuber
Iterativní metody obecneˇ
Iterativní
metody
• obecneˇ platí, že volíme libovolné ~
x (0) a dále
napoˇcítáváme
~x (k +1) =
kde volba
B(k )~x (k ) + ~c (k ),
B(k ) a ~c (k ) záleží na použité metodeˇ
• dále požadujeme, aby pro všechna k platilo
B(k )~x ∗ + ~c (k ),
kde ~x ∗ je pˇresné ˇrešení soustavy A~x = ~b
~x ∗ =
6 / 62
Tomáš
Oberhuber
Iterativní
metody
Iterativní metody obecneˇ
Theorem 1
Iterativní metoda tvaru
~x (k +1) =
ˇ
splnující
~x ∗ =
B(k )~x (k ) + ~c (k ),
B(k )~x ∗ + ~c (k ),
konverguje k ~x ∗ pro libovolné ~x (0) práveˇ když,
lim
k →∞
B(k )B(k −1) . . . B(0) = θ.
Dukaz.
˚
Platí
~x (k ) − ~x ∗ =
=
=
=
B(k −1)~x (k −1) + ~c (k −1) − B(k −1)~x ∗ − ~c (k −1)
~∗
B(k −1)~x(k −1) − B(k −1)
x
(k −1) ~ (k −1)
∗
B
x
− ~x
B(k −1)B(k −2) . . . B(0) ~x (0) − ~x ∗
7 / 62
Tomáš
Oberhuber
Iterativní
metody
Citlivost a stabilita iterativních
metod
• pokud behem
ˇ
ˇ
ˇ tˇreba
výpoˇctu dojde k nejaké
chybe,
vinou zaokrouhlení, mužeme
˚
se na výsledný vektor
˚ být libovolné a
dívat jako na nové ~x (0) , které muže
metoda bude konvergovat i tak
• chyby ve výpoˇctech se tedy nekumulují, ale eliminují
• jde o tzv. samoopravující vlastnost iterativních metod
• tím pádem již dále nemusíme studovat citlivot
iterativních metod na výpoˇcetní chyby
8 / 62
Tomáš
Oberhuber
Iterativní metody obecneˇ
Iterativní
metody
ˇ na:
Iterativní metody se delí
B(k ) a ~c (k ) jsou konstantní tj.
~x (k +1) = B~x (k ) + ~c
• nestacionární – B(k ) a ~c (k ) jsou ruzné
˚
pro každé k
• stacionární –
Výhoda stacionárních metod je, že pracujeme jen s jednou
maticí. Jsou proto jednodušší na implementaci a také na
numerickou analýzu. Dále budeme studovat pouze
stacionární metody.
9 / 62
Tomáš
Oberhuber
Iterativní
metody
Stacionární iterativní metody
Theorem 2
Metoda daná vztahem
~x (k +1) =
ˇ
a splnující
~x ∗ =
B~x (k ) + ~c
B~x ∗ + ~c ,
konverguje k ~x ∗ pro libovolné ~x (0) práveˇ když platí,
lim
k →∞
Bk = θ.
Dukaz.
˚
ˇ 1a
Plyne z vety
lim
k →∞
B(k −1) . . . B(0) = klim
Bk .
→∞
10 / 62
Tomáš
Oberhuber
Stacionární iterativní metody
Iterativní
metody
Theorem 3
Metoda daná vztahem
~x (k +1) =
ˇ
a splnující
~x ∗ =
B~x (k ) + ~c ,
B~x ∗ + ~c ,
konverguje k ~x ∗ pro libovolné ~x (0) práveˇ když platí, že
ρ ( ) < 1.
B
Dukaz.
˚
Plyne z podmínek pro konvergenci posloupnosti
Bk .
11 / 62
Tomáš
Oberhuber
Stacionární iterativní metody
Iterativní
metody
Theorem 4
Postaˇcující podmínkou pro to, aby metoda daná vztahem
~x (k +1) =
ˇ
a splnující
~x ∗ =
B~x (k ) + ~c ,
B~x ∗ + ~c ,
B
ˇ
konvergovala k ~x ∗ pro libovolné ~x (0) je, že k k < 1 v nejaké
maticové normeˇ k·k.
Dukaz.
˚
Plyne z podmínek pro konvergenci posloupnosti
Bk .
12 / 62
Tomáš
Oberhuber
Iterativní
metody
Stacionární iterativní metody
Theorem 5
Aposteriorní odhady chyb pro stacionární iterativní
metody: Pro metodu danou vztahem
~x (k +1) =
ˇ
a splnující
~x ∗ =
B~x (k ) + ~c ,
B~x ∗ + ~c ,
A
kde ~x ∗ je rˇešením soustavy lineárních rovnic ~x = ~b, platí
následující odhady chyby aproximace rˇešení:
(k )
• ~
x − ~x ∗ ≤ −1 ~x (k ) − ~b = −1 ~r (k ) ,
(k )
• ~
x − ~x ∗ ≤ ( − )−1 k k ~x (k −1) − ~x (k ) ,
A A
I B
B
A
pˇri použití souhlasné maticové normy s danou vektorovou
normou.
Remark 6
Zde jsme použili definici rezidua v k -té iteraci
~r (k ) =
A~x (k ) − ~b.
13 / 62
Tomáš
Oberhuber
Stacionární iterativní metody
Iterativní
metody
Remark 7
• aposteriorní odhady chyb jsou duležité
ˇ kdy
˚
pro zjišt
ení,
výpoˇcet metody zastavit a to sice tehdy, je-li ~r (k ) nebo ~x (k −1) − ~x (k ) dostateˇcneˇ malé
• pojem "dostateˇcneˇ malé" se cˇ asto vztahuje k normeˇ
matice
napˇr.
A nebo pravé strany ~b rovnice, tj. požadujeme
(k ) ~r <
~ b
14 / 62
Tomáš
Oberhuber
Iterativní
metody
Stacionární iterativní metody
Theorem 8
Apriorní odhad chyby pro stacionární iterativní metody:
ˇ
ˇ Pak pro
Bud’ kBk < 1v nejaké
maticové norme.
∞
posloupnost ~x (k ) k =1 generovanou pˇredpisem
~x (k +1) =
ˇ
a splnující
~x ∗ =
B~x (k ) + ~c ,
B~x ∗ + ~c ,
platí
"
~ (k ) ~ ∗ k ~ (0) x − x ≤ k k x +
B
#
~c ,
1−k k
B
kde použitá maticová norma je souhlasná s použitou
vektorovou normou.
15 / 62
Tomáš
Oberhuber
Stacionární iterativní metody
Iterativní
metody
• nyní již víme, za jakých podmínek stacionární metoda
ˇ
konverguje k nejakému
vektoru ~x ∗ a známé odhady
chyb
• dále se budeme zabývat otázkou, jak volit matici
vektor ~c , aby vztah
~x ∗ =
Ba
B~x ∗ + ~c ,
platil pro ˇrešení ~x ∗ soustavy lineárních rovnic
A~x = ~b
16 / 62
Tomáš
Oberhuber
Iterativní
metody
Metoda postupných aproximací
• pˇri návrhu nejjednodušší iterativní metody chceme najít
ˇ
nejaký
zpusob,
˚
jak zlepšit urˇcitou aproximaci ˇrešení ~x (k )
• potˇrebujeme tedy odhadnout, jaké chyby se s danou
aproximací dopouštíme
• jediný odhad, který máme k dispozici, je reziduum
~r (k ) = ~x (k ) − ~b
A
• zkusme tedy metodu
A
I A) ~x (k ) + ~b
~x (k +1) = ~x (k ) −~r (k ) = ~x (k ) − ~x (k ) + ~b = ( −
• vidíme, že metoda splnuje
ˇ
podmínku
I A) ~x ∗ + ~b
~x ∗ = ( −
• je tedy
B
=
⇔
A~x ∗ = ~b
I − A,
~c = ~b
17 / 62
Tomáš
Oberhuber
Iterativní
metody
Metoda postupných aproximací
• tuto metodu nazveme metodou postupných
aproximací
• snadno vidíme, že platí
Theorem 9
Metoda postupných aproximací pro rovnici
regulární maticí ) ve tvaru
A
A~x = ~b (s
I A) ~x (k ) + ~b,
~x (k +1) = ( −
konverguje pro každé ~x (0) k rˇešení zadané rovnice práveˇ
když
ρ ( − ) < 1.
I A
Je-li
I Ak < 1,
k −
ˇ
pro nejakou
maticovou normu k·k, pak tato metoda také
konverguje k rˇešení soustavy pro libovolné ~x (0) .
18 / 62
Tomáš
Oberhuber
Metoda postupných aproximací
Iterativní
metody
Theorem 10
ˇ
Bud’ p (x) polynom promenné
x. Bud’
Pak p (λ) ∈ σ (p ( )).
A
A ∈ C n,n , λ ∈ σ (A).
Dukaz.
˚
Lze ukázat snadno s využitím Jordanova tvaru matice.
19 / 62
Tomáš
Oberhuber
Iterativní
metody
Metoda postupných aproximací
Theorem 11
A
Bud’ hermitovská a PD. Pak metoda postupných
aproximací konverguje práveˇ když platí
I A > 0.
2 >
Dukaz.
˚
A
I A
I A
Je-li hermitovská a má-li být ρ ( − ) < 1, pak musí být
ˇ
σ ( − ) ⊂ (−1, 1) a podle pˇredchozí vety
A
I
σ ( ) ∈ (0, 2) ⇔ 2 > A > 0.
Remark 12
Takových matic ale moc není.
20 / 62
Tomáš
Oberhuber
ˇ
Pˇredpodmínená
metoda
postupných aproximací
Iterativní
metody
• abychom zlepšili konvergenci metody postupných
ˇ
aproximací, použijeme tzv. pˇredpodmínení
~x = ~b vynásobíme vhodnou regulární maticí
• soustavu
H
A
• dostaneme soustavu
HA~x = H~b
• ˇrešení se tím nezmení
ˇ
21 / 62
Tomáš
Oberhuber
ˇ
Pˇredpodmínená
metoda
postupných aproximací
Iterativní
metody
• metoda postupných aproximací má nyní mít tvar
~x (k +1) = ~x (k ) − ~r (k )
= ~x (k ) − ~x (k ) + ~b
~x (k ) − ~b
= ~x (k ) −
=
HA
H
H A
(I − HA) ~x (k ) + H~b
• a samozˇrejmeˇ platí podmínka
I HA) ~x ∗ + H~b = ~x ∗ − HA~x ∗ + H~b = ~x ∗
( −
22 / 62
Tomáš
Oberhuber
Iterativní
metody
ˇ
Pˇredpodmínená
metoda
postupných aproximací
I HA) ~x (k ) + H~b vidíme, že
B = I − HA,
~c = H~b
• ze vztahu ~
x (k +1) = ( −
• platí tedy následující vety
ˇ
23 / 62
Tomáš
Oberhuber
ˇ
Pˇredpodmínená
metoda
postupných aproximací
Iterativní
metody
Theorem 13
ˇ metoda postupných aproximací pro rovnici
Pˇredpodmínená
~x = ~b (s regulární maticí ) s pˇredpodmínením
ˇ
ve tvaru
A
A
H
I HA) ~x (k ) + H~b,
~x (k +1) = ( −
konverguje pro každé ~x (0) k rˇešení zadané rovnice práveˇ
když
ρ( −
) < 1.
I HA
Je-li
I HAk < 1,
k −
ˇ
pro nejakou
maticovou normu k·k, pak tato metoda také
konverguje k rˇešení soustavy pro libovolné ~x (0) .
24 / 62
Tomáš
Oberhuber
ˇ
Pˇredpodmínená
metoda
postupných aproximací
Iterativní
metody
Theorem 14
A
ˇ metoda
Bud’ hermitovská a PD. Pak pˇredpodmínená
postupných aproximací konverguje, pokud platí
W + W∗ > A > 0,
W H
kde
= −1 . Konvergence je navíc monotónní vzhledem
k vektorové normeˇ k·kA .
Remark 15
ˇ je
Pro metodu postupných aproximací bez pˇredpodmínení,
= a dostáváme již dokázané kritérium
W I
I A > 0.
2 >
ˇ nám ale neˇríká nutnou podmínku!!!
Tato veta
25 / 62
Tomáš
Oberhuber
ˇ
Pˇredpodmínená
metoda
postupných aproximací
Iterativní
metody
H
+ H~b dostáváme, že
• nyní vyvstává otázka, jak volit matici
• ze vztahu
pˇri volbeˇ
~x (k +1)
H=A
=
−1
~x (k )
je
−
HA
~x (k +1) = ~x (k ) −
= ~x (k ) −
~x (k )
HA~x (k ) + H~b
A−1A~x (k ) + A−1~b
= ~x (k ) − ~x (k ) + ~x ∗
= ~x ∗
• metoda by tak konvergovala po jedné iteraci
• získat
A−1 ale umíme jen pomocí GEM, cˇ emuž jsme
ˇ vyhnout
se chteli
• umení
ˇ je najít matici tak, aby co nejlépe
aproximovala −1 ale byla snadno dosažitelná
A
H
26 / 62
Tomáš
Oberhuber
Iterativní
metody
ˇ
Pˇredpodmínená
metoda
postupných aproximací
• víme, že napoˇcítání LU rozkladu je ekvivalentní znalosti
inverzní matice
• oblíbenou a úˇcinnou metodou pˇredpodmínení
ˇ je
neuplný LU rozklad, kdy prvky LU rozkladu, které jsou
v absolutní hodnoteˇ malé zanedbáváme
• neúplný LU rozklad funguje jako pojítko mezi pˇrímými a
iterativními metodami
• pokud napoˇcítám dost pˇresný LU rozklad, získám
rychlou konevergenci, ale strávím hodneˇ cˇ asu
výpoˇctem samotného rozkladu
• toto pˇredpodmínení
ˇ se ale pˇríliš nepoužívá v kombinaci
se stacionárními metodami
• následneˇ si ukážeme nekolik
ˇ
stacionárních metod a
ukážeme si, jak pro neˇ vypadá matice
H
• tyto matice se pak cˇ asto používají jako pˇredpodmínení
ˇ
ˇ
u nekterých
jiných iterativních metod
27 / 62
Tomáš
Oberhuber
Richardsonovy iterace
Iterativní
metody
• tato metoda je dána pˇredpisem
A
~x (k +1) = ~x (k ) − θ( ~x (k ) − ~b),
pro θ ∈
C
• parametru θ se ˇríká relaxacní
ˇ parametr
• snadno vidíme, že
H
B
I
= θ
=
−θ
~c = θ~b
I
A
28 / 62
Tomáš
Oberhuber
Iterativní
metody
Richardsonovy iterace
Theorem 16
A
Je-li hermitovská a PD, pak metoda Richardsonových
iterací konverguje práveˇ když
I A > 0.
2
>
θ
Konvergence je navíc monotónní vzhledem k vektorové
normeˇ k·kA .
Dukaz.
˚
W
ˇ 14, kde
Dúkaz postaˇcující podmínky plyne z vety
=
Že jde o nutnou podmínku lze ukázata podobneˇ jako u
metody postupných aproximací.
Remark 17
Pro hermitovskou a PD matici
θ<
1
θ
I.
A metoda konverguje pro
2
.
ρ( )
A
29 / 62
Tomáš
Oberhuber
Jacobiho metoda
Iterativní
metody
• navrhl ji Carl Gustav Jakob Jacobi (1845)
• u Jacobiho metody napoˇcítáváme i-tou složku vektoru
~x (k +1) z vektoru ~x (k ) tak, aby byla splnena
ˇ
i-tá rovnice
soustavy tj.
(k +1)
xi
=
1 (k )
(k )
bi − ai1 x1 − ai2 x2 − . . .
aii
(k )
(k )
(k )
−ai,i−1 xi−1 − ai,i+1 xi+1 − . . . − ain xn
,
pro i = 1, . . . , n.
30 / 62
Tomáš
Oberhuber
Iterativní
metody
Jacobiho metoda
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
31
32
33
34
35
36
37
bool jacobiMethod ( const M a t r i x& A ,
const V e c t o r& b ,
V e c t o r& x ,
double eps )
{
double normB = norm ( b ) ;
double r = eps + 1 . 0 ;
Vector y ;
while ( r > eps )
{
/ ∗ ∗∗∗
∗ Jacobiho i t e r a c e
∗/
f o r ( i n t i = 0 ; i < n ; i ++ )
{
double s = 0 . 0 ;
f o r ( i n t j = 0 ; j < n ; j ++ )
i f ( j != i )
s = s + A[ i ] [ j ] ∗ x [ j ] ;
y [ i ] = ( b [ i ] − s ) / A[ i ] [ i ] ;
}
/ ∗ ∗∗∗
∗ Vypocet r e z i d u a
∗/
r = 0.0;
f o r ( i n t i = 0 ; i < n ; i ++ )
{
double Ax_i = 0 . 0 ;
f o r ( i n t j = 0 ; j < n ; j ++ )
Ax_i = Ax_i + A [ i ] [ j ] ∗ y [ j ] ;
r = r + ( Ax_i − b _ i ) ∗ ( Ax_i − b _ i ) ;
x[ i ] = y[ i ];
}
r = s q r t ( r ) / normB ;
}
}
31 / 62
Tomáš
Oberhuber
Jacobiho metoda - numerická
analýza
Iterativní
metody
• matici
kde
•
•
•
A pˇrepíšeme ve tvaru
A = D − L − R,
D je diagonála matice A
L jsou záporneˇ vzaté prvky matice A pod diagonálou
R jsou záporneˇ vzaté prvky matice A nad diagonálou
32 / 62
Tomáš
Oberhuber
Iterativní
metody
Jacobiho metoda - numerická
analýza
• pro numerickou analýzu lze tuto metodu maticoveˇ
zapsat takto
~x (k +1) =
=
=
D−1 ~b + (D − A) ~x (k )
D−1 (L + R) ~x (k ) + D−1~b
I − D−1A ~x (k ) + D−1~b
h
i
33 / 62
Tomáš
Oberhuber
Iterativní
metody
Jacobiho metoda - numerická
analýza
• vidíme tedy, že
B
=
~c =
H
=
I − D−1A
D−1~b
D−1
34 / 62
Tomáš
Oberhuber
Iterativní
metody
Jacobiho metoda - numerická
analýza
Theorem 18
Nutnou a postaˇcující podmínkou pro to, aby Jacobiho
metoda konvergovala pro libovolné ~x (0) je
−1
ρ
( + ) < 1.
D L R
Postaˇcující podmínkou je
−1
( +
D L R) < 1,
v libovolné maticové normeˇ k·k.
35 / 62
Tomáš
Oberhuber
Iterativní
metody
Jacobiho metoda - numerická
analýza
Definition 19
Maticí s pˇrevládající diagonálou nazveme takovou matici,
pro kterou platí
X
aij < |aii | ,
j=1,...,n;j6=i
pro všechna i = 1, . . . n.
36 / 62
Tomáš
Oberhuber
Jacobiho metoda - numerická
analýza
Iterativní
metody
Theorem 20
A
Má-li matice pˇrevládající diagonálu, pak Jacobiho metoda
konverguje k rˇešení soustavy rovnic ~x = ~b pro libovolnou
volbu ~x (0) .
A
Theorem 21
A
Je-li matice hermitovská a PD, pak Jacobiho metoda
konverguje k rˇešení soustavy rovnic ~x = ~b pro libovolnou
volbu ~x (0) , práveˇ když platí 2 > > 0 tj. a 2 − jsou
pozitivneˇ definitní. Konvergence je navíc monotónní
vzhledem k vektorové normeˇ k·kA .
D A
A
A D A
Dukaz.
˚
ˇ 14, kde
Dukaz
˚
postaˇcující podmínky plyne z vety
W = H−1 = D. Že jde o podmínku nutnou lze dokázat
podobneˇ jako u metody postupných aproximací, není to ale
triviální.
37 / 62
Tomáš
Oberhuber
Iterativní
metody
Jacobiho metoda - pˇríklad
Example 22
A




 
−2 1
0
−1
1
=  1 −2 1  , ~b =  0  , ~x =  1  ,
0
1 −2
−1
1
k
0
1
2
3
4
5
6
...
~x k
(0, 0, 0)
( 12 , 0, 21 )
( 21 , 21 , 21 )
( 43 , 21 , 43 )
( 43 , 43 , 43 )
( 87 , 43 , 87 )
( 87 , 87 , 87 )
...
~r k
(1, 0, 1)
(0, 1, 0)
( 12 , 0, 21 )
(0, 12 , 0)
( 14 , 0, 41 )
(0, 14 , 0)
( 18 , 0, 81 )
...
38 / 62
Tomáš
Oberhuber
Iterativní
metody
Gaussova-Seidelova metoda
• roku 1823 jí popsal C. F. Gauss v jednom dopise svému
studentovi Ch.L.Gerlingovi
• roku 1874 jí publikoval P. L. von Seidel
• narozdíl od Jacobiho metody, tato metoda využívá k
(k +1)
již napoˇcítané
napoˇcítání nové složky vektoru ~xi
(k +1)
(k +1)
složky tohoto vektoru tj. ~x1
, . . . , ~xi−1
• metoda je tedy dána pˇredpisem
(k +1)
xi
=
1 (k +1)
(k +1)
bi − ai1 x1
− ai2 x2
− ...
aii
(k +1)
−ai,i−1 xi−1
(k )
(k )
− ai,i+1 xi+1 − . . . − ain xn
,
pro i = 1, . . . , n.
39 / 62
Tomáš
Oberhuber
Iterativní
metody
Gaussova-Seidelova metoda
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
31
32
33
34
35
bool gaussSeidelMethod ( const M a t r i x& A ,
const V e c t o r& b ,
V e c t o r& x ,
double eps )
{
double normB = norm ( b ) ;
double r = eps + 1 . 0 ;
while ( r > eps )
{
/ ∗ ∗∗∗
∗ Gaussova−S e i d e lo v a i t e r a c e
∗/
f o r ( i n t i = 0 ; i < n ; i ++ )
{
double s = 0 . 0 ;
f o r ( i n t j = 0 ; j < n ; j ++ )
i f ( j != i )
s = s + A[ i ] [ j ] ∗ x [ j ] ;
x [ i ] = ( b [ i ] − s ) / A[ i ] [ i ] ;
}
/ ∗ ∗∗∗
∗ Vypocet r e z i d u a
∗/
r = 0.0;
f o r ( i n t i = 0 ; i < n ; i ++ )
{
double Ax_i = 0 . 0 ;
f o r ( i n t j = 0 ; j < n ; j ++ )
Ax_i = Ax_i + A [ i ] [ j ] ∗ x [ j ] ;
r = r + ( Ax_i − b _ i ) ∗ ( Ax_i − b _ i ) ;
}
r = s q r t ( r ) / normB ;
}
}
40 / 62
Tomáš
Oberhuber
Iterativní
metody
Gaussova-Seidelova metoda numerická analýza
• pro numerickou analýzu lze tuto metodu maticoveˇ
zapsat takto
~x (k +1) =
D~x (k +1) − L~x (k +1)
~ (k +1)
x
=
=
~x (k +1) =
~x (k +1) =
D−1 ~b + L~x (k +1) + R~x (k )
~b + R~x (k )
i
h
(D − L)−1 ~b + R~x (k )
(D − L)−1 R~x (k ) + (D − L)−1~b
h
i
I − (D − L)−1A ~x (k ) +
(D − L)−1~b
h
i
41 / 62
Tomáš
Oberhuber
Iterativní
metody
Gaussova-Seidelova metoda numerická analýza
• vidíme tedy, že
B
~c =
H
D − L)−1R
(D − L)−1~b
(D − L)−1
= (
=
42 / 62
Tomáš
Oberhuber
Iterativní
metody
Gaussova-Seidelova metoda numerická analýza
Theorem 23
Nutnou a postaˇcující podmínkou pro to, aby
Gaussova-Seidelova metoda konvergovala pro libovolné
~x (0) je
ρ ( − )−1
< 1.
D L R
Postaˇcující podmínkou je
−1 (
−
)
< 1,
D L R
ˇ
v nekteré
maticové normeˇ k·k.
43 / 62
Tomáš
Oberhuber
Gaussova-Seidelova metoda numerická analýza
Iterativní
metody
Theorem 24
A
Má-li matice pˇrevládající diagonálu, pak
Gaussova-Seidelova metoda konverguje k rˇešení soustavy
~x = ~b pro libovolnou volbu ~x (0) .
A
Theorem 25
A
Je-li matice hermitovská a pozitivneˇ definitní, pak
Gaussova-Seidelova metoda konverguje k rˇešení soustavy
~x = ~b pro libovolnou volbu ~x (0) . Konvergence je navíc
monotónní vzhledem k vektorové normeˇ k·kA .
A
Dukaz.
˚
ˇ 14, kde W = D − L a
Plyne z vety
W + W∗ = 2D − L − L∗ = D + A > A > 0.
44 / 62
Tomáš
Oberhuber
Gaussova-Seidelova metoda numerická analýza
Iterativní
metody
Remark 26
A
D A
Bud’ PD matice taková, že 2 − není PD. Potom
Jacobiho metoda nekonverguje, ale Gaussova-Seidelova
ano. Gaussova-Seidelova metoda je považována ze lepší
ˇ
metodu, než Jacobiho, i když v nekterých
pˇrípadech
Jacobiho metoda konverguje lépe nebo konverguje tehdy,
když Gaussova-Seidelova metoda nekonverguje.
45 / 62
Tomáš
Oberhuber
Gaussova-Seidelova metoda pˇríklad
Iterativní
metody
Example 27
A




 
−2 1
0
−1
1
=  1 −2 1  , ~b =  0  , ~x =  1  ,
0
1 −2
−1
1
Jacobi
Gauss − Seidel
~x k
~r k
~x k
~r k
0
(0, 0, 0) (1, 0, 1)
(0, 0, 0)
(1, 0, 1)
1
( 21 , 0, 12 ) (0, 1, 0)
( 12 , 41 , 58 )
( 34 , 58 , 0)
13
3
2
( 21 , 12 , 12 ) ( 12 , 0, 21 )
( 58 , 85 , 16
)
( 38 , 16
, 0)
3 1 3
1
13 13 29
3 3
3
( 4 , 2 , 4 ) (0, 2 , 0)
( 16 , 16 , 32 )
( 16 , 32 , 0)
29 61
3 3
4
( 43 , 34 , 34 ) ( 14 , 0, 41 )
( 29
( 32
, 64 , 0)
32 , 32 , 64 )
7 3 7
1
61 61 125
3
3
5
( 8 , 4 , 8 ) (0, 4 , 0)
( 64 , 64 , 128 ) ( 64 , 128
, 0)
7 7 7
1
125 125 253
3
1
3
6
( 8 , 8 , 8 ) ( 8 , 0, 8 ) ( 128 , 128 , 256 ) ( 128 , 256
, 0)
...
...
...
...
...
k
46 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda
Iterativní
metody
• David M. Young, Jr., 1950
• je známa zejména pod zkratkou SOR = Successive
Over Relaxation method
• jde o modifikaci Gaussovy-Seidelovy metody
• metoda funguje velmi dobˇre zejména na lineární
systémy pocházející z metody koneˇcných diferencí pro
ˇrešení parabolických nebo eliptických parciálních
diferenciálních rovnic
47 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda
Iterativní
metody
• Gaussovu-Seidelovu metodu mužeme
˚
pˇrepsat do tvaru
~xi(k +1) = ~xi(k ) + ∆~xi(k ) ,
kde
(k )
∆~xi
=
1 (k +1)
(k +1)
(k +1)
−ai1~x1
− ai2~x2
− . . . − ai,i−1~xi−1
aii
(k )
(k )
(k )
−aii ~x − ai,i+1~x
− . . . − ain~xn + ~bi
i
i+1
48 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda
Iterativní
metody
• SOR metoda je pak definována jako
~xi(k +1) = ~xi(k ) + ω∆~xi(k ) ,
• ω je tzv. relaxaˇcní parametr
• je-li ω = 1, dostáváme Gaussovu-Seidelovu metodu
• podobnou modifikaci lze snadno provést i pro Jacobiho
metodu
49 / 62
Tomáš
Oberhuber
Iterativní
metody
Super-relaxaˇcní metoda
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
31
32
33
34
35
bool sorMethod ( const M a t r i x& A ,
const V e c t o r& b ,
V e c t o r& x ,
double eps ,
double omega )
{
double normB = norm ( b ) ;
double r = eps + 1 . 0 ;
while ( r > eps )
{
/ ∗ ∗∗∗
∗ SOR i t e r a c e
∗/
f o r ( i n t i = 0 ; i < n ; i ++ )
{
double s = 0 . 0 ;
f o r ( i n t j = 0 ; j < n ; j ++ )
s = s + A[ i ] [ j ] ∗ x [ j ] ;
x [ i ] += omega ∗ ( b [ i ] − s ) / A [ i ] [
}
/ ∗ ∗∗∗
∗ Vypocet r e z i d u a
∗/
r = 0.0;
f o r ( i n t i = 0 ; i < n ; i ++ )
{
double Ax_i = 0 . 0 ;
f o r ( i n t j = 0 ; j < n ; j ++ )
Ax_i = Ax_i + A [ i ] [ j ] ∗ x [ j ] ;
r = r + ( Ax_i − b _ i ) ∗ ( Ax_i − b _ i ) ;
}
r = s q r t ( r ) / normB ;
i
];
}
}
50 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda numerická analýza
Iterativní
metody
• maticoveˇ lze metodu zapsat jako
~x (k +1) = ~x (k ) + ω
D−1 L~x (k +1) − D~x (k ) + R~x (k ) + ~b
tj.
I
~x (k +1) = ( − ω(
tj.
Bω
Hω
~cω
D − ωL)−1A)~x (k ) + ω(D − ωL)−1~b,
I D L A
D L
D + ω R) ,
D L
D L
=
− ω( − ω )−1
= ( − ω )−1 ((1 − ω)
= ω( − ω )−1 ,
= ω( − ω )−1~b.
51 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda numerická analýza
Iterativní
metody
Theorem 28
Nutnou a postaˇcující podmínkou pro to, aby SOR metoda
konvergovala pro libovolné ~x (0) je
ρ(
Bω ) < 1.
Postaˇcující podmínkou je
k
Bω k < 1,
ˇ
v nekteré
maticové normeˇ k·k.
52 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda numerická analýza
Iterativní
metody
Theorem 29
Pro libovolné ω ∈
R platí
ρ (Bω ) ≥ |ω − 1| ,
a tedy SOR metoda nekonverguje pro ω ≤ 0 nebo ω ≥ 2.
Theorem 30
A
Má-li matice pˇrevládající diagonálu, pak SOR metoda
konverguje k rˇešení soustavy ~x = ~b pro libovolnou volbu
~x (0) pokud je 0 < ω ≤ 1.
A
Dukaz.
˚
Zatím neznám.
53 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda numerická analýza
Iterativní
metody
Theorem 31
A
(Ostrowski): Je-li matice hermitovská a pozitivneˇ definitní,
pak SOR metoda konverguje k rˇešení soustavy ~x = ~b pro
libovolnou volbu ~x (0) práveˇ když 0 < ω < 2. Konvergence je
navíc monotónní vzhledem k vektorové normeˇ k·kA .
A
Dukaz.
˚
ˇ 14:
Z vety
W
1
=
ω
D−L⇒W+W
∗
=
2
−1 D+A>A>0
ω
54 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda numerická analýza
Iterativní
metody
Remark 32
Budeme zkoumat konvergenci SOR metody v závislosti na
parametru ω. Kromeˇ konvergence samotné nás bude
zajímat i její rychlost, tj. budeme chtít najít ω0 takové, aby
bylo ρ ( ω ) resp. k ω k co nejmenší.
B
B
Remark 33
Vliv parametru ω na rychlost konvergence budeme
analyzovat pro speciální typ matic tzv. dvoucyklických,
shodneˇ uspoˇrádaných. Dá se ukázat, že tyto matice
vznikají práveˇ pˇri rˇešení eliptických nebo parabolických
rovnic pomocí metody koneˇcných diferencí.
55 / 62
Tomáš
Oberhuber
Iterativní
metody
Super-relaxaˇcní metoda numerická analýza
Definition 34
I
ˇ maticí ij pro (1 ≤ i, j ≤ n)
Elementární permutacní
nazveme takovou cˇ tvercovou matici n-tého 0ˇrádu, která má
všechny diagonální prvky rovny 1, kromeˇ prvku˚ na i-tém a
j-tém ˇrádku, které jsou rovny 0. Mimodiagonální prvky jsou
nulové kromeˇ prvku˚ na pozicích (i, j) a (j, i), které jsou
rovny 1.

Iij











=











1
..






















.
1
0 ... ... ... 1
..
..
. 1
.
..
..
..
.
.
.
..
.
.
1 ..
1 ... ... ... 0
1
..
.
1
ˇ maticí nazveme každou matici, kterou lze
Permutacní
vyjádˇrit jako souˇcin libovolného poˇctu elementárních
permutaˇcních matic.
56 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda numerická analýza
Iterativní
metody
• násobením libovolné obdélnikové matice elementární
permutaˇcní maticí zleva reps. zprava odpovídá
prohození i-tého a j-tého ˇrádku resp. sloupce
• snadno vidíme, že
Iij = ITij = I−1
ij
a tedy elementární permutaˇcní matice je symetrická a
ortogonální
• ... a tedy i permutaˇcní matice je symetrická a
ortogonální
57 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda numerická analýza
Iterativní
metody
Definition 35
C
ˇ
Ctvercová
matice se nazývá slabeˇ cyklická s indexem
2, jestliže existuje permutaˇcní matice taková, že matice
T je blokového tvaru
θ
1
T
=
θ
2
P
PCP
PCP
M
M
58 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda numerická analýza
Iterativní
metody
Definition 36
A D L R
A
Necht’ = − − . Matice se nazývá dvoucyklická,
je-li odpovídající Jacobiho matice J = −1 ( + ) slabeˇ
cyklická s indexem 2.
Definition 37
B
D L R
A
Dvoucyklická matice se nazývá shodneˇ uspoˇrádaná,
jestliže vlastní cˇ íslo matice
α
D−1L + α1 D−1R
nazávisí na volbeˇ cˇ ísla α 6= 0.
59 / 62
Tomáš
Oberhuber
Super-relaxaˇcní metoda numerická analýza
Iterativní
metody
Theorem 38
A
Necht’ matice je dvoucyklická, shodneˇ uspoˇrádaná.
Necht’ λ 6= 0 je vlastní cˇ íslo matice ω a ω 6= 0. Necht’ µ
ˇ
splnuje
rovnici
(λ + ω − 1)2 = ω 2 µ2 λ.
B
B
Pak µ je vlastní cˇ íslo Jacobiho matice J a naopak, je-li µ
vlastní cˇ íslo matice J , pak pˇríslušné λ je vlastním cˇ íslem
matice ω . Navíc platí, že pro
B
B
ωopt =
2
q
1 + 1 − ρ(
BJ )2
,
B
nabývá ρ ( ω ) své minimum a SOR metoda tedy
konverguje nejrychleji.
60 / 62
Tomáš
Oberhuber
Shrnutí
Iterativní
metody
Podmínky na matici
A pro konvergenci metod:
s pˇrevládající diagonálou
Jacobi
Gauss-Seidel
SOR
!
!
!
hermitovská a ...
2 > >0
A>0
A>0
D A
61 / 62
Tomáš
Oberhuber
Shrnutí
Iterativní
metody
• porovnání pˇrímých a iteracních
ˇ
metod pro rˇešení
lineárních soustav rovnic
• podle cˇ eho byste se rozhodli, kterou tˇrídu metod
použijete?
• konvergence stacionárních metod, apriorní a
aposteriorní odhady chyb
• Jacobiho, Gaussova-Seidelova a SOR metoda
• odvození, maticový tvar (matice
B, H), konvergence
62 / 62
Download

Iterativní metody