Numerické metody
9. přednáška, 19. dubna 2018
Jiří Zelinka
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 1 / 24
Opakování
Iterační metody řešení systémů lineárních rovnic
Systém
Ax = b
převedeme na
x = Tx + g
xk+1
= Txk
+ g, k = 0, 1, . . .
Hlavní věta o konvergenci iteračního procesu
Posloupnost xk ∞
k=0
určená iteračním procesem x = Tx + g
konverguje pro každou počáteční aproximaci x0
∈ Rn
právě
tehdy, když ρ(T) < 1, přičemž
lim
k→∞
xk
= x∗
, x∗
= Tx∗
+ g
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 2 / 24
Jacobiova iterační metoda
Ax = b, A =





a11 · · · · · · a1n
a21 a2n
...
...
an1 · · · · · · ann





Matici A zapišme ve tvaru
A = D + L + U,
kde
D =



a11 0
...
0 ann


 ,
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 3 / 24
L =





0 0
a21
...
...
...
...
an1 · · · an,n−1 0





,
U =





0 a12 · · · a1n
... ...
...
... an−1,n
0 0





.
D je diagonální matice, L je dolní trojúhelníková matice
s nulami na diagonále a U je horní trojúhelníková matice
s nulami na diagonále.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 4 / 24
Ax = (D + L + U)x = b
Dx = −(L + U)x + b.
Pokud aii = 0, i = 1, . . . , n, je matice D regulární
a z předchozí rovnice lze vypočítat
x = −D−1
(L + U)x + D−1
b.
D−1
=



1
a11
0
...
0 1
ann


 ,
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 5 / 24
Maticový tvar Jacobiovy iterační metody
Jacobiova iterační matice: TJ = −D−1
(L + U)
xk+1
= TJxk
+ D−1
b,
TJ = (tij ), tij = −
aij
aii
pro i = j, tii = 0 pro i = 1, . . . , n.
TJ =












0 −
a12
a11
· · · −
a1n
a11
−
a21
a22
0 −
a2n
a22
...
...
...
−
an1
ann
−
an2
ann
· · · 0












, D−1
b =












b1
a11
b2
a22
...
bn
ann












