Numerické metody 9. přednáška, 21. dubna 2016 Jiří Zelinka Jiří Zelinka Numerické metody 9. prednáška, 21. dubna 2016 1/24 Opakování Iteracní metody resení systémů lineárních rovnic Systém Ax = b převedeme na x = Tx + g x k+l = TxK + g, k = 0,1,... Hlavní veta o konvergenci iteracního procesu Posloupnost {xk}™ určená iteračním procesem x = Tx + g konverguje pro každou počáteční aproximaci x° g W právě tehdy, když p(T) < 1, přičemž lim xk = x*, x* = Tx* + g k^řoc Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 2/24 Jacobiova iterační metoda Ax = b, a = / 3n 321 Matici A zapišme ve tvaru A = D-L-U 3\„ \ 32 n 'nn kde / au D = \ 0 0 \ ' nn J Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 3/24 L = U = í o — 321 ■ V -a„i ■ ( O -312 V o a„,„-i O / -ain \ 3/7—1,n O D je diagonálni 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. Jirí Zelinka Numerické metody 9. přednáška, 21. dubna 2016 4/24 Ax = {D-L -U)x = b Dx = (L+U)x + b. Pokud a,-,- /O, / = 1,..., n, je matice D regulární a z předchozí rovnice lze vypočítat x = D-X(L+ U)x + D~lb. -i /jit D-1 = \ 0 o \ i. / Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 5/24 Maticový tvar Jacobiovy iteracní metody Jacobiova iterační matice: Tj — D~1(L+ U) x k+l _ = Tjx* + D-1 b, Tj = (tij), tjj = -0 pro i ^ j, ta = 0 pro i = 1,..., n. í Tj = \ 2nl ]nn 2n2 ]nn 0 ai2 31" ^ 321 322 0 ^2/7 322 0 0"^ = 322 Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 6/24 Realizace výpočtu: Z první rovnice vypočteme xi: 1 anxi+3i2X2H-----\-alnxn = 6X x1 = —(61-312*2-•. .-3i„xn) Xf+1 = —(bí - 312X2* - ... - 3i„x£), z druhé rovnice vypočteme x2: *2 = -(P2 - 321*1 ~ a23*3 ~ • • • ~ 3l#i*J5 322 obecně z /-té rovnice vypočteme x,\ . . // // 7=1 až z /7-té rovnice vypočteme xn, a na pravé straně takto získaného systému jsou prvky matice T/. < n > < fl > < , > < , > == Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 7/24 Veta o konvergenci Jacobiovy iteracní metody: Posloupnost {xk}™ generovaná metodou x^1 = Tjxk + D~lb konverguje pro každou počáteční aproximaci x° g Rn právě tehdy, když g(Tj) < 1. Odhad chyby: x — x OO — T k co 1 - T OO x — X o oo Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 8/24 Silné řádkové sumační kriterium: Nechť matice A je ryze řádkově diagonálně dominantní, tj n au > \au i = 1.....n. 7=1 Pak Jacobiova iterační metoda konverguje pro každou počáteční aproximaci x° g IR". Silné sloupcové sumační kriterium: Nechť matice A je ryze sloupcově diagonálně dominantní, tj n 2kk > k — 1..... n. i //c Pak Jacobiova iterační metoda konverguje pro každou počáteční aproximaci x° e IR". 4 n „ 1= s} <\(y Jiří Zelinka Numerické metody 9. prednáška, 21. dubna 2016 9/24 Gaussova-Seidelova iterační metoda Z první rovnice vypočteme xi: aiix1+a12x2-\-----^alnxn = bľ xx = —{bi-a12x2-.. .-a^x,,) x*+1 — —(61 — ai2X^ — ... — a^x^), z druhé rovnice vypočteme x2, pro xi použijme novou iteraci x k+l 2 — -(^2 321*1 ^3X3 — • • • — a22 ze třetí rovnice vypočteme x3, pro xi a x2 použijme novou iteraci: *3+1 — —(63 — a3ix^+1 — a32X2+1 — 334X4 ... — a^x^), ^33 Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 10 / 24 Obecně /-i x"+1 = - Z_j o.. J j=i " £« xfc + El . .... "II dll 7=1+1 / = 1,..., n. Maticový zápis: Ax = b (D-L-U)x = b {D-L)x = Í7x + b. a,v 0, / = 1,..., /?, ^> matice D — /_ je regulární a x = (D - L)"1^ + (D - L)"1^ Položme Tg = (D — L)~lU, Gaussova-Seidelova iterační metoda je tvaru -i x k+l = TGxk + g, g = (D-L)~1b. Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 11 /24 Veta Posloupnost {xk}™ generovaná Gaussovou-Seidelovou iterační metodou xk+l = (D - Ľ)~lUxk + (D - L)'1 b. konverguje pro každou počáteční aproximaci x° g W právě tehdy, když g(Tc) < 1. Silné řádkové sumační kriterium: Nechť matice A je ryze řádkově diagonálně dominantní, tj. n '// > 1 = 1 • • • /7. Pak Gaussova-Seidelova iterační metoda konverguje pro každou počáteční aproximaci x° e IR". Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 12/24 Silné sloupcové sumační kriterium: Nechť matice A je ryze sloupcově diagonálně dominantní, tj n 2kk > k — 1..... n. i //c Pak Gaussova-Seidelova iterační metoda konverguje pro každou počáteční aproximaci x° e IR". Příklad Geometrický význam Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 13/24 Veta (Stein-Rosenberg) Nechť pro prvky matice A platí a,y < 0 pro všechna / ^ j a a,-,- > 0, / = 1,..., n. Pak platí právě jedno z následujících tvrzení: • 0 < q(TG) < q(Tj) < 1 • 1 < q(Tj) < g{TG) • q(Tj) = q(Tg) = 0 o g(Tj) = q(Tg) = 1. To znamená, že konvergují-li obě metody, Gaussova-Seidelova metoda konverguje rychleji. 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, 21. dubna 2016 14/24 Relaxační metody Modifikace Gaussovy-Seidelovy metody, oj - relaxační parametr xi k+l _ = (1 - u)xf + oj au bi- i-i 7=1 n k+1 j=i+l SijXj Relaxační metodu lze maticově zapsat takto x^1 = (D - - u)D + ojU]xk + oj{D - ojL)~1b TU = (D- ojL)-1^! - lo)D + loU] Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 15/24 Hodnoty parametru uj\ Pro 0 < uj < 1 se iterační metody nazývají metodami dolní relaxace. Tyto metody jsou vhodné v případě, že Gaussova-Seidelova metoda nekonverguje. Pro uj = 1 je relaxační metoda totožná s Gaussovou- Seidelovou metodou. Pro 1 < uj se metody nazývají metodami horní relaxace, nebo častěji SOR metodami (SOR = Successive Over-Relaxation). Tyto metody lze užít ke zrychlení konvergence Gaussovy-Seidelovy metody. Věta (Kahan). Nechť a,-,- 0, / = 1,..., n. Pak q{Tu) > \cj - 1 . Důsledek: Má smysl uvažovat jen uj g (0,2). Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 16/24 Veta (Ostrowski-Reich). Pro pozitivně definitní matici A platí g(Tu) < 1 pro všechna uj e (0,2). Třídiagonální matice: A = {a;j);J=n aij = 0, pro |/ -j\ > 1 Veta. Nechť A je třídiagonální pozitivně definitní matice. Pak q{Tg) = q2(Tj) < 1 a optimální hodnota relaxačního parametru je dána vztahem ^ — ^opt — l + y/l-r(Tj)' Při této volbě je q{Tu) = |1 — uj Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 17/24 Cykly v iteračních metodách Například pro systémy xi + kx2 = bi xi — kx2 = b2 Jacobiova metoda: cyklus délky 4 Gaussova-Seidelova metoda: cyklus délky 2 Relaxační metody: cykly různých délek Xi + x2 = 1 qxi + x2 = 1. Pro u = 2: -1 -2 2g 4q-l To = Ai,2 = cos (^±/sin q = -(1+cos cp) cp = 2ttI/p, 0 < / < p/2 ^> existuje cyklus délky p. Body cyklu leží na elipse se středem v hledanéiu řešení. ^ Jiří Zelinka 1= ^) c\ o Numerické metody 9. přednáška, 21. dubna 2016 18/24 Iterační metody pro systémy nelineárních rovnic fi{xi,... ,xm) = o Kořen systému: uspořádaná /77-tice reálných číše C = (£l •)••••) £>m)i která tomuto systému vyhovuje. Vektorový tvar: F(x) = 0, xGr,o = (Or..,0)rg Rm. Jiří Zelinka Numerické metody 9. přednáška, 21. dubna 2016 19/24 Systém převedeme na ekvivalentní rovnici x = G(x), xgR m neboli *i = • • • ,xm) Xm = gm{xU...,Xm) a budeme hledat pevný bod zobrazení G : Rm —R Definujme nyní v prostoru K"7 metriku: g(x,y) = max |x; - y; Prostor IRm s takto definovanou metrikou je úplným metrickým prostorem. Nyní lze pro vyšetřování konvergence iteračního procesu xk+l = G(x*), k = 0,1, 2,... užít Banachovy věty o pevném bodě. < n > < _ „ < = , < =, Jiří Zelinka Numerické metody 9. prednáška, 21. dubna 2016 20/24 Veta Nechť zobrazení G : rm rm je kontrakce na rm, q(G(x), G(y)) < qg{x,y), 0 < q < 1. Pak pro každou počáteční aproximaci x°