Numerické metody
9. přednáška, 26. dubna 2019
Jiří Zelinka
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 1 / 38
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, 26. dubna 2019 2 / 38
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, 26. dubna 2019 3 / 38
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, 26. dubna 2019 4 / 38
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, 26. dubna 2019 5 / 38
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, 26. dubna 2019 6 / 38
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, 26. dubna 2019 7 / 38
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, 26. dubna 2019 8 / 38
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, 26. dubna 2019 9 / 38
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, 26. dubna 2019 10 / 38
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, 26. dubna 2019 11 / 38
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, 26. dubna 2019 12 / 38
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, 26. dubna 2019 13 / 38
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, 26. dubna 2019 14 / 38
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, 26. dubna 2019 15 / 38
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, 26. dubna 2019 16 / 38
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, 26. dubna 2019 17 / 38
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, 26. dubna 2019 18 / 38
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, 26. dubna 2019 19 / 38
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, 26. dubna 2019 20 / 38
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, 26. dubna 2019 21 / 38
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, 26. dubna 2019 22 / 38
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, 26. dubna 2019 23 / 38
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, 26. dubna 2019 24 / 38
Řešení systémů lin. rovnic – přímé metody
Základní pojmy
Ax = b, A =



a11 · · · a1n
...
...
an1 · · · ann


, x =



x1
...
xn


, b =



b1
...
bn


.
A – matice soustavy, regulární. Řešení: ˜x = A−1
b
Rozšířená matice soustavy:
(A | b) =





a11 · · · a1n b1
a21 · · · a2n b2
...
...
...
an1 · · · ann bn





Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 25 / 38
A – dolní trojúhelníková: aij = 0 pro i < j.
A – horní trojúhelníková: aij = 0 pro i > j.
A – pásová, jestliže existují p, q, 1 < p, q < n taková, že
aij = 0, jestliže i + p ≤ j nebo j + q ≤ i, šířka pásu
w = p + q − 1.
A – třídiagonální pro p = q = 2.
A – ryze řádkově diagonálně dominantní
|aii | >
n
j=1
j=i
|aij |, i = 1, . . . , n.
A – ryze sloupcově diagonálně dominantní – podobně
Věta:
Ryze řádkově (sloupcově) diagonálně dominantní matice, je
regulární.
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 26 / 38
Gaussova eliminační metoda
Úprava soustavy na soustavu s horní trojúhelníkovou maticí R:
(A | b) −→ R | ˜b
Pak provádíme tzv. zpětný chod – počítáme řešení od poslední
složky k první.
Elementární úpravy a matice úprav
Nemění řešení, každá úprava má inverzi.
1. násobení řádku nenulovou konstantou c
Ic =






1 0 · · · · · · · · · 0
0 1 0 · · · · · · 0
.
.
.
.
.
.
.
.
.
.
.
. c
.
.
.
..
.
.
.
.
0 0 · · · · · · · · · 1






, I−1
c =






1 0 · · · · · · · · · 0
0 1 0 · · · · · · 0
.
.
.
.
.
.
.
.
.
.
.
. 1/c
.
.
.
..
.
.
.
.
0 0 · · · · · · · · · 1






Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 27 / 38
2. výměna řádků i, k
Pi,k =






















1 0 · · · 0
0 1 0
...
1 0 0 · · · 0 0 0 · · · 0
0 0 0 · · · 0 1 0 · · · 0
0 0 1 0 0 · · · 0
...
...
...
...
...
...
...
0 0 1 0 0 · · · 0
0 1 0 · · · 0 0 0 · · · 0
0 0 0 · · · 0 0 1
...
...
0 0 1






















, P−1
i,k = Pi,k
Pi,k – permutační matice
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 28 / 38
3. přičtení c násobku i-tého řádku ke k-tému, i < k
Gi,k,c =
i
k














1 0 · · · 0
0 1 · · · 0
...
...
1
...
c 1
...
0 1














.
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 29 / 38
G−1
i,k,c =
i
k














1 0 · · · 0
0 1 · · · 0
...
...
1
...
−c 1
...
0 1














.
Ve skutečnosti pro převod na trojúhelníkovou matici stačí 2. a
3. elementární úprava.
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 30 / 38
Postup při Gaussově eliminaci
(1) výměna 1. a k-tého řádku (v případě potřeby)
(A | b) −→ A(1)
| b(1)
, A(1)
| b(1)
= P1,k · (A | b)
(1’) vynulování prvního sloupce pod hlavní diagonálou
A(1 )
| b(1 )
= G1· A(1)
| b(1)
, G1 =