.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 6 / 24
Realizace výpočtu:
Z první rovnice vypočteme x1:
a11x1+a12x2+· · ·+a1nxn = b1 ⇒ x1 =
1
a11
(b1−a12x2−. . .−a1nxn)
xk+1
1 =
1
a11
(b1 − a12xk
2 − . . . − a1nxk
n ),
z druhé rovnice vypočteme x2:
xk+1
2 =
1
a22
(b2 − a21xk
1 − a23xk
3 − . . . − a1nxk
n ),
obecně z i-té rovnice vypočteme xi :
xk+1
i = −
n
j=1
j=i
aij
aii
xk
j +
bi
aii
,
až z n-té rovnice vypočteme xn, a na pravé straně takto
získaného systému jsou prvky matice TJ.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 7 / 24
Věta o konvergenci Jacobiovy iterační metody:
Posloupnost xk ∞
k=0
generovaná metodou
xk+1
= TJxk
+ D−1
b konverguje pro každou počáteční
aproximaci x0
∈ Rn
právě tehdy, když (TJ) < 1.
Odhad chyby:
x∗
− xk
∞ ≤
TJ
k
∞
1 − T ∞
x1
− x0
∞.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 8 / 24
Silné řádkové sumační kriterium:
Nechť matice A je ryze řádkově diagonálně dominantní, tj.
|aii | >
n
j=1
j=i
|aij |, i = 1, . . . , n.
Pak Jacobiova iterační metoda konverguje pro každou
počáteční aproximaci x0
∈ Rn
.
Silné sloupcové sumační kriterium:
Nechť matice A je ryze sloupcově diagonálně dominantní, tj.
|akk| >
n
i=1
i =k
|aik|, k = 1, . . . , n.
Pak Jacobiova iterační metoda konverguje pro každou
počáteční aproximaci x0
∈ Rn
.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 9 / 24
Gaussova-Seidelova iterační metoda
Z první rovnice vypočteme x1:
a11x1+a12x2+· · ·+a1nxn = b1 ⇒ x1 =
1
a11
(b1−a12x2−. . .−a1nxn)
xk+1
1 =
1
a11
(b1 − a12xk
2 − . . . − a1nxk
n ),
z druhé rovnice vypočteme x2, pro x1 použijme novou iteraci:
xk+1
2 =
1
a22
(b2 − a21xk+1
1 − a23xk
3 − . . . − a1nxk
n ),
ze třetí rovnice vypočteme x3, pro x1 a x2 použijme novou
iteraci:
xk+1
3 =
1
a33
(b3 − a31xk+1
1 − a32xk+1
2 − a34xk
4 . . . − a1nxk
n ),
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 10 / 24
Obecně
xk+1
i = −
i−1
j=1
aij
aii
xk+1
j −
n
j=i+1
aij
aii
xk
j +
bi
aii
, i = 1, . . . , n.
Maticový zápis:
Ax = b ⇒ (D + L + U)x = b
(D + L)x = −Ux + b.
aii = 0, i = 1, . . . , n, ⇒ matice D + L je regulární a
x = −(D + L)−1
Ux + (D − L)−1
b.
Položme TG = −(D + L)−1
U, Gaussova-Seidelova iterační
metoda je tvaru
xk+1
= TG xk
+ g, g = (D + L)−1
b.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 11 / 24
Věta
Posloupnost xk ∞
k=0
generovaná Gaussovou-Seidelovou
iterační metodou xk+1
= TG xk
+ g. konverguje pro každou
počáteční aproximaci x0
∈ Rn
právě tehdy, když (TG ) < 1.
Silné řádkové sumační kriterium:
Nechť matice A je ryze řádkově diagonálně dominantní, tj.
|aii | >
n
j=1
j=i
|aij |, i = 1, . . . , n.
Pak Gaussova-Seidelova iterační metoda konverguje pro
každou počáteční aproximaci x0
∈ Rn
.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 12 / 24
Silné sloupcové sumační kriterium:
Nechť matice A je ryze sloupcově diagonálně dominantní, tj.
|akk| >
n
i=1
i =k
|aik|, k = 1, . . . , n.
Pak Gaussova-Seidelova iterační metoda konverguje pro
každou počáteční aproximaci x0
∈ Rn
.
Věta
Nechť A je pozitivně definitní matice. Pak Gaussova-Seidelova
metoda konverguje pro každou počáteční aproximaci.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 13 / 24
Věta (Stein-Rosenberg)
Nechť pro prvky matice A platí aij ≤ 0 pro všechna i = j
a aii > 0, i = 1, . . . , n. Pak platí právě jedno z následujících
tvrzení:
0 < (TG ) < (TJ) < 1
1 < (TJ) < (TG )
(TJ) = (TG ) = 0
(TJ) = (TG ) = 1.
To znamená, že konvergují-li obě metody, Gaussova-Seidelova
metoda konverguje rychleji.
Konec opakování
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 14 / 24
Relaxační metody
Modifikace Gaussovy–Seidelovy metody, ω – relaxační
parametr
xk+1
i = (1 − ω)xk
i +
ω
aii
bi −
i−1
j=1
aij xk+1
j −
n
j=i+1
aij xk
j .
Relaxační metodu lze maticově zapsat takto
xk+1
= (D − ωL)−1
[(1 − ω)D + ωU]xk
+ ω(D − ωL)−1
b
Tω = (D − ωL)−1
[(1 − ω)D + ωU]
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 15 / 24
Hodnoty parametru ω:
Pro 0 < ω < 1 se iterační metody nazývají metodami dolní
relaxace. Tyto metody jsou vhodné v případě,
že Gaussova-Seidelova metoda nekonverguje.
Pro ω = 1 je relaxační metoda totožná s GaussovouSeidelovou
metodou.
Pro 1 < ω se metody nazývají metodami horní relaxace,
nebo častěji SOR metodami (SOR
= Successive Over-Relaxation). Tyto metody
lze užít ke zrychlení konvergence GaussovySeidelovy
metody.
Věta (Kahan).
Nechť aii = 0, i = 1, . . . , n. Pak
(Tω) ≥ |ω − 1|.
Důsledek: Má smysl uvažovat jen ω ∈ (0, 2).
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 16 / 24
Věta (Ostrowski-Reich).
Pro pozitivně definitní matici A platí (Tω) < 1 pro všechna
ω ∈ (0, 2).
Třídiagonální matice:
A = (aij )n
i,j=1, aij = 0, pro |i − j| > 1
Věta.
Nechť A je třídiagonální pozitivně definitní matice. Pak
(TG ) = 2
(TJ) < 1 a optimální hodnota relaxačního
parametru je dána vztahem
ω = ωopt =
2
1 + 1 − 2(TJ)
.
Při této volbě je (Tω) = |1 − ω|.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 17 / 24
Cykly v iteračních metodách
Například pro systémy
x1 + kx2 = b1
x1 − kx2 = b2
Jacobiova metoda: cyklus délky 4
Gaussova–Seidelova metoda: cyklus délky 2
Relaxační metody: cykly různých délek
x1 + x2 = 1
qx1 + x2 = 1.
Pro ω = 2:
T2 =
−1 −2
2q 4q − 1
, λ1,2 = cos ϕ±i sin ϕ, q =
1
2
(1+cos ϕ)
ϕ = 2πl/p, 0 < l < p/2 ⇒ existuje cyklus délky p.
Body cyklu leží na elipse se středem v hledaném řešení.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 18 / 24
Iterační metody pro systémy nelineárních rovnic
f1(x1, . . . , xm) = 0
...
fm(x1, . . . , xm) = 0
Kořen systému: uspořádaná m-tice reálných čísel
ξ = (ξ1, . . . , ξm),
která tomuto systému vyhovuje.
Vektorový tvar:
F(x) = o, x ∈ Rm
, o = (0, . . . , 0)T
∈ Rm
.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 19 / 24
Systém převedeme na ekvivalentní rovnici
x = G(x), x ∈ Rm
neboli
x1 = g1(x1, . . . , xm)
...
xm = gm(x1, . . . , xm)
a budeme hledat pevný bod zobrazení G : Rm
→ Rm
.
Definujme nyní v prostoru Rm
metriku:
(x, y) = max
1≤i ≤m
|xi − yi |.
Prostor Rm
s takto definovanou metrikou je úplným
metrickým prostorem. Nyní lze pro vyšetřování konvergence
iteračního procesu xk+1
= G(xk
), k = 0, 1, 2, . . . užít
Banachovy věty o pevném bodě.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 20 / 24
Věta
Nechť zobrazení G : Rm
→ Rm
je kontrakce na Rm
,
(G(x), G(y)) ≤ q (x, y), 0 ≤ q < 1.
Pak pro každou počáteční aproximaci x0
∈ Rm
je posloupnost
xk ∞
k=0
, k = 0, 1, 2, . . ., xk
= G(xk−1
), konvergentní v Rm
a lim
k→∞
xk
= ξ, kde ξ je jediný pevný bod zobrazení G.
Věta
Nechť ξ ∈ Rm
je pevný bod rovnice x = G(x). Nechť funkce
gi , i = 1, . . . , m, mají spojité parciální derivace pro všechna
x ∈ Ω(ξ, r), Ω(ξ, r) = {x| (x, ξ) ≤ r}. Nechť dále platí
∂gi (x)
∂xj
≤
q
m
, i, j = 1, . . . , m,
0 ≤ q < 1 a nechť x0
∈ Ω(ξ, r). Pak všechny iterace xk ∞
k=0
určené vztahem xk+1
= G(xk
) leží v množině Ω(ξ, r)
a lim
k→∞
xk
= ξ.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 21 / 24
Jednodušší předpoklad:
max
1≤i ≤m
m
j=1
∂gi (x)
∂xj
≤ q < 1.
Seidelova matoda
Pro výpočet xk+1
i použijeme již vypočtených hodnot xk+1
1 , . . . ,
xk+1
i−1 , tj.
xk+1
1 = g1(xk
1 , xk
2 , . . . , xk
m)
xk+1
2 = g2(xk+1
1 , xk
2 , . . . , xk
m)
xk+1
3 = g3(xk+1
1 , xk+1
2 , . . . , xk
m)
...
xk+1
m = gm(xk+1
1 , . . . , xk+1
m−1, xk
m).
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 22 / 24
Newtonova metoda pro systémy nelineárních rovnic
F(x) = o, F ∈ C2
(O(ξ))
JF (x) =






