Faculty of Nuclear Sciences and Physical Engineering Czech Technical University in Prague

Tomáš Oberhuber Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda

Tomáš Oberhuber Faculty of Nuclear Sciences and Physical Engineering Czech Technical University in Prague

GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

1 / 65

Tomáš Oberhuber Iterativní metody Stacionární metody

1 Iterativní metody 2 Stacionární metody

Metoda postupných aproximací

3 Metoda postupných aproximací

ˇ Pˇredpodmínení

ˇ 4 Pˇredpodmínení

Richardsonovy iterace Jacobiho metoda

5 Richardsonovy iterace

GaussovaSeidelova metoda

6 Jacobiho metoda

Superrelaxaˇcní metoda

7 Gaussova-Seidelova metoda

Shrnutí a porovnání pˇrímých a iterativních metod

8 Super-relaxaˇcní metoda 9 Shrnutí a porovnání pˇrímých a iterativních metod 2 / 65

Tomáš Oberhuber

Iterativní metody

Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

A

Budeme se zabývat metodami pro ˇrešení úlohy ~x = ~b s ˇ regulární maticí . Chceme najít efektivnejší metody, než je GEM. (k ) ∞ • iterativní metody generují posloupnost ~ , která x k =1 konverguje k ˇrešení soustavy

A

lim ~x (k ) = ~x ∗

k →∞

• prakticky tak (obecne) ˇ nikdy nedostaneme pˇresné

ˇrešení, ale mužeme ˇ r) libovolneˇ pˇresnou ˚ dostat (témeˇ aproximaci

2 / 65

Tomáš Oberhuber

Iterativní metody obecneˇ

Iterativní metody Stacionární metody Metoda postupných aproximací

• obecneˇ platí, že volíme libovolné ~ x (0) a dále

napoˇcítáváme ~x (k +1) =

ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

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 ∗ =

3 / 65

Tomáš Oberhuber Iterativní metody

Iterativní metody obecneˇ Theorem 1 Iterativní metoda tvaru

Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení

~x (k +1) = ˇ splnující ~x ∗ =

lim

k →∞

Jacobiho metoda

Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

B(k )~x ∗ + ~c (k ),

konverguje k ~x ∗ pro libovolné ~x (0) práveˇ když,

Richardsonovy iterace


B(k )~x (k ) + ~c (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 ∗ 4 / 65


Citlivost a stabilita iterativních metod

Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda

• 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


5 / 65

Tomáš Oberhuber

Iterativní metody obecneˇ

Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda

ˇ na: Iterativní metody pro ˇrešení lineárních soustav 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.


6 / 65


Stacionární iterativní metody Theorem 2 Metoda daná vztahem ~x (k +1) =

Metoda postupných aproximací ˇ Pˇredpodmínení

ˇ a splnující ~x ∗ =



B~x ∗ + ~c ,

konverguje k ~x ∗ pro libovolné ~x (0) práveˇ když platí,

GaussovaSeidelova metoda Superrelaxaˇcní metoda

B~x (k ) + ~c

lim

k →∞

Bk = θ.

Dukaz. ˚ ˇ 1a Plyne z vety lim

k →∞

B(k −1) . . . B(0) = klim Bk . →∞ 7 / 65

Tomáš Oberhuber

Stacionární iterativní metody

Iterativní metody Stacionární metody

Theorem 3


Metoda daná vztahem ~x (k +1) =

ˇ Pˇredpodmínení Richardsonovy iterace


Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

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 . 8 / 65

Tomáš Oberhuber



Theorem 4


Postaˇcující podmínkou pro to, aby metoda daná vztahem ~x (k +1) =




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 . 9 / 65

Tomáš Oberhuber Iterativní metody Stacionární metody Metoda postupných aproximací

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 ) ,


pˇri použití souhlasné maticové normy s danou vektorovou normou.


Remark 6

Jacobiho metoda

A A I B B

A

Zde jsme použili definici rezidua v k -té iteraci ~r (k ) =

A~x (k ) − ~b. 10 / 65

Tomáš Oberhuber



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


11 / 65

Tomáš Oberhuber Iterativní metody Stacionární metody Metoda postupných aproximací

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) =

Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod


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. 12 / 65