1 · · · · · · 0
l21 1
...
...
...
ln1 0 · · · 1





,
lk1 = −
a
(1)
k1
a
(1)
11
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 31 / 38
(i) výměna i-tého. a k-tého řádku (v případě potřeby)
A(i)
| b(i)
= Pi,k · A(i−1 )
| b(i−1 )
(i’) vynulování i-tého sloupce pod hlavní diagonálou
A(i )
| b(i )
= Gi · A(i)
| b(i)
,
Gi =











1 · · · 0
...
... 1
li+1,i
...
...
...
0 ln,i 1











, lki = −
a
(i)
ki
a
(i)
ii
i = 2, . . . , n − 1
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 32 / 38
Gaussova eliminace bez výměny řádků:
R | ˜b = Gn−1 · . . . G2 · G1 · (A | b)
tedy
R = Gn−1 · . . . G2 · G1 · A
odkud
G−1
1 G−1
2 . . . G−1
n−1R = A.
Matice Gi jsou dolní trojúhelníková, tedy G−i
1 jsou také dolní
trojúhelníkové, takže
A = L · R, L = G−1
1 G−1
2 . . . G−1
n−1.
L – dolní trojúhelníková matice.
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 33 / 38
G−1
i =











1 · · · 0
...
... 1
−li+1,i
...
...
...
0 −lni 1











,
pak
L = G−1
1 G−1
2 . . . G−1
n−1 =





1 · · · · · · 0
−l21
...
...
...
...
...
−ln1 · · · −ln,n−1 1





.
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 34 / 38
LR rozklad
A = L · R: LR (též LU) rozklad matice A
Použití při řešení soustavy: substituce Rx = y, řešíme
soustavu Ly = b s dolní trojúhelníkovou maticí, pak Rx = y s
horní trojúhelníkovou maticí.
LR rozklad s výměnou řádků
Provádíme LR rozklad matice P · A, kde P je vhodná
permutační matice, která provádí výměnu řádků.
Při praktickém výpočtu, pokud narazíme na potřebu vyměnit
řádky, postupujeme takto:
Pokud bychom vyměnili řádky předem, v už vypočítané části
matice L by tyto řádky byly vyměněny, zbytek by se nezměnil.
Proto můžeme vyměnit řádky ve vypočítané části matice L.
Dále v pomocném vektoru p na začátku nastaveném na
(1, 2, . . . , n)T
zaznamenáme výměnu řádků, tedy vyměníme v
něm stejné řádky.
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 35 / 38
Příklad
2x1 + 4x2 − x3 = −5
x1 + x2 − 3x3 = −9
4x1 + x2 + 2x3 = 9
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 36 / 38
Výběr vedoucího prvku (pivota) – částečný
Při úpravě i-tého sloupce najdeme mezi prvky a
(i−1)
ii , . . . , a
(i−1)
ni
prvek s maximální absoulutní hodnotou (např. a
(i−1)
ki ), pak
vyměníme i-tý a k-tý řádek.
Výběr vedoucího prvku (pivota) – úplný
Hledáme prvek s maximální absoultní hodnotou mezi a
(i−1)
jk ,
i ≤ j ≤ n, i ≤ k ≤ n, pak vyměníme příslušný řádek a
spoupec, čímž se změní pořadí proměnných.
Příklad:
LR rozklad matice
A =
0, 0001 1
1 1
.
se zaokrouhlováním na 3 platné číslice.
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 37 / 38
Věta
Nechť všechny hlavní minory matice A ∈ Mn jsou různé od
nuly, tj.
a11 = 0,
a11 a12
a21 a22
= 0,
a11 a12 a13
a21 a22 a23
a31 a32 a33
= 0, . . . det A = 0.
Pak matici A lze rozložit na součin dolní a horní trojúhelníkové
matice.
Důsledky
Nechť matice A je pozitivně definitní. Pak GEM lze
provést bez výměny řádků a sloupců.
Nechť matice A je ryze řádkově diagonálně dominantní.
Pak GEM lze provést bez výměny řádků a sloupců.
Jiří Zelinka Numerické metody 9. přednáška, 26. dubna 2019 38 / 38