∂f1(x)
∂x1
· · ·
∂f1(x)
∂xm
...
...
...
∂fm(x)
∂x1
· · ·
∂fm(x)
∂xm






Taylorův rozvoj:
F(x + h) = F(x) + JF (x) · h + O( h 2
) · (1, . . . , 1)T
x = xk
, xk+1
= xk
+ h ⇒ h = xk+1
− xk
Zanedbáme chybový člen,
F(x + h) = o ⇒ JF (xk
)(xk+1
− xk
) = −F(xk
)
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 23 / 24
xk+1
= xk
− J−1
F (xk
)F(xk
)
Iterační funkce
G(x) = x − J−1
F (x)F(x)
Věta
Nechť ξ je kořenem rovnice F(x) = o. Nechť JF (x) je
regulární matice se spojitými prvky v okolí O(ξ) bodu ξ,
přičemž
J−1
F (x) ∞
≤ K, K = konst.,
pro všechna x z tohoto okolí. Nechť funkce fi , i = 1, . . . , m,
mají spojité druhé parciální derivace v O(ξ).
Posloupnost xk ∞
k=0
určená Newtonovou metodou konverguje
ke kořenu ξ za předpokladu, že počáteční aproximace x0
leží
dostatečně blízko ξ. Řád metody je roven dvěma.
Jiří Zelinka Numerické metody 9. přednáška, 19. dubna 2018 24 / 24