Tomáš Oberhuber



• 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


13 / 65

Tomáš Oberhuber Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda

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

I A) ~x (k ) + ~b

• vidíme, že metoda splnuje ˇ podmínku

I A) ~x ∗ + ~b


A

~x (k +1) = ~x (k ) −~r (k ) = ~x (k ) − ~x (k ) + ~b = ( −

~x ∗ = ( − • je tedy

B

=

⇔

A~x ∗ = ~b

I − A,

~c = ~b

14 / 65


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

I A) ~x (k ) + ~b,

~x (k +1) = ( −


A~x = ~b (s

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) . 15 / 65

Tomáš Oberhuber


Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda

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.


16 / 65


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 >

ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda

Dukaz. ˚

A I A

A

I

σ ( ) ∈ (0, 2) ⇔ 2 > A > 0.


I A

Je-li hermitovská a má-li být ρ ( − ) < 1, pak musí být ˇ σ ( − ) ⊂ (−1, 1) a podle pˇredchozí vety

Remark 12 Takových matic ale moc není. 17 / 65

Tomáš Oberhuber

ˇ Pˇredpodmínená metoda postupných aproximací


• 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í ˇ


18 / 65

Tomáš Oberhuber



• 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 ) −

Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace

=

Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda

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 ∗

( −


19 / 65



Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda

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 ˇ


20 / 65

Tomáš Oberhuber



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


GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

H

I HA) ~x (k ) + H~b,

~x (k +1) = ( −


A

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) . 21 / 65

Tomáš Oberhuber



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 22 / 65

Tomáš Oberhuber



• ze vztahu

pˇri volbeˇ


~x (k +1)

H=A

Jacobiho metoda

je

−

HA

~x (k )

HA~x (k ) + H~b A−1A~x (k ) + A−1~b

= ~x (k ) − ~x (k ) + ~x ∗ = ~x ∗



=

−1

~x (k )

~x (k +1) = ~x (k ) − = ~x (k ) −



H + H~b dostáváme, že

• nyní vyvstává otázka, jak volit matici

• 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

23 / 65

Tomáš Oberhuber



• tato metoda je dána pˇredpisem

A


~x (k +1) = ~x (k ) − θ( ~x (k ) − ~b),

ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda

pro θ ∈

R

• parametru θ se ˇríká relaxacní ˇ parametr • snadno vidíme, že

H B

I

= θ = −θ ~c = θ~b

I

A


24 / 65


Richardsonovy iterace Theorem 16

A

Je-li hermitovská a PD, pak metoda Richardsonových iterací konverguje práveˇ když

I A > 0.

2 > θ

Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace

Konvergence je navíc monotónní vzhledem k vektorové normeˇ k·kA .

Jacobiho metoda

Dukaz. ˚


ˇ 14, kde Dukaz ˚ postaˇcující podmínky plyne z vety = Že jde o nutnou podmínku lze ukázata podobneˇ jako u metody postupných aproximací.


W

Remark 17 Pro hermitovskou a PD matici θ<

1 θ

I.

A metoda konverguje pro

2 . ρ( )

A

25 / 65

Tomáš Oberhuber

Jacobiho metoda


• 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. 26 / 65

Tomáš Oberhuber 1 2 Iterativní 3 metody 4 5 Stacionární 6 metody 7 8 Metoda 9 postupných 10 11 aproximací 12 ˇ 13 Pˇredpodmínení 14 Richardsonovy 15 16 iterace 17 18 Jacobiho 19 metoda 20 21 Gaussova22 Seidelova 23 metoda 24 25 26 Super27 relaxaˇcní 28 metoda 29 30 Shrnutí a 31 porovnání 32 pˇrímých a 33 iterativních 34 35 metod 36 37

Jacobiho metoda 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 ; } }

27 / 65

Tomáš Oberhuber

Jacobiho metoda - numerická analýza

Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení

• matici

Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda

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


28 / 65




• 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


29 / 65




