3. Iterační metody řešení soustav Tento učební text byl podpořen z Operačního programu Praha - Adaptabilita
Irena Sýkorová
Soustavy lineárních rovnic Uvažujeme soustavu lineárních rovnic
A¯ x = ¯b,
(1)
kde matice soustavy A = (aij ) je reálná regulární matice řádu n, ¯b = (b1 , b2 , . . . , bn )′ je sloupcový (∗) (∗) (∗) vektor pravých stran, hledané řešení označíme x ¯(∗) = (x1 , x2 , . . . , xn )′ . Soustavu (1) vyjádříme ve tvaru x ¯ = Mx ¯ + v¯, (2) kde M je vhodná matice, v¯ je sloupcový vektor. Iterační metoda spočívá v tom, že zvolíme nějaké počáteční přiblížení x ¯(0) a podle vzorce x ¯(k+1) = M x ¯(k) + v¯,
k = 0, 1, . . .
(3)
počítáme posloupnost vektorů x ¯(k) , která konverguje k řešení původní soustavy x ¯(∗) . Výpočet ukončíme, až k¯ x(k+1) − x ¯(k) k < ε nebo po daném počtu iterací. Podmínky konvergence: a) nutná a postačující Metoda (3) konverguje pro každé počáteční přiblížení x ¯(0) , právě když spektrální poloměr iterační matice ρ(M ) < 1. b) postačující Jestliže pro nějakou normu matice M platí kM k < 1, pak metoda (3) konverguje pro každé počáteční přiblížení x ¯(0) . Je-li splněná postačující podmínka, můžeme odhadnout velikost chyby v k-tém kroku: k¯ x∗ − x ¯(k) k ≤
kM k k¯ x(k) − x ¯(k−1) k 1 − kM k
nebo k¯ x∗ − x ¯(k) k ≤
kM kk k¯ x(1) − x ¯(0) k. 1 − kM k
Popíšeme dvě základní iterační metody. K tomu bude potřeba vyjádřit si matici soustavy A ve tvaru A = D + L + U, (4) kde D je diagonální matice, L je ostře (s nulami na diagonále) dolní trojúhelníková a U je ostře horní trojúhelníková matice, tj.
a11 a21 . ..
an1
a12 a22 .. .
... ... .. .
an2
...
a1n a2n = .. .
ann
a11 0 = ... 0
2
0 a22 .. .
... ... .. .
0
...
0 0 0 a21 + . .. . .. an1 ann
0 0 .. .
... ... .. .
an2
...
0 0 0 0 +. .. . ..
0
0
a12 0 .. .
... ... .. .
0
...
a1n a2n , .. . 0
Jacobiova metoda Matici soutavy A vyjádříme ve tvaru Pak soustavu (1) můžeme upravit A¯ x = ¯b (D + L + U ) x ¯ = ¯b D¯ x = −(L + U ) x ¯ + ¯b x ¯ = −D−1 (L + U ) x ¯ + D−1¯b, tedy ve vzorci (2) volíme MJ = −D−1 (L + U ),
v¯J = D−1¯b.
Zvolíme tedy x ¯(0) a počítáme pro k = 0, 1, 2, . . . x ¯(k+1) = −D−1 (L + U ) x ¯(k) + D−1¯b
(J)
a rozepsáním po souřadnicích dostáváme předpis vhodný pro počítání (k+1)
xi
=−
1 aii
n X
(k)
aij xj
+
j=1, j6=i
bi , aii
i = 1, . . . , n,
k = 0, 1, . . .
(4)
Podmínky konvergence: a) ρ(MJ ) < 1 ⇔ Jacobiova metoda konverguje. Vlastní čísla matice MJ lze získat řešením rovnice det(L + λD + U ) = 0. b) Existuje norma, pro kterou kMJ k < 1 ⇒ Jacobiova metoda konverguje. c) A je ostře diagonálně dominantní ⇒ Jacobiova metoda konverguje. Příklad 1: Je dána soustava lineárních rovnic 4x1 + x2 − 2x3 = 4 2x1 + 5x2 + 2x3 = 5
(5)
x1 + x2 − 5x3 = 4. Dokážeme, že Jacobiova metoda bude konvergovat, zvolíme počáteční přiblížení x ¯(0) = (0, 0, 0)′ a vypočítáme první dvě iterace. Řešení: Matice soustavy
4 A = 2 1
1 −2 5 2 1 −5
je ostře diagonálně dominantní, např. pro řádky platí |4| |5| | − 5|
> | 1 | + | − 2| > |2| + |2| > | 1 | + | 1 |.
Je tedy splněná postačující podmínka konvergence c), Jacobiova metoda konverguje.
3
Soustavu (5) upravíme tak, že z první rovnice vyjádříme x1 , ze druhé x2 a ze třetí x3 . x1 = − 14 ( x2 = x3 =
x2 − 2x3 − 4)
− 51 (2x1 1 5 ( x1 +
+ 2x3 − 5) − 4).
x2
Vezmeme vektor x ¯(0) a podle vzorce (4) počítáme jednotlivé souřadnice vektoru x ¯(1) : (0)
(1)
(0)
x1 = − 41 (x2 − 2x3 − 4) = − 14 (0 − 0 − 4) = 1 (1)
(0)
(0)
x2 = − 15 (2x1 + 2x3 − 5) = − 15 (0 + 0 − 5) = 1 (1)
x3 =
1 5(
(0)
(0)
x1 + x2 − 4) =
1 5 (0
+ 0 − 4) = − 54 .
První iterace je tedy x ¯(1) = (1, 1, − 45 )′ = (1; 1; −0, 8)′ . Pomocí tohoto vektoru počítáme druhou iteraci: (2)
(1)
(1)
x1 = − 14 (x2 − 2x3 − 4) = − 14 (1 + 2 · (1)
(2)
(1)
4 5
x2 = − 51 (2x1 + 2x3 − 5) = − 51 (2 − 2 · (2)
x3 =
1 5(
(1)
(1)
x1 + x2 − 4) =
1 5 (1
− 4) = 4 5
7 20
− 5) =
23 25
+ 1 − 4) = − 52 .
7 23 , 25 , − 25 )′ = (0, 35; 0, 92; −0, 4)′ . Druhá iterace je x ¯(2) = ( 20
Pro porovnání uvedeme přesné řešení soustavy x ¯∗ = ( 12 , 1, − 12 )′ = (0, 5; 1; −0, 5)′ a provedeme odhad chyby, např. pro normu k · k∞ k¯ x∗ − x ¯(2) k∞ ≤
kMJ k∞ 0, 8 k¯ x(k) − x ¯(k−1) k∞ = · 0, 65 = 2, 6. 1 − kMJ k∞ 1 − 0, 8
Odhad je v tomto případě velmi hrubý, protože k¯ x∗ − x ¯(2) k∞ = 0, 15. MATLAB – Jacobiova metoda A – matice soustavy, b – vektor pravých stran, x(0) – počáteční přiblížení výpočet matice MJ >> U=triu(A,1) >> L=tril(A,-1) >> D=A-U-L >> MJ=(-1)*inv(D)*(U+L) výpočet vektoru vJ >> vJ=inv(D)*b výpočet jednotlivých iterací x(1) , x(2) , atd. >> x1=MJ*x0+vJ >> x2=MJ*x1+vJ atd.
4
Gaussova–Seidelova metoda Jestliže při výpočtu i-té složky vektoru (k + 1)-ní aproximace, tj. čísla xk+1 , využijeme už vypoi k+1 čítané složky xk+1 , . . . , x , dostaneme tzv. Gaussovu-Seidelovu metodu. Tu dostaneme ze vzorce 1 i−1 (2) volbou MG = −(D + L)−1 U, v¯G = (D + L)−1¯b a můžeme ji vyjádřit v maticovém tvaru x ¯(k+1) = −(D + L)−1 U x ¯(k) + (D + L)−1¯b. Rozepsáním po souřadnicích dostáváme n i−1 X 1 X bi (k+1) (k) (k+1) xi =− + aij xj aij xj + , aii j=1 a ii j=i+1
i = 1, . . . , n,
k = 0, 1, . . .
(6)
Podmínky konvergence: a) ρ(MG ) < 1 ⇔ G-S metoda konverguje. Vlastní čísla matice MG lze získat řešením rovnice det(λL + λD + U ) = 0. b) Existuje norma, pro kterou kMG k < 1 ⇒ G-S metoda konverguje. c) A je ostře diagonálně dominantní ⇒ G-S metoda konverguje. d) A je symetrická a pozitivně definitní ⇒ G-S metoda konverguje. Příklad 2: Pro soustavu (5) z předchozího příkladu vypočítáme dvě iterace G-S metodou. Řešení: Protože matice soustavy A je ostře diagonálně dominantní, je splněná postačující podmínka konvergence c) a G-S metoda konverguje. Soustavu upravíme stejně jako u Jacobiovy metody x1 = − 14 ( x2 = x3 =
x2 − 2x3 − 4)
− 51 (2x1 1 5 ( x1 +
+ 2x3 − 5) − 4),
x2
ale při výpočtu užijeme vzorec (6), tj. při výpočtu druhé a třetí složky vektoru x ¯(1) už využijeme složky vypočítané v předchozích rovnicích. (0)
(1)
(0)
x1 = − 14 (x2 − 2x3 − 4) = − 41 (0 − 0 − 4) = 1 (1)
(1)
(0)
x1 + 2x3 − 5) = − 51 (2 + 0 − 5) = x2 = − 15 (2x (1)
1 5(
x3 =
(1)
(1)
(1) (1) x1 + x2 − 4) =
1 5 (1
+
3 5
3 5
− 4) = − 12 25 .
′ ′ První iterace je tedy x ¯(1) = (1, 53 , − 12 25 ) = (1; 0, 6; −0, 48) . Stejným způsobem počítáme druhou iteraci: (1)
(2)
(1)
x1 = − 41 (x2 − 2x3 − 4) = − 14 ( 53 + 2 · (2)
(2)
(1)
x1 + 2x3 − 5) = − 15 (2 · x2 = − 51 (2x (2)
x3 =
1 5(
(2)
(2)
x1 + x2 − 4) =
1 61 5 ( 100
61 100
+
12 25
− 4) =
−2·
474 500
12 25
61 100
− 5) =
474 500
− 4) = − 1221 2500 .
1221 ′ 61 474 , 500 , − 2500 ) = (0, 61; 0, 948; −0, 4884)′ . Druhá iterace je x ¯(2) = ( 100
Chybu budeme odhadovat podobně jako u Jacobiovy metody k¯ x∗ − x ¯(2) k∞ ≤
kMG k∞ 0, 75 k¯ x(k) − x ¯(k−1) k∞ = · 0, 39 = 1, 17. 1 − kMG k∞ 1 − 0, 75
Skutečná chyba je menší, protože k¯ x∗ − x ¯(2) k∞ = 0, 11. Pokud konvergují obě metody – Jacobiova i Gaussova–Seidelova, Gaussova–Seidelova metoda konverguje rychleji.
5
MATLAB – Gaussova–Seidelova metoda A – matice soustavy, b – vektor pravých stran, x(0) – počáteční přiblížení výpočet matice MG >> U=triu(A,1) >> L=tril(A,-1) >> D=A-U-L >> MG=(-1)*inv(D+L)*U výpočet vektoru vG >> vG=inv(D+L)*b výpočet jednotlivých iterací x(1) , x(2) , atd. >> x1=MG*x0+vG >> x2=MG*x1+vG atd.
Soustavy nelineárních rovnic – Newtonova metoda Uvažujeme soustavu nelineárních rovnic F (¯ x) = o¯,
(7)
neboli f1 (x1 , . . . , xn ) = 0 f2 (x1 , . . . , xn ) = 0 ... fn (x1 , . . . , xn ) = 0 a předpokládáme, že funkce f1 , . . . , fn mají spojité druhé parciální derivace. Metodu odvodíme pro dvě rovnice o dvou neznámých f1 (x1 , x2 ) = 0 f2 (x1 , x2 ) = 0. V okolí c¯ = (c1 , c2 )′ nahradíme funkce f1 , f2 Taylorovým polynomem prvního stupně ∂f1 (c1 , c2 )(x1 − c1 ) + ∂x1 ∂f2 f2 (x1 , x2 ) ≈ f2 (c1 , c2 ) + (c1 , c2 )(x1 − c1 ) + ∂x1 f1 (x1 , x2 ) ≈ f1 (c1 , c2 ) +
∂f1 (c1 , c2 )(x2 − c2 ) = 0 ∂x2 ∂f2 (c1 , c2 )(x2 − c2 ) = 0. ∂x2
Tuto soustavu zapíšeme v maticovém tvaru
∂f1 (c , c ) ∂x1 1 2 ∂f 2 (c1 , c2 ) ∂x1
6
∂f1 ! Ã ! (c1 , c2 ) Ã f1 (c1 , c2 ) x1 − c1 ∂x2 =− x −c ∂f2 f2 (c1 , c2 ) 2 2 (c1 , c2 ) ∂x2
(8)
Označíme-li F ′ tzv. Jacobiovu matici ∂f1 ∂f1 ∂x1 ∂x2 F′ = ∂f ∂f2 2 ∂x1 ∂x2 soustava (8) pak je
F ′ (¯ c) d¯ = −F (¯ c)
d¯ = x ¯ − c¯,
a
a
¯ x ¯ = c¯ + d.
Chceme-li vyřešit soustavu nelineárních rovnic Newtonovou metodou, postupujeme tak, že zvolíme nějaký počáteční vektor x ¯(0) a pro k = 0, 1, 2, . . . počítáme ve dvou krocích: i) nejprve určíme vektor d¯(k) jako řešení soustavy F ′ (¯ x(k) ) d¯(k) = −F (¯ x(k) ),
(9)
ii) položíme x ¯(k+1) = x ¯(k) + d¯(k) . Pro každý vektor x ¯(k) musí být matice F ′ (¯ x(k) ) regulární, aby soustava (9) měla právě jedno řešení. Pokud tato podmínka v některé iteraci není splněná, je potřeba zvolit jiný vektor x ¯(0) a počítat znova. Newtonova metoda konverguje rychle, ale jen když x ¯(0) je dostatečně blízko přesného řešení. Příklad 3: Pro soustavu nelineárních rovnic x21 + x22 = 4 x2 = x31 + 1 vypočítáme dvě iterace Newtonovou metodou s počátečním vektorem x ¯(0) = (1, 2)′ . Řešení: Soustavu vyjádříme ve vektorovém tvaru ¶ µ ¶ ¶ µ 2 µ 0 x1 + x22 − 4 f1 (x1 , x2 ) = o¯, = = F ′ (¯ x) = 0 −x31 + x2 − 1 f2 (x1 , x2 )
kde x ¯=
µ
x1 x2
¶
a určíme Jacobiovu matici
∂f1 ∂x 1 F ′ (¯ x) = ∂f 2
∂x1
(1)
∂f1 Ã 2x1 ∂x2 = ∂f −3x2 2
2x2 1
1
!
.
∂x2
(1)
Při výpočtu první iterace x ¯(1) = (x1 , x2 )′ : (0) (0) i) nejprve určíme vektor d¯(0) = (d1 , d2 )′ tak, že vyřešíme soustavu F ′ (¯ x(0) ) d¯(0) = −F (¯ x(0) ), po dosazení x ¯(0) = (1, 2)′ µ
2·1 −3 · 12
2·2 1
¶µ
(0)
d1 (0) d2
¶
=−
µ
12 + 22 − 4 −13 + 2 − 1
¶
⇒
µ
2 4 −3 1
¶µ
(0)
d1 (0) d2
¶
=−
µ ¶ 1 0
7
1 3 ′ a odtud dostaneme d¯(0) = (− 14 , − 14 ).
ii) ve druhém kroku první iterace vypočítáme x ¯(1) = x ¯(0) + d¯(0) , tedy (1)
x ¯
=
à ! 1 2
+
Ã
1 − 14 3 − 14
!
=
à 13 ! 14 25 14
≈
Ã
0, 9286 1, 7857
! (2)
(2)
!
⇒
Stejným způsobem postupujeme i při výpočtu druhé iterace x ¯(2) = (x1 , x2 )′ : (1) (1) i) nejprve určíme vektor d¯(1) = (d , d )′ tak, že vyřešíme soustavu 1
2
′
(1)
F (¯ x
) d¯(1) = −F (¯ x(1) ),
25 ′ po dosazení x ¯(1) = ( 13 14 , 14 ) Ã ! Ã (1) ! 2 · 25 2 · 13 d1 14 14 2 −3 · ( 13 14 )
1 Ã
⇒
13 7
− 507 196
=−
25 2 2 ( 13 14 ) + ( 14 ) − 4
25 3 −( 13 14 ) + 14 − 1 Ã 5 ! ! Ã (1) ! 25 d1 7 98 =− 41 (1) − 1 d2 2744
(1)
d2
Ã
11 ′ 105 , − 1171 ). a odtud dostaneme d¯(1) = (− 11161
ii) ve druhém kroku druhé iterace vypočítáme x ¯(2) = x ¯(1) + d¯(1) , tedy x ¯(2) =
à 13 ! 14 25 14
+
Ã
105 − 11161 11 − 1171
!
=
à 1319 ! 1435 4876 2745
≈
Ã
MATLAB – Newtonova metoda F – dané funkce (sloupcový vektor), x(0) – počáteční přiblížení označení symbolických proměnných >> syms x,y; zadání počátečního vektoru x(0) a funkce F >> x0=[1; 2] >> F=[xb 2+yb 2-4; y-xb 3-1] výpočet jacobiovy matice >> J=jacobian(F,[x,y]) (0)
(0)
dosazení x = x1 , y = x2 do jacobiovy matice a vektoru F >> J0=subs(J,[x,y],x0) >> F0=subs(F,[x,y],x0) výpočet první aproximace x(1) (d(0) = −inv(J0) ∗ F 0) >> x1=x0-inv(J0)*F0 výpočet druhé aproximace x(2) (d(1) = −inv(J1) ∗ F 1) >> J1=subs(J,[x,y],x1) >> F1=subs(F,[x,y],x1) >> x2=x1-inv(J1)*F1 atd.
8
0, 9192 1, 7763
!
.
