Pokročilé numerické metody I Systémy lineárních rovnic Jiří Zelinka Jiří Zelinka Systémy lineárních rovnic 1 / 40 Iterační metody pro systémy lineárních rovnic Systém Ax = b převedeme na x = Tx + g ˆx – řešení ˆx = (E − T)−1 g za předpokladu, že E − T je regulární. x0 ∈ Rn – libovolná počáteční aproximace. Posloupnost xk ∞ k=0 určená rekurentně vztahem xk+1 = Txk + g, k = 0, 1, . . . se nazývá iterační posloupnost a matice T se nazývá iterační matice. Jiří Zelinka Systémy lineárních rovnic 2 / 40 Problémy: 1 Jak zvolit iterační matici T, tj. jakým způsobem převést systém Ax = b na systém x = Tx + g? 2 Za jakých předpokladů posloupnost xk ∞ k=0 konverguje pro libovolnou počáteční aproximaci k přesnému řešení ˆx? x1 = Tx0 + g, x2 = Tx1 + g = T(Tx0 + g) + g = T2 x0 + (T + E)g, x3 = Tx2 + g = T3 x0 + (T2 + T + E)g, ... xk+1 = Tk+1 x0 + (Tk + Tk−1 + . . . + E)g. Definice Řekneme, že matice H je konvergentní, jestliže lim k→∞ Hk = O, kde O je nulová matice, konvergence je bodová. Jiří Zelinka Systémy lineárních rovnic 3 / 40 Věta 1 Následující tvrzení jsou ekvivalentní: 1 H je konvergentní matice. 2 lim k→∞ ∥Hk ∥ = 0 pro nějakou přidruženou maticovou normu. 3 ρ(H) < 1 (ρ(H) je spektrální poloměr H). 4 lim k→∞ Hk x = o pro libovolný vektor x ∈ Rn . Lemma Nechť ρ(T) < 1. Pak E − T je regulární a platí (E − T)−1 = E + T + T2 + . . . Jiří Zelinka Systémy lineárních rovnic 4 / 40 Věta 2 (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 Důsledek Nechť pro nějakou přidruženou maticovou normu platí ∥T∥ < 1. Pak posloupnost xk ∞ k=0 generovaná iteračním procesem x = Tx + g konverguje k řešení ˆx = (E − T)−1 g pro každou počáteční aproximaci x0 ∈ Rn . Dále platí ∥ˆx − xk ∥ ≤ ∥T∥k ∥ˆx − x0 ∥, ∥ˆx − xk ∥ ≤ ∥T∥k 1 − ∥T∥ ∥x1 − x0 ∥. Jiří Zelinka Systémy lineárních rovnic 5 / 40 Jacobiova iterační metoda Ax = b, A =      a11 · · · · · · a1n a21 a2n ... ... an1 · · · · · · ann      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 ), Jiří Zelinka Systémy lineárních rovnic 6 / 40 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. Na pravé straně jsou vždy všechny prvky Ax kromě diagonálního, kterým se rovnice vydělí. Jiří Zelinka Systémy lineárních rovnic 7 / 40 Maticový zápis: Matici A zapišme ve tvaru A = D + L + U, kde D =    a11 0 ... 0 ann    , Jiří Zelinka Systémy lineárních rovnic 8 / 40 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 Systémy lineárních rovnic 9 / 40 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 Systémy lineárních rovnic 10 / 40 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 Systémy lineárních rovnic 11 / 40 Věta 3 (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 Systémy lineárních rovnic 12 / 40 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 Systémy lineárních rovnic 13 / 40 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 Systémy lineárních rovnic 14 / 40 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 Systémy lineárních rovnic 15 / 40 Věta 4 Posloupnost xk ∞ k=0 generovaná Gaussovou-Seidelovou iterační metodou xk+1 = −(D + L)−1 Uxk + (D + L)−1 b. 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 Systémy lineárních rovnic 16 / 40 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 5 Nechť A je pozitivně definitní matice. Pak Gaussova-Seidelova metoda konverguje pro každou počáteční aproximaci. Jiří Zelinka Systémy lineárních rovnic 17 / 40 Věta 6 (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. Jiří Zelinka Systémy lineárních rovnic 18 / 40 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 Systémy lineárních rovnic 19 / 40 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 7 (Kahan) Nechť aii ̸= 0, i = 1, . . . , n. Pak ϱ(Tω) ≥ |ω − 1|. Důsledek: Má smysl uvažovat jen ω ∈ (0, 2). Jiří Zelinka Systémy lineárních rovnic 20 / 40 Věta 8 (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 9 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 Systémy lineárních rovnic 21 / 40 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 Systémy lineárních rovnic 22 / 40 Metody založené na minimalizaci kvadratické formy Věta 10 Jestliže A je pozitivně definitní matice, pak řešení systému Ax = b je ekvivalentní minimalizaci kvadratické funkce Q(x) = 1 2 xT Ax − xT b. Tato kvadratická funkce má jediné minimum, kterého nabývá v řešení systému Ax = b, tj. pro ˆx = A−1 b. Důsledek min Q(x) = Q(ˆx) = −1 2 bT A−1 b. Jiří Zelinka Systémy lineárních rovnic 23 / 40 Geometrická interpretace Nechť ˆx je řešení systému Ax = b a nechť x = ˆx + s. Pak Q(x)+ 1 2 ˆxT b = 1 2 (ˆx+s)T A (ˆx+s)−(ˆx+s)T b+ 1 2 ˆxT b = 1 2 sT As . Plocha o rovnici 1 2 sT As = konst. je hyperelipsoid v proměnných s1,. . . ,sn se středem v s = o, tj. v ˆx, tedy i rovnice Q(x) = konst. představuje hyperelipsoid, směry jeho os jsou dány vlastními vektory matice A, délky poloos vlastními čísly. Minimalizace funkce Q: vybereme počáteční aproximaci x1 a pak určíme x2 tak, že zvolíme nějaký směr v1 a vzdálenost t1 ve směru v1. Obecně xk+1 = xk + tkvk, i = 1, 2, . . . Vektory vk se nazývají směrové vektory. Jiří Zelinka Systémy lineárních rovnic 24 / 40 Metoda maximálního spádu (steepest descent) x1 je počáteční aproximace řešení ˆx, položme r1 = b − Ax1: reziduový vektor. Hledáme takový vektor v1, ∥v1∥ = ∥b − Ax1∥, pro který δ δt Q(x1 + tv1) t=0 = min. Výsledek: v1 = r1 = b − Ax1. Dále hledáme minimum Q ve směru v1 = r1 od bodu x1, tedy minimum kvadratické funkce Q(x1 + tr1) v bodě t1. Výsledek: t1 = rT 1 r1 rT 1 Ar1 , takže x2 = x1 + rT 1 r1 rT 1 Ar1 r1 a obecně tk+1 = rk T rk rk T Ark , xk+1 = xk + rk T rk rk T Ark rk . Jiří Zelinka Systémy lineárních rovnic 25 / 40 Vlastnosti iteračního procesu: 1 rk+1 = rk − rk T rk rk T Ark Ark = rk − tkArk 2 (rk+1, rk) = rT k rk+1 = 0 3 Q(xk+1) < Q(xk), k = 1, 2, . . . Rychlost konvergence Buď ˆx řešení systému Ax = b. Pak přímým výpočtem dostaneme Q(ˆx) = −1 2 bT ˆx = −1 2 bT A−1 b. Věta 11 Označme E(xk) = Q(xk) − Q(ˆx) = 1 2 (xk − ˆx)T A(xk − ˆx). Pak E(xk+1) ≤ λmax −λmin λmax +λmin 2 E(xk) Jiří Zelinka Systémy lineárních rovnic 26 / 40 -10 -5 0 5 10 -10 -5 0 5 10 Metoda maxmálního spádu Jiří Zelinka Systémy lineárních rovnic 27 / 40 Metoda sdružených (konjugovaných) gradientů Nechť p1, . . . , pn tvoří bázi, ˆx je řešení Ax = b, ˆx = n i=1 ci pi . Pokud jsou pi ortogonální, pak ci = ( ˆx,pi ) (pi ,pi ) . Problém: neznáme ˆx. Výchozí rovnost násobíme maticí A: Aˆx = b = n i=1 ci Api . Kdyby platilo (Apj , pi ) = pT i Apj = 0 pro i ̸= j, pak ci = (b, pi ) (Api , pi ) . Proto definujme nový skalární součin (x, y)A = (Ax, y), vektory, pro něž (x, y)A = 0, se nazývají A-ortogonální (A-sdružené, A-konjugované). Cílem metody je sestrojit A-ortogonální vektory p1, . . . , pn podobně jako u Gramova–Schmidtova ortogonalizačního procesu. Jiří Zelinka Systémy lineárních rovnic 28 / 40 Modifikovaný Gramův–Schmidtův ortogonalizační proces: Z dané báze r1, . . . , rn vytváříme novou, A-ortogonální bázi. Položme p1 = r1 a dále pk = rk + k−1 i=1 βk,i pk, pak βk,i = − (rk, pi )A (pi , pi )A . Jako bázi rk bereme postupně rezidua, v každém kroku počítáme koeficienty α i β. Snadno se ověří (viz metoda maximálního spádu), že pro iteraci xk a směr daný vektorem p má funkce Q minimum v xk+1 = xk + αkp, pro αk = pT rk pT Ap = (rk, p) (p, p)A . Jiří Zelinka Systémy lineárních rovnic 29 / 40 Počáteční krok x1 – počáteční iterace, p1 = r1 = b − Ax1. Následně x2 = x1 + α1p1 pro α1 = (r1,p1) (p1,p1)A , potom r2 = b − Ax2 = r1 − α1Ap1 a p2 = r2 + β2,1p1 pro β2,1 = − (r2,p1)A (p1,p1)A . Obecný krok Pro iteraci xk, reziduum rk = b − Axk a směrový vektor pk položíme xk+1 = xk + αkpk, kde αk = (rk ,pk ) (pk ,pk )A . rk+1 = b − Axk+1 = rk − αkApk, pk+1 = rk+1 + k i=1 βk+1,i pk pro βk+1,i = − (rk+1,pi )A (pi ,pi )A . Jiří Zelinka Systémy lineárních rovnic 30 / 40 Věta 12 (Vlastnosti definovaných vektorů) 1. L(r1, . . . , rk) = L(p1, . . . , pk) pro k = 1, 2, . . . . 2. Vektory p1, p2, . . . jsou A-ortogonální. 3. Vektor rk − rk+1 je kolmý k pi pro i ̸= k. Důkaz: rk − rk+1 = αkApk, takže pro i ̸= k dostáváme (pi , rk − rk+1) = αk(pi , Apk) = αk(pi , pk)A = 0. 4. (rk, pk)A = (pk, pk)A, ∀k. Důkaz indukcí: Pro k = 1 platí, nechť platí pro k, pak (rk+1, pk+1)A = (pk+1− k i=1 βk+1,i pi , pk+1)A = (pk+1, pk+1)A Jiří Zelinka Systémy lineárních rovnic 31 / 40 5. Vektory r1, r2, . . . jsou ortogonální. Důkaz indukcí: Nejprve dokážeme (r2, r1) = 0: (r2, r1) = (r1, r1)−α1(Ap1, r1) = (r1, r1)− (r1, p1) (p1, p1)A (p1, r1)A = 0, neboť p1 = r1. Dále nechť jsou kolmé vektory r1, . . . , rk pro k ≥ 2. Pak (rk+1, rk) = (rk, rk) − αk(Apk, rk) = (rk, rk)− − (rk ,pk ) (pk ,pk )A (pk, rk)A(viz 4.) = (rk, rk) − (rk, rk + i