• vidíme tedy, že

B

=

~c =

H

=

I − D−1A D−1~b D−1


30 / 65



Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

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.

31 / 65




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.


32 / 65

Tomáš Oberhuber



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í. 33 / 65


Jacobiho metoda - pˇríklad Example 22


A

     1 −2 1 0 −1 =  1 −2 1  , ~b =  0  , ~x =  1  , 1 0 1 −2 −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 ) ... 34 / 65

Tomáš Oberhuber

Gaussova-Seidelova metoda


• roku 1823 jí popsal Carl Friedrich Gauss v jednom

dopise svému studentovi Ch.L.Gerlingovi • roku 1874 jí publikoval Phillip Ludwig von Seidel


35 / 65

Tomáš Oberhuber Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda

Gaussova-Seidelova metoda • 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.


36 / 65

Tomáš Oberhuber 1 2 Iterativní 3 metody 4 5 Stacionární 6 metody 7 8 Metoda 9 postupných 10 aproximací 11 ˇ 12 Pˇredpodmínení 13 Richardsonovy 14 15 iterace 16 17 Jacobiho 18 metoda 19 20 Gaussova21 Seidelova 22 metoda 23 24 Super25 relaxaˇcní 26 metoda 27 28 Shrnutí a 29 porovnání 30 31 pˇrímých a 32 iterativních 33 metod 34 35

Gaussova-Seidelova metoda 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 ; } }

37 / 65

Tomáš Oberhuber Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

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

38 / 65


Gaussova-Seidelova metoda numerická analýza


• vidíme tedy, že

B

~c =

H

D − L)−1R (D − L)−1~b (D − L)−1

= ( =


39 / 65




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.

40 / 65

Tomáš Oberhuber



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. 41 / 65

Tomáš Oberhuber



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.


42 / 65

Tomáš Oberhuber

Gaussova-Seidelova metoda pˇríklad

Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení

Example 27

A



     −1 1 −2 1 0 =  1 −2 1  , ~b =  0  , ~x =  1  , 0 1 −2 −1 1

Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

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 ) , 0) 2 ( 21 , 12 , 12 ) ( 12 , 0, 21 ) ( 58 , 85 , 16 ( 38 , 16 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 , 0) 5 ( 8 , 4 , 8 ) (0, 4 , 0) ( 64 , 64 , 128 ) ( 64 , 128 7 7 7 1 1 125 125 253 3 3 6 ( 8 , 8 , 8 ) ( 8 , 0, 8 ) ( 128 , 128 , 256 ) ( 128 , 256 , 0) ... ... ... ... ... k

43 / 65

Tomáš Oberhuber

Super-relaxaˇcní metoda


• 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

44 / 65

Tomáš Oberhuber



• 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


45 / 65

Tomáš Oberhuber



• 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


46 / 65

Tomáš Oberhuber 1 2 Iterativní 3 metody 4 5 Stacionární 6 metody 7 8 Metoda 9 postupných 10 aproximací 11 12 ˇ 13 Pˇredpodmínení 14 Richardsonovy 15 iterace 16 17 Jacobiho 18 19 metoda 20 21 Gaussova22 Seidelova 23 metoda 24 25 Super26 relaxaˇcní 27 metoda 28 29 Shrnutí a 30 porovnání 31 pˇrímých a 32 iterativních 33 34 metod 35 36

Super-relaxaˇcní metoda 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 ++ ) i f ( j != i ) 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

];

} }

47 / 65

Tomáš Oberhuber

Super-relaxaˇcní metoda numerická analýza


• 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.

48 / 65

Tomáš Oberhuber



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.

49 / 65

Tomáš Oberhuber



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. 50 / 65

Tomáš Oberhuber



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: Dukaz ˚ postaˇcující podmínky z vety 1 2 ∗ = − ⇒ + = −1 D+A>A>0 ω ω

W

D L W W

51 / 65

Tomáš Oberhuber



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í.

52 / 65

Tomáš Oberhuber Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení

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 ˇ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. 53 / 65

Tomáš Oberhuber


Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace

• 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í

54 / 65

Tomáš Oberhuber



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


55 / 65

Tomáš Oberhuber



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.

56 / 65

Tomáš Oberhuber



Theorem 38

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


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 λ.

ωopt =

2 q 1 + 1 − ρ(

BJ )2

,

B

nabývá ρ ( ω ) své minimum a SOR metoda tedy konverguje nejrychleji. 57 / 65

Tomáš Oberhuber

Shrnutí


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


58 / 65

Tomáš Oberhuber

Výhody iterativních metod


• jsou-li vhodneˇ použity, jsou mnohem rychlejší než

pˇrímé metody založené na GEM • GEM má složitost n3 , pro velké matice je témeˇ ˇr

nepoužitelná • 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 • v praxi se cˇ asto podaˇrí, že dobrou aproximaci ˇrešení získáme po pár iteracích nehledeˇ na velikost matice


59 / 65


Výhody iterativních metod Definition 39 ˇ ˇ Rídkou maticí nazveme takovou, u které je vetšina prvku˚ nulových.


Zdroj: Matrix market 60 / 65


Výhody iterativních metod Iterativní metody dokáží lépe využít vlastností ˇrídkých matic.

Stacionární metody Metoda postupných aproximací

• pˇri šikovném uložení ˇrídkých matic muže ˚ mít násobení


• modifikace GEM pro obecné ˇrídké matice je možná, ale

Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda

matice a vektoru ješteˇ menší nároˇcnost, než n2 algoritmicky velmi složitá • pˇri aplikaci GEM na ˇrídkou matici muže ˚ vzniknout

mnoho nových nenulových prvku˚ a narust ˚ tak nároky ˇ na pamet’


• iterativní metody nemení ˇ samotnou maticí


• v nekterých ˇ situacích ani nemáme matici

jsou snažší na implementaci

A, díky tomu

A

explicitneˇ ˇ pak je GEM nemožné použít uloženou v pameti,

61 / 65

Tomáš Oberhuber

Výhody iterativních metod


• iterativní metody dokážou velice efektivneˇ využít

znalosti pˇribližného ˇrešení, konvergují pak velice rychle • v pˇrípadeˇ GEMu je nám jakýkoliv odhad ˇrešení k

niˇcemu


62 / 65

Tomáš Oberhuber

Nevýhody iterativních metod


• vetšina ˇ z nich nedá ani teoreticky pˇresné ˇrešení po

koneˇcném poˇctu kroku˚ • 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ˇ • vyšetˇrování konvergence je cˇ asto velice složité • modifikovaná GEM funguje pro libovolnou regulární

matici


63 / 65

Tomáš Oberhuber Iterativní metody Stacionární metody Metoda postupných aproximací ˇ Pˇredpodmínení Richardsonovy iterace Jacobiho metoda GaussovaSeidelova metoda Superrelaxaˇcní metoda Shrnutí a porovnání pˇrímých a iterativních metod

Výpoˇcetní pˇríklad Pomocí programu stationary-solver otestujte jednotlivé metody. Argumenty pˇríkazové ˇrádky jsou • -input-file nebo -i - soubor s maticí soustavy • -method - stacionární metoda

(richardson/jacobi/sor) • -relaxation - relaxaˇcní parametr metody • -matrix-format - maticový formát

(dense/ellpack) • -initial-value - nastavení vektoru ~ x (0) • -max-iterations - maximální poˇcet iterací • -convergence-residue - zastavující kritérium • výhoda stacionárních metod oproti GEM jsou videt ˇ na

PD matici gr_30_30.mtx • špatná konvergence na PD matici bcsstk01.mtx • indefinitní matice pores_1.mtx 64 / 65

Tomáš Oberhuber

Shrnutí


• konvergence stacionárních metod, apriorní a

aposteriorní odhady chyb • Jacobiho, Gaussova-Seidelova a SOR metoda • odvození, maticový tvar (matice

B, H), konvergence

• 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?


65 / 65

Faculty of Nuclear Sciences and Physical Engineering Czech Technical University in Prague

Recommend Documents