NUMERICKÉ METODY VÍTĚZSLAV VESELÝ PODPORA VÝUKY NA POČÍTAČOVÉ SÍTI http://www.math.muni.cz/~vesely a) MATLAB: Algoritmy: Příkazy pro připojení a výpis obsahu knihoven: d = jmathews('chap x') chap x . . . adresář korespondující s číslem kapitoly z učebnice J. H. Mathews: NUMERICAL METHODS d = vesely('nummet') . . . výpis podknihoven autora d = vesely('nummet/podknihovna') . . . výpis algoritmů v podknihovně d . . . výstup=úplná cesta do adresáře knihovny (podknihovny) Příklady: zktest([ ],'?') . . . výpis témat zktest([ ],'téma') . . . zadání příkladů zvoleného tématu. Pomocí příkazu zktest lze studovat i vzorová řešení, která jsou k některým příkladům připojena, provádět aktivní přezkoušení aj. (podrobnosti jsou v help zktest). b) Jazyky C, PASCAL, FORTRAN: Algoritmy: V knihovně d = vesely('nummet') lze nalézt podknihovny mathews.c, mathews.pas a mathews.for, které obsahují analogické algoritmy jako knihovna MATLABu jmathews zpracované po řadě v jazycích C, PASCAL a FORTRAN. Date: 3. března 2000. 1 2 VÍTĚZSLAV VESELÝ TYPICKÉ SCHÉMA NUMERICKÉ ANALÝZY ŘEŠENÉHO PROBLÉMU Obecné schéma Příklad 1. REÁLNÝ (FYZIKÁLNÍ) PROBLÉM Koule vyrobená z materiálu o hustotě = 0, 638 g/cm3 a poloměru r [cm] je ponořena do vody. Zjistěte závislost hloubky ponoru d [cm] v závislosti na r pro 8 r 12 [cm]. Reálný problém idealizujeme pomocí vhodného fyzikálně-matematického modelu a dopouštíme se tak chyby formulace problému. Předpokládáme ˇ ideálně homogenní materiál ˇ ideálně čistou vodu o hustotě 1 g/cm3 ˇ zanedbáváme vliv povrchového napětí vody 2. MATEMATICKÝ PROBLÉM Problém formulujeme jako explicitní nebo implicitní funkční vztah: y = F(x) nebo f(x, y) = 0, kde x reprezentuje vstupní data z vhodného definičního oboru a y jsou příslušná výstupní data. Spočtěme objem vody vytlačené ponořenou částí koule: MV (d) = d 0 r2 (t) dt = d 0 (r2 - (r - t)2 ) dt = d 0 (2rt - t2 ) dt = (3rd2 - d3 )/3. Podle Archimedova zákona hmotnost koule se musí rovnat hmotnosti vytlačené vody: Mk := 4r3 /3 = 1 MV (d). Obdrželi jsme tak matematický problém nalézt d splňující (implicitní) kubickou rovnici: d3 - 3rd2 + 4r3 = 0 pro x = [r, ] [8, 12] × {}, y = d. K nalezení numerického řešení je zpravidla nutné provést numerickou aproximaci (diskretizaci) matematického problému. Vzniklá chyba se nazývá chybou numerické aproximace nebo také diskretizační chybou. x ; x = [r, ] a y ; y = d, kde r := [r0, r1, . . . , rN ], d := [d0, d1, . . . , dN ], rj = 8 + j(12 - 8)/N. Nekonečný definiční obor [8, 12] × {} původní rovnice tedy nahradíme konečným diskrétním oborem pro j = 0, 1, . . . , N. 3. NUMERICKÝ PROBLÉM Po diskretizaci obdržíme: y = F(x) nebo f(x, y) = 0, kde x je vektor vstupních dat y je příslušný vektor výstupních dat. Pro každé j = 0, 1, . . . , N řešíme vzhledem k dj rovnici: d3 j - 3rjd2 j + 4r3 j = 0. Zvolíme vhodnou numerickou metodu (postup) řešení numerického problému a provedeme její algoritmizaci. Ve složitějších případech je snaha problém rozložit na jednodušší problémy se známým řešením (výstup jednoho je vstupem jiného). Zvolená metoda rovněž vnáší do řešení tzv. chybu metody. Nevhodně zvolená metoda může vést k znehodnocení výsledků v důsledku kumulace například zaokrouhlovacích chyb. Můžeme si vybrat například z těchto dvou num- erických postupů: a) spočteme analytické řešení pomocí Car- tanových vzorců. b) užijeme některou iterační metodu na řešení nelineárních rovnic, kdy počáteční odhad řešení se postupně zpřesňuje až do dosažení požadované přesnosti. NUMERICKÉ METODY 3 4. ALGORITMUS ŘEŠENÍ Počítáme y = F(x), resp. řešíme f(x, y) v konečném počtu kroků pro zadaný vstupní vek- tor x. Některé z nich mohou samy představovat dílčí algoritmy řešící některý již známý numer- ický problém. a) postupně dosazujeme do Cartanových vzorců rj pro j = 0, 1, . . . , N a vypočítáváme příslušná reálná řešení dj. b) zkonstruujeme konvergentní iterační schéma d(n) = (d(n-1) ), n = 1, 2, . . . pro vhodný počáteční odhad d(0) takové, že d(n) d pro n . Za řešení považujeme d(M) d takové, že d - d(M) < ve vhodné normě . Algoritmus kódujeme ve vhodném pro- gramovacím jazyku, což po jeho vykonání na počítači vede k výsledkům zatíženým navíc zaokrouhlovacími chybami v důsledku konečné aritmetiky počítače a chybami (ne- dokonale změřených) vstupních dat. Hustota je zjištěna fyzikálním měřením s omezenou přesností. 5. IMPLEMENTACE + VÝPOČET Programová realizace kroků algoritmu v konkrétním programovacím jazyku (C,Pascal, Fortran, MATLAB) na konkrétním typu počítače. Spočtení přibližných hodnot dj, j = 0, 1, . . . , N zatížených vlivem všech výše popsaných chyb. OBECNÉ DOPORUČENÍ PRO VOLBU NUMERICKÉ METODY Při transformaci matematického problému na numerický se vždy snažíme maximálně využít všech informací o hledaném řešení. Takto získané algoritmy jsou sice speciálnější, ale na rozdíl od univerzálně použitelných jsou obvykle efek- tivnější (rychlejší, přesnější). a) užití Cartanových vzorců představuje speciální metodu pro řešení kubické rovnice (dosazení do vzorce je jednoduché, rychlé a přesné). b) užití iterační metody představuje uni- verzální postup použitelný pro širokou třídu nelineárních rovnic (pomalá: počet iterací závisí na požadované přesnosti). NUMERICKÁ STABILITA V důsledku nevhodné kumulace chyb může algoritmus, resp. numerická metoda dávat zcela zne- hodnocené výsledky, kdy malá chyba na vstupu má za následek velkou chybu na výstupu. Definice. Algoritmus (resp. numerická metoda) se nazývá numericky stabilní, jestliže > 0 > 0 : x - x < y - y < při provádění výpočtu se zaokrouhlovacími chybami (resp. pouze s chybou numerické metody bez zaokrouhlovacích chyb). V opačném případě se algoritmus (resp. numerická metoda) nazývá numericky nestabilní a je třeba zvolit jiný algoritmus (resp. jinou numerickou metodu). 4 VÍTĚZSLAV VESELÝ 1. Úvod 1.1. Pomocná tvrzení. Označení . N. . . množina všech přirozených čísel N0 = N {0} . . . množina všech nezáporných celých čísel Z. . . množina všech celých čísel R. . . množina všech reálných čísel C. . . množina všech komplexních čísel int(x) . . . celá část čísla x (a, b) . . . otevřený interval [a, b] . . . uzavřený interval f(x), x S . . . reálná funkce f definovaná na množině S R I(x0, x1, . . . , xn) = {x | min(x0, x1, . . . , xn) < x < max(x0, x1, . . . , xn)} I[x0, x1, . . . , xn] = {x | min(x0, x1, . . . , xn) x max(x0, x1, . . . , xn)} s := v, resp. v =: s . . . označení výrazu v symbolem s. Definice (Limita funkce). L je limitou funkce f(x), x S v bodě x0 > 0 > 0 : |x - x0| < |f(x) - L| < . Píšeme lim xx0 f(x) = L lim h0 f(x0 + h) = L Definice (Spojitost funkce). Řekneme, že f(x) je spojitá v bodě x0, jestliže lim xx0 f(x) = f(x0) lim h0 f(x0 + h) = f(x0) Značíme C(S), C[a, b] . . . množina všech funkcí spojitých v každém bodě x S , resp. x [a, b]. Definice (Limita posloupnosti a součet nekonečné řady). L je limitou posloupnosti {xn} n=1 > 0 N = N() N : |xn - L| < n N. Píšeme lim n xn = L lim n (L - xn) = 0 xn L pro n . {n} n=1, n = L - xn . . . posloupnost chyb. sn = n k=1 xk . . .částečný součet nekonečné řady limn sn = s := k=1 xk . . . součet nekonečné řady. Věta 1.1. f(x), x S je spojitá v x0 {xn} n=1, xn S, limn xn = x0 platí limn f(xn) = f(x0). Věta 1.2 (Bolzanova věta). f C[a, b], L I[f(a), f(b)] c [a, b] : f(c) = L. Věta 1.3 (Weierstrasseova věta). f C[a, b] M1 = minx[a,b] f(x), M2 = maxx[a,b] f(x) a x1, x2 [a, b] : M1 = f(x1), M2 = f(x2). Důsledek 1.4. f C[a, b] nabývá na [a, b] všech hodnot mezi svou nejmenší a největší hodnotou (plyne z 1.2 a 1.3). Definice (Derivace funkce). Necht' f(x) je definována v nějakém okolí bodu x0: x (x0 - , x0 + ), > 0. NUMERICKÉ METODY 5 Pak limitu (pokud existuje) lim xx0 f(x) - f(x0) x - x0 =: f (x0) lim h0 f(x0 + h) - f(x0) h =: f (x0) nazýváme derivací funkce f(x) v bodě x0. Značíme Cn (S), Cn [a, b], n N, n 0 . . . množina všech funkcí majících v každém bodě x S , resp. x [a, b] spojité všechny derivace až do řádu n (C0 := C). Věta 1.5. f (x0) existuje f(x) je spojitá v bodě x0. Věta 1.6 (Rolleova věta). f C[a, b], f (x) existuje na (a, b), f(a) = f(b) = 0 c (a, b) : f (c) = 0. Důsledek 1.7 (Zobecněná Rolleova věta). Necht' f C[a, b], f (x), . . . , f(n) (x) existují na (a, b) a f(x0) = = f(xn) = 0 pro x0, . . . , xn [a, b]. Pak c (a, b) : f(n) (c) = 0. Důkaz. Indukcí vzhledem k n (bez újmy na obecnosti lze předpokládat x0 < x1 < < xn): 1) Pro n = 1 plyne tvrzení z 1.6 pro a = x0, b = x1. 2) Necht' n > 1 a tvrzení platí pro n - 1. Pak ci (xi-1, xi), i = 1, . . . , n : f (ci) = 0 dle 1.6. Podle indukčního předpokladu c (a, b) : f (n-1) (c) = f(n) (c) = 0. Věta 1.8 (Věta o extrémech diferencovatelné funkce). f C[a, b], f (x) existuje na (a, b), f(x1) = minx[a,b] f(x), f(x2) = maxx[a,b] f(x) xi = a nebo xi = b nebo f (xi) = 0 pro i = 1, 2. Věta 1.9 (Lagrangeova věta o střední hodnotě pro derivace). f C[a, b] a existuje vlastní nebo nevlastní f (x) na (a, b) c (a, b): f (c) = f(b) - f(a) b - a . Věta 1.10 (Věta o primitivní funkci). f C[a, b] funkce F(x) C1 [a, b] nazývaná primitivní funkce nebo také antiderivace funkce f(x) taková, že F (x) = f(x), přičemž platí b a f(x) dx = F(b) - F(a). Zejména tedy můžeme psát d dx x a f(t) dt = f(x). Věta 1.11 (1. Věta o střední hodnotě pro integrály). f C[a, b] c (a, b): b a f(x) dx = f(c)(b - a). Věta 1.12 (2. Věta o střední hodnotě pro integrály). f C[a, b], g(x) integrovatelná a neměnící znaménko na [a, b] c (a, b): b a f(x)g(x) dx = f(c) b a g(x) dx. Poznámka (Riemannova suma). Necht' f C[a, b], a = x0 < x1 < < xn = b, tk (xk-1, xk), xk = xk - xk-1. Pak n k=1 f(tk)xk n ---- b a f(x) dx. n k=1 f(tk)xk je tzv. Riemannova suma, pomocí níž lze b a f(x) dx aproximovat s libovolnou přesností, jestliže zvolíme n dostatečně velké. 6 VÍTĚZSLAV VESELÝ Věta 1.13 (Taylorova věta). f Cn+1 [a, b], x0 [a, b] pro x [a, b] = (x) I(x0, x) : f(x) = Pn(x) + Rn(x) Pn(x) = f(x0) + f (x0) 1! (x - x0) + f (x0) 2! (x - x0)2 + + f(n) (x0) n! (x - x0)n Rn(x) = f(n+1) () (n + 1)! (x - x0)n+1 . Důsledek 1.14. Polynom P(x) stupně n je jednoznačně určen svými derivacemi P(x0), P (x0), . . . , P(n) (x0) v libo- volném bodě x0. Důkaz. Položíme-li f = P, pak fn+1 (x) 0 Rn(x) 0 pro každé x0. Věta 1.15 (Taylorova věta pro funkce více proměnných). Necht' funkce f(x) := f(x1, x2, . . . , xm) má v bodě x0 = [x0 1, x0 2, . . . , x0 m] a jeho nějakém okolí spojité parciální derivace až do řádu n + 1 včetně. Pak pro libovolný bod x tohoto okolí existuje = [x0 1 + 1(x1 - x0 1), x0 2 + 2(x2 - x0 2), . . . , x0 m + m(xm - x0 m)], 0 < i < 1 tak, že f(x) = Pn(x) + Rn(x) Pn(x) = f(x0) + df(x0) 1! + d2 f(x0) 2! + + dn f(x0) n! Rn(x) = dn+1 f() (n + 1)! , kde dk f(x0) = x1 (x1 - x0 1) + x2 (x2 - x0 2) + + xn (xm - x0 m) k f(x0) pro k = 1, 2, . . .. Důsledek 1.16. Speciálně pro funkci f o m proměnných, která má spojité všechny parciální derivace až do řádu 2 v nějakém okolí bodu x0 = [x0 1, x0 2, . . . , x0 m] platí pro každé x = [x1, x2, . . . , xm] z tohoto okolí: f(x1, x2, . . . , xm) f(x0) + f(x0) x1 (x1 - x0 1) + f(x0) x2 (x2 - x0 2) + + f(x0) xn (xm - x0 m). NUMERICKÉ METODY 7 1.2. Vyhodnocení polynomu -- Hornerovo schéma. (Příklady v MATLABu: zktest([ ],'horner')) P(x) = anxn + an-1xn-1 + + a0, P(z) = ? Přímý výpočet: z2 = zz z3 = z2 z . . . zn = zn-1 z n - 1 násobení a1z a2z2 . . . anzn n násobení a0 + a1z + . . . + anzn n sečítání Celkem 2n - 1 násobení a n sečítání Princip Hornerova schématu: (. . . (( n-1 anz + an-1)z + an-2)z + + a1)z + a0 Celkem n násobení a n sečítání = teoretické minimum Věta 1.17 (Hornerovo schéma). P(x) = P1(x)(x - z) + P0, P1(x) = bnxn-1 + bn-1xn-2 + + b1, b0 := P0 = P(z), kde bn = an bk = ak + zbk+1, k = n - 1, n - 2, . . . , 1, 0 n násobení a n sečítání. Důkaz. P1(x)x + P0 = bnxn + bn-1xn-1 + . . . + b2x2 + b1x + b0 -zP1(x) = - zbnxn-1 - . . . - b3zx2 - b2zx - b1z P(x) = anxn + an-1xn-1 + . . . + a2x2 + a1x + a0 Odtud porovnáním koeficientů obdržíme požadované vztahy pro bk. Algoritmus. P : an an-1 an-2 . . . a1 a0 = ×z+ = ×z+ = . . . = ×z+ = bn bn-1 bn-2 . . . b1 b0 = P(z) Důsledek 1.18 (Rozšířené Hornerovo schéma). P0(x) := P(x) = Pn(z)(x - z)n + + P2(z)(x - z)2 + P1(z)(x - z) + P0(z), kde Pk(x) = Pk+1(x)(x - z) + Pk(z), k = 0, 1, . . . , n - 1. Důkaz. Rozklad z věty 1.17 aplikujeme n-krát postupně na P0(x), P1(x), . . . , Pn-1(x): P0(x) = (. . . ((Pn(z)(x - z) + Pn-1(z)) Pn-1(x) (x - z) + Pn-2(z)) Pn-2(x) ... (x - z) + + P1(z)) P1(x) (x - z) + P0(z). Důsledek 1.19 (Výpočet k-té derivace pomocí rozšířeného Hornerova schématu). P (k) 0 (z) = k!Pk(z) pro k = 0, 1, . . . , n. Důkaz. Pj(z)(x - z)j (k) x=z = 0 pro j = k k!Pk(z) pro j = k . 8 VÍTĚZSLAV VESELÝ Algoritmus. P0 : an an-1 an-2 . . . a2 a1 a0 = ×z+ = ×z+ = . . . = ×z+ = ×z+ = P1 : b (1) n b (1) n-1 b (1) n-2 . . . b (1) 2 b (1) 1 b (1) 0 = P0(z) = ×z+ = ×z+ = . . . = ×z+ = P2 : b (2) n b (2) n-1 b (2) n-2 . . . b (2) 2 b (2) 1 = P1(z) ... . . . ... = ×z+ = Pn : b (n) n b (n) n-1 = Pn-1(z) = Pn(z) NUMERICKÉ METODY 9 1.3. Analýza chyb. (Příklady v MATLABu: zktest([ ],'numchyby')) Definice 1.20. Necht' p je aproximace čísla p (p p). Pak Ep = p - p nazýváme absolutní chybou a Rp = p-p p , p = 0 nazýváme relativní chybou této aproximace. Píšeme též p = p p, pokud |Ep| p (tj. p [p - p, p + p]). Příklad 1.21. Rp je objektivnější indikátor, nebot' chybu vztahuje relativně k velikosti čísla: (1) x = 3, 141592, x = 3, 14 Ex = 0, 001592, Rx = 0, 000507 (2) y = 1 000 000 velké , y = 999 996 Ey = 4 velká , Ry = 0, 000004 malá (3) z = 0, 000012 malé , z = 0, 000009 Ez = 0, 000003 malá , Rz = 0, 25 velká Definice 1.22. Říkáme, že p aproximuje p na d významných cifer, jestliže d je největší nezáporné celé číslo takové, že |Rp| = |p-p| |p| < 10-d 2 = 5 × 10-(d+1) Příklad 1.23. Určíme počet významných cifer pro aproximace z příkladu 1.21: (1) Rx = 0, 000507 10-3 2 d 3 (x = 3, 14 . . . 3 platné cifry dle očekávání) (2) 0, 0000005 < Ry = 0, 000004 < 0, 000005 d = 5 (3) 0, 05 < Rz = 0, 25 < 0, 5 d = 0 DRUHY CHYB A) Chyba formulace problému (reálný problém ; matematický problém) Tuto chybu nelze zpravidla kvantifikovat, pouze můžeme prověřit, zda spočtené výsledky jsou v souladu s realitou. B) Chyba numerické aproximace (matematický problém ; numerický problém ; algoritmus) Tato chyba se podle okolností nazývá ­ Diskretizační chyba: např. při výpočtu b a f(x) dx aproximujeme f(x) po částech kon- stantní funkcí (viz Riemannovu sumu v odst. 1.1). ­ Chyba useknutím: vzniká v situacích, kdy nějaký nekonečný limitní proces ukončíme po konečném počtu kroků, např. nekonečnou řadu n=1 xn usekneme na její částečný součet N n=1 xn. Obecně lze říci, že chyby tohoto druhu vznikají při nahrazení komplikovaného matem- atického výrazu jednodušší formulí. Snažíme se přitom nalézt horní odhad absolutní nebo relativní chyby, které jsme se dopustili. Za typický příklad může posloužit přibližný výpočet funkční hodnoty f(x) Pn(x) pomocí Taylorova rozvoje 1.13 v okolí nějakého vhodného blízkého bodu x0. Potom Ef(x) = f(x) - Pn(x) = Rn(x) = f(n+1) () (n + 1)! (x - x0)n+1 , I(x, x0). Pokud známe horní odhad |f(n+1) (t)| L pro t I(x, x0), spočteme snadno horní odhad pro Ef(x) : |Ef(x)| L (n+1)! |x - x0|n+1 =: f(x). Příklad 1.24. eu = 1 + u + u2 2! + u3 3! + . . . je Taylorův rozvoj v bodě u = 0. Pro u = x2 , x 0, 1 2 , n = 4 dostáváme L = maxu[0, 1 4 ] |eu | = e 1 4 (totiž d5 du5 eu = eu ) a odtud ex2 1 + x2 + x4 2! + x6 3! + x8 4! , |Eex2 | 4 e|x10 | 5! 4 e 5!210 = 1, 045 × 10-5 pro x 0, 1 2 . 10 VÍTĚZSLAV VESELÝ Pomocí této aproximace pak přibližně spočteme integrál: p = 0, 544987104184 = 1/2 0 ex2 dx x + x3 3 + x5 2!5 + x7 3!7 + x9 4!9 1 2 0 = 0, 544986720817 = p Ep = p - p = 0, 000000383367, |Ep| < 1 2 4 e 5!210 = 4 e 5!211 . = 5, 225 × 10-6 =: p. 10-6 2 < Rp . = 7, 03442 × 10-7 < 10-5 2 aproximace na 5 významných cifer Rp Rp = Ep p p p . = 9, 587 × 10-6 aproximace na 4 významné cifry (pesimističtější než skutečnost) POZOR! Aproximace Rp Rp může být nebezpečně zkreslující pro p (a tedy i pro p) blízké k nule. Podobně jako v předchozím příkladě je často funkce f(h) nahrazena aproximací f(h), přičemž je znám horní odhad chyby ve tvaru M|hn |, M > 0. Toto vede k následující definici. Definice 1.25. (1) Necht' f(h) f(h) a M > 0 a n N tak, že |f(h) - f(h)| M|hn | pro dostatečně malá h. Pak říkáme, že f(h) aproximuje f(h) s řádem aproximace O(hn ) a píšeme f(h) = f(h) + O(hn ) pro h 0 . (2) Jestliže limn an = a, limn rn = 0 a K > 0 tak, že |an -a| < K|rn| pro dostatečně velké n, pak říkáme, že an konverguje k a s řádem konvergence O(rn) a píšeme an = a + O(rn) pro n . Poznámka 1.26. Předchozí definici lze ještě dále zobecnit v následujícím smyslu: Necht' a je konečné nebo a = , pak píšeme (1) f(x) = O(g(x)) pro x a M > 0 a > 0 : |f(x)| M|g(x)| pro x, |x - a| < Interpretace: v okolí bodu a se f(x) neliší významně od g(x). (2) f(x) = o(g(x)) pro x a limxa f(x) g(x) = 0 Interpretace: v okolí bodu a konverguje f(x) k nule rychleji, než g(x). Věta 1.27. Necht' f(h) = f(h) + O(hn ) pro h 0 a g(h) = g(h) + O(hm ) pro h 0 a r = min(n, m). Pak (1) f(h) g(h) = f(h) g(h) + O(hr ). (2) f(h) g(h) = f(h) g(h) + O(hr ), pokud f a g jsou v nějakém okolí nuly ohraničené. (3) 1 g(h) = 1 g(h) + O(hm ), pokud v nějakém okolí nuly je |g(h)| > 0 pro vhodné . (4) f(h) g(h) = f(h) g(h) + O(hr ), pokud pokud f a 1 g jsou v nějakém okolí nuly ohraničené. Důkaz. Necht' pro dosti malá h je |f(h) - f(h)| M1|h|n a |g(h) - g(h)| M2|h|m . Potom (1) |f(h) g(h) - (f(h) g(h))| |f(h) - f(h)| + |g(h) - g(h)| M1|h|n + M2|h|m = (M1|h|n-r + M2|h|m-r ) M |h|r . (2) |fg - fg| = |fg - fg + fg - fg| |g||f - f| + |f|g| - g| |g| L1 M1|h|n + |f| L2 M2|h|m (L1M1|h|n-r + L2M2|h|m-r ) M |h|r , kde skutečně f(h) je ohraničená pro dosti malá h (|f(h)| L2), nebot' |f(h) - f(h)| M2|h|n |f(h) - f(h)| 0 pro h 0 f(h)- < f(h) < f(h)+ pro dosti malá h, přičemž f(h) je ohraničená dle předpokladu. NUMERICKÉ METODY 11 (3) Platí 1/|gg| 2/2 , nebot' |g| = |g| - |g| + |g| ||g| - |g|| + |g| |g - g| + |g| /2 + |g|, kde skutečně pro dosti malá h je |g - g| /2 vzhledem k |g(h) - g(h)| 0 pro h 0. Pak |1/g - 1/g| = (1/|gg|) 2/2 |g - g| (2M2/2 ) M |h|m . (4) je přímým důsledkem (2) a (3). Věta 1.28 (Řád aproximace Taylorovým polynomem). Necht' f Cn+1 [a, b]; x0, x0 + h [a, b], potom f(x0 + h) = n k=0 f(k) (x0) k! hk + O(hn+1 ) . Důkaz. Vyjádření plyne z věty 1.13, jestliže položíme h := x - x0. Pak |Rn(x0 + h)| = |f(x0 + h) - Pn(x0 + h)| = fn+1 ((h)) (n+1)! |h|n+1 , kde (h) I[x0, x0 + h] a f(n+1) (x) je na tomto intervalu spojitá a tedy i ohraničená podle 1.3. Příklad 1.29 (O-aritmetika). Pro eh = 1 + h + h2 2! + h3 3! + O(h4 ) a cos(h) = 1 - h2 2! + h4 4! + O(h6 ) spočteme formálně řád aproximace pro jejich součet a součin. a) součet: eh + cos(h) = 1 + h + h2 2! + h3 3! + O(h4 ) + 1 - h2 2! + h4 4! + O(h6 ) = = 2 + h + h3 3! + h4 4! O(h4) + O(h4 ) + O(h6 ) O(h4) O(h4) . b) součin: eh cos(h) = 1 + h + h2 2! + h3 3! 1 - h2 2! + h4 4! 1+h+ h3 3! - h3 2! +O(h4) + 1 + h + h2 2! + h3 3! O(h6 ) O(h6) + + 1 - h2 2! + h4 4! O(h4 ) O(h4) + O(h4 )O(h6 ) O(h10) O(h4) = = 1 + h - h3 3 + O(h4 ) + O(h6 ) + O(h4 ) O(h4) . Vidíme, že tytéž výsledky bylo možno obdržet přímou aplikací obecné věty 1.27. C) Chyby vstupních dat (vstup ; výpočet programem) Zdrojem těchto chyb jsou nepřesnosti vzniklé při pořizování vstupních dat. V některých případech je znám apriori jejich horní odhad (například chyba kalibrace měřícího přístroje), jindy nikoliv (např. chyby lidského činitele). D) Strojová chyba (algoritmus ; výpočet programem) Zdrojem chyb během výpočtu je nepřesné zobrazování čísel v paměti počítače jako důsledek její konečné velikosti. (1) Strojová reprezentace celých čísel na n bitů (počítání modulo 2n ) (i) bez znaménka (unsigned integer) a 0 rozsah: 0 a 2n - 1 = (11 . . . 1 n )2 12 VÍTĚZSLAV VESELÝ (ii) se znaménkem (signed integer) a libovolné rozsah: (10 . . . 0 n )2 = -2n-1 a 2n-1 - 1 = (01 . . . 1 n )2. Současně vidíme, že prvý bit určuje znaménko: 1=minus, 0=plus. Je zřejmé, že v rámci uvedených rozsahů jsou celočíselné výpočty abso- lutně přesné, zatímco mimo ně naopak zcela chybné. Příklad 1.30 (n = 3). (i) 0 a 7 : mod 8 0 1 2 3 4 5 6 7 0 1 2 3 . . . (ii) -4 a 3 : 0 -7 -6 -5 mod 8 -4 -3 -2 -1 0 1 2 3 . . . K číslu a najdeme v modulární aritmetice číslo opačné -a snadno z rovnice: a + a (11...1)2 +1 = 2n 0 mod 2n . Odtud dostáváme -a = a + 1, kde a vznikne z a negací po bitech. Například a = 2 = (010)2 2 = (101)2 2 + 2 = (111)2 -2 = 2 + 1 = (110)2 = 6 mod 8. V počítačích se celá čísla zobrazují zpravidla v těchto přesnostech: n = 8 (1 bajt) . . . (un)signed char (1 znak) n = 16 (2 bajty) . . . (un)signed short integer (poloviční přesnost) n = 32 (4 bajty) . . . (un)signed integer (jednoduchá přesnost) n = 64 (8 bajtů) . . . (un)signed long integer (dvojnásobná přesnost) (2) Strojová reprezentace reálných čísel na n bitů Necht' q 2 značí základ číselné soustavy. V počítačích pracujeme s čísly nejčastěji v soustavách q = 2, 8, 16. Přesná reálná čísla se reprezentují v tzv. semilogaritmickém tvaru pohyblivé řádové čárky (normalizovaná mantisa + exponent): p = mantisa d1, d2d3 . . . dkdk+1 . . . ×qe , kde e Z je exponent a 1 d1 q - 1, 0 dj q - 1 (pro j > 1) jsou cifry man- tisy. Zejména v případě q = 2 je d1 = 1, takže tento bit je možno využít pro zobrazení znaménka mantisy (1=minus, 0=plus). Strojová reálná čísla se ukládají pouze s konečným počtem k cifer mantisy. Obdržíme tak přibližnou reprezentaci, která vznikne bud' pouhým odseknutím (chopping) přebývajících cifer nebo se navíc poslední k-tá cifra dk zaokrouhlí (rounding). Současně se také vhodně omezí rozsah exponentu: -emin e emax. Dostáváme tedy tyto aproxi- mace: a) p = flchop(p) = d1, d2d3 . . . dk × qe , s absolutní chybou aproximace 0 |p - p| < qe-(k-1) , b) p = flround(p) = d1, d2d3 . . . dk-1dk × qe , dk = round(dk, dk+1 . . . ) s absolutní chybou aproximace 0 |p - p| 1 2 qe-(k-1) . Příklad 1.31. (a) q = 10, k = 6 p = 22 7 = 3, 1428571428 flchop(p) = 3, 14285 × 100 flround(p) = 3, 14286 × 100 (b) q = 2, k = 5 p = 0,1 = (1, 10011001)2 × 2-4 flchop(p) = (1, 1001)2 × 2-4 flround(p) = (1, 1010)2 × 2-4 Tento příklad je současně ilustrací čísla, které má konečný počet cifer v dekadické soustavě, ale nekonečný počet cifer v binární soustavě a není tedy v paměti počítače zobrazeno přesně. V počítačích se reálná čísla zobrazují zpravidla v těchto přesnostech: NUMERICKÉ METODY 13 a) jednoduchá přesnost (4 bajty): n = 32 = 24 bitů mantisy + 8 bitů exponentu Rozsah exponentu: -27 -128 e 27 - 1 127 Dekadický rozsah: 2, 938736 × 10-39 až 1, 701412 × 1038 , kde 2, 938736 × 10-39 . = 1 × 2-128 a 1, 701412 × 1038 . = (1, 11 . . . 1)2 2 ×2127 . = 2 × 2127 = 2128 Dekadická přesnost mantisy: 2-23 . = 1, 2 × 10-7 cca 7 dekadických cifer přesnosti, což však vzhledem k příkladu 1.31(b) neznamená, že každé číslo s nejvýše 7 dekadickými ciframi musí být zobrazeno přesně. b) dvojnásobná přesnost (8 bajtů): n = 64 = 53 bitů mantisy + 11 bitů exponentu Rozsah exponentu: -210 -1024 e 210 - 1 1023 . Dekadický rozsah: 5, 562684646268003 × 10-309 až 8, 988465674311580 × 10307 , kde 5, 562684646268003 × 10-309 . = 1 × 2-1024 a 8, 988465674311580 × 10307 . = (1, 11 . . . 1)2 2 ×21023 . = 2 × 21023 = 21024 . Dekadická přesnost mantisy: 2-52 . = 2, 2 × 10-16 cca 16 dekadických cifer přesnosti. Běžně používaná binární reprezentace dle normy IEEE (např. počítače třídy PC) je v poněkud modifikovaném tvaru: p = flIEEE(p) = s e10e9e8 . . . e0d2, d3 . . . d53 , kde ˇ s1d2, d3 . . . d53 . . . mantisa se znaménkovým bitem s (1=minus, 0=plus) a dvěma binárními místy před řádovou čárkou (d1 = 1 a d2), ˇ e = (e10e9e8 . . . e0)2 . . . exponent v 11-bitové binární reprezentaci se znaménkem podle (1)(ii) v symetrickém rozsahu -(210 - 1) e 210 - 1. Tedy e10 = 1 odpovídá nezáporným hodnotám exponentu a e10 = 0 záporným hodnotám, přičemž případ e = -210 byl vyloučen, nebot' pro něj jsou všechny bity ex- ponentu e10e9e8 . . . e0 nulové a vznikla by kolize s vyjádřením čísla 0, která je dána nulovými hodnotami všech bitů IEEE reprezentace (jinak totiž tyto nulové bity by určovaly kladné číslo (10, 0 . . . 0)2 × 2-210 ). ŠÍŘENÍ CHYB PŘI VÝPOČTU A ZTRÁTA PŘESNOSTI Zdroje chyb: ˇ chyby vstupních dat (nepřesná měření), ˇ strojové chyby čísel v pohyblivé řádové čárce. Budeme vyšetřovat šíření chyb při provádění aritmetických operací. Příklad 1.32. Necht' p = p + Ep a q = q + Eq, pak dostáváme pro a) součet a rozdíl Epq = Ep Eq a Rpq = EpEq pq . Důkaz. Epq = (p q) - (p q) = p - p (q - q) = Ep Eq. b) násobení Epq = pEq + qEp + EpEq pEq + qEp a Rpq = Epq pq Rp + Rq . 14 VÍTĚZSLAV VESELÝ Důkaz. Epq = pq - pq = (p + Ep)(q + Eq) - pq = pEq + qEp + EpEq, kde EpEq 0. Pak Rpq pEq+qEp pq = p p Eq q + q q Ep p Rp + Rq, jestliže položíme p p 1 a q q 1. Důsledek 1.33. |Ep| p, |Eq| q a) |Epq| |Ep| + |Eq| p + q. b) |Epq| |pEq + qEp| |p|q + |q|p. Pro odhad absolutní a relativní chyby obecné operace dané funkčním předpisem lze užít následující obecnou větu. Věta 1.34. Necht' q = f(p1, p2, . . . , pn) =: f(p) a pi = pi + Ei pro i = 1, 2, . . . , n. Položme p = [p1, p2, . . . , pn] a p = [p1, p2, . . . , pn]. Má-li f spojité všechny parciální derivace až do řádu 2 v nějakém okolí O(p), p O(p), potom Eq = Ef(p) n i=1 f(p) pi Ei a Rq = Rf(p) n i=1 ln |f(p)| pi Ei , pokud |f| > 0 v O(p). Důkaz. Užijeme Taylorovu aproximaci z 1.16 v okolí p: Ef(p) = f(p) - f(p) n i=1 f(p) pi (pi - pi) = n i=1 f(p) pi Ei. Přitom zanedbání vyšších členů Taylorova rozvoje je oprávněné, nebot' v nich vystupují součiny typu EiEj . . . . Použijeme-li místo f funkce g(p1, p2, . . . , pn) := ln |f(p1, p2, . . . , pn)|, pak g(p) pi = f(p) pi f(p) g(p) - g(p) n i=1 f(p) pi Ei f(p) Ef(p) f(p) = Rf(p) Rf(p). Důsledek 1.35. |Ei| i pro i = 1, 2, . . . , n Eq = Ef(p) n i=1 f(p) pi i . Poznámka 1.36. Vztahy z příkladu 1.32 dostaneme ihned z věty 1.34, jestliže položíme f(p, q) = p+q, resp. f(p, q) = pq. Podobně odvodíme ještě vztahy pro podíl, jestliže f(p, q) = p q : Ep q 1 q Ep - p q2 Eq = qEp - pEq q2 |Ep q | |p|q + |q|p q2 = pq q2 , Rp q Rp q qEp - pEq q2 q p = Ep p - Eq q = Rp - Rq Rp - Rq . Poznámka 1.37 (Ztráta přesnosti při odečítání blízkých čísel). Podle 1.32a) je Rp-q |Ep|+|Eq| |p-q| , kde p - q 0 |Rp-q| . Situaci demonstruje následující příklad: p = 3, 1415926536, q = 3, 1415957341 mají 11 významných cifer, tj. podle 1.22 platí |Rp|, |Rq| < 5 × 10-12 . Avšak p - q = -0, 0000030805 má jen 4 významné cifry, tj. |Rp-q| < 5 × 10-5 . Vskutku: |Ep| = p|Rp| a |Eq| = q|Rq| |Ep + Eq| |Ep| + |Eq| < max(p, q) × 2 × 5 × 10-12 = q × 10-11 |Rp-q| = |Ep+Eq| |p-q| < q×10-11 3,0805×10-6 < 3,15×10-11 3×10-6 = 1,05 × 10-5 < 5 × 10-5 . Příklad 1.38. Tento příklad ilustruje jednostrannou kumulaci chyb při opakovaném sečítání na počítači pracujícím s přesností na 9 dekadických cifer (absolutní chyby se sečítají podle 1.32a) ): 100000 k=1 0,1 = 9999,9947. Důvodem je skutečnost, že číslo 0,1 má podle 1.31(b) nekonečně mnoho binárních cifer a jeho strojová reprezentace je tedy pouze přibližná. NUMERICKÉ METODY 15 2. Řešení nelineárních rovnic (Algoritmy v MATLABu: vesely('nummet/iterace'), jmathews('chap 2')) Formulace problému: Je dána rovnice f(x) = 0, resp. g(x) = x, x [a, b]. (2.1) Nalezněte (alespoň jedno) řešení, tj. p [a, b] takové, že f(p) = 0, resp. g(p) = p. Řešení p rovnice (2.1) nazýváme kořenem funkce f, resp. pevným bodem funkce g. Poznámka. f(x) = 0 g(x) = x, kde g(x) = x + f(x) nebo g(x) = x - f(x). Postup řešení: a) Separace kořenů: Nalezneme dělení a = a0 < a1 < < aK = b tak, že v každém subin- tervalu (ak, ak+1), k = 0, 1, . . . , K - 1 leží právě jedno řešení rovnice (2.1). Postup: zvolíme dostatečně hustou sít' xn = a+nx, n = 0, 1, . . . , N, x = b-a N a vykreslíme graf funkce f(x), resp. g(x) - x na intervalu [a, b]. Například v MATLABu stačí provést: dx=(b-a)/N; x=a:dx:b; plot(x,f(x)); Z grafu odhadneme přibližně pozice dělících bodů ak, které separují průsečíky grafu s osou x. Musíme ovšem dávat pozor, aby x bylo dostatečně malé. Jinak by se mohlo stát, že mineme kořeny ležící ve vzdálenosti menší než x. Místo hrubé vizuální lokalizace kořenů lze postupovat i exaktněji. Například knihovna jmathews nabízí funkci approot.m (+ demo a2 4.m), která vrací přibližné pozice kořenů r1, . . . , rM . Je zde použit následující jednoduchý algoritmus pro nalezení rj: rj = 1 2 (xn+1 + xn) pro některé n 1, je-li splněna některá z těchto podmínek (yn := f(xn)): (i) yn-1yn < 0 (funkce mění znaménko) nebo (ii) |yn| < a (yn - yn-1)(yn+1 - yn) < 0 (funkční hodnota v xn je malá a tečna mění znaménko -- snaha zachytit dvojnásobný kořen, který je bodem dotyku grafu s osou x). Po provedené separaci pak v každém subintervalu [ak, ak+1] zpřesňujeme pozici (jediného) kořene vhodnou iterační metodou (viz b)). Nadále tedy můžeme předpokládat, že p [a, b] je jediné řešení rovnice (2.1). b) Nalezení řešení p [a, b]: Zkonstruujeme vhodnou konvergentní posloupnost {pn} n=0, pn p, kde p0 je počáteční odhad ad a). Pro dostatečně velké N pak můžeme pN považovat za dostatečně přesnou aproximaci hledaného kořene. Definice 2.1. Řekneme, že kořen p rovnice f(x) = 0 má násobnost M (M N), jestliže existuje funkce h(x) taková, že f(x) = (x - p)M h(x) a 0 = |h(p)| < +. Lemma 2.2. f(x) CM [a, b], f(p) = f (p) = = f(M-1) (p) = 0, f(M) (p) = 0 p je kořenem f(x) o násobnosti M. Kritéria pro ukončení iteračního procesu Necht' |En| n a |Rn| n jsou horní odhady po řadě absolutní a relativní chyby. (KAE) Test, zda absolutní chyba je menší než zadaná přesnost > 0 (a) n < , pokud známe n, jinak (b) |En| |pn - pn-1| < (nemusí být spolehlivé, pokud f(pn-1)f(pn) > 0, kdy kořen neleží v I[pn-1, pn]) (KRE) Test, zda relativní chyba je menší než zadaná přesnost > 0 (a) n < , pokud známe n, jinak (b) |Rn| 2|pn-pn-1| |pn|+|pn-1| < , kde |pn|+|pn-1| 2 |p| (nemusí být spolehlivé, pokud f(pn-1)f(pn) > 0, kdy kořen neleží v I[pn-1, pn] nebo pn 0) (KY) Test, zda funkční hodnota je menší než zadaná tolerance > 0: |f(pn)| < (KD) Detekce dvojnásobného kořene |f(pn-1)|, |f(pn)|, |f(pn+1)| < (malé funkční hodnoty v okolí pn) a současně (f(pn) - f(pn-1))(f(pn+1) - f(pn)) < 0 (tečna mění znaménko). Lokalizace kořenů s vyšší násobností než 2 je problematická. 16 VÍTĚZSLAV VESELÝ Výše uvedená kritéria lze vhodně kombinovat užitím logických spojek disjunkce |, případně konjunkce &, například (KAE) | (KY), případně (KAE) & (KY). 2.1. Řád konvergence. Aitkenova metoda zrychlení konvergence. Definice 2.3 (Řád konvergence). Necht' pn p pro n = 0, 1, 2, . . . , en := p - pn. Jestliže existují A > 0, R > 0 tak, že lim n |p - pn+1| |p - pn|R = lim n |en+1| |en|R = A (tzv. asymptotická chybová konstanta), (2.2) pak říkáme, že {pn} n=0 konverguje k p s řádem konvergence R. Speciální případy: lineární konvergence (R = 1), kvadratická konvergence (R = 2). Poznámka. Pro velké n (2.2) |en+1| A|en|R , tj. se zvětšujícím se R roste rychlost konvergence en 0. Například pro R = 2 a |en| 10-2 můžeme očekávat |en+1| 10-4 . Pro některé posloupnosti přitom R nemusí být celé číslo. Definice 2.4. Pro {pn} n=0 nazýváme {k pn} n=0 posloupností jejích k-tých diferencí (k N), jestliže pn := pn+1 - pn, 2 pn := (pn) = pn+2 - pn+1 pn+1 -(pn+1 - pn pn ) = pn+2 - 2pn+1 + pn, ... k pn := (k-1 pn) . Věta 2.5 (Zrychlení konvergence Aitkenovou extrapolací). Necht' {pn} n=0, pn p. Jestliže platí (1) p - pn = 0 pro dosti velká n, (2) A, |A| < 1 tak, že limn p-pn+1 p-pn = A (alespoň lineární konvergence, nebot' může být A = 0), potom posloupnost {qn} n=0, kde qn := pn - (pn)2 2pn = pn - (pn+1 - pn)2 pn+2 - 2pn+1 + pn = pnpn+2 - p2 n+1 pn+2 - 2pn+1 + pn (2.3a) konverguje k p rychleji než {pn} n=0 v tom smyslu, že limn|p-qn p-pn | = 0, přičemž 2 pn = 0 pro dosti velká n. Poznámka 2.6 (Přibližné odvození vztahu (2.3a)). Dle (2) platí pro velká n přibližně: p - pn+1 p - pn A p - pn+2 p - pn+1 (p - pn+1)2 (p - pn)(p - pn+2) p2 - 2ppn+1 + p2 n+1 p2 - ppn+2 - ppn + pnpn+2. Odtud po úpravě obdržíme p pnpn+2 - p2 n+1 pn+2 - 2pn+1 + pn , což odpovídá (2.3a). Aitkenova extrapolace se někdy nazývá Aitkenův 2 -proces. Většinou alternativně používáme zrychlenou posloupnost {pn} n=0 (n 2): pn := qn-2 = pn-2pn - p2 n-1 pn - 2pn-1 + pn-2 = pn - (pn - pn-1)2 pn - 2pn-1 + pn-2 = pn - (pn-1)2 2pn-2 . (2.3b) Její formální výhodou je, že pro výpočet pn používá pouze hodnoty pj pro j n a nikoliv budoucí hodnoty jako v (2.3a). Příklady (MATLAB). zktest([ ],'itaitken') Algoritmus (MATLAB). vesely('nummet/iterace'): aitken.m NUMERICKÉ METODY 17 2.2. Metody dělení intervalu. Předpoklady: f(x) C[a, b], f(a)f(b) < 0, f(p) = 0. Definice 2.7. Necht' : R × R R je libovolná funkce dvou proměnných taková, že x1 < (x1, x2) < x2 platí pro každé x1, x2 R, x1 < x2. Konstruujeme konečnou nebo nekonečnou posloupnost intervalů {[an, bn]}n tak, že [a0, b0] [a1, b1] [an, bn] . . . , p (an, bn) pro n = 0, 1, . . . , takto: (1) a0 := a, b0 := b (2) Jestliže [an, bn] je již zkonstruováno, položíme pn := (an, bn). Mohou nastat tři případy: () je-li f(pn) = 0, pak p = pn je kořen a konstrukce končí; jinak položíme: () an+1 := an, bn+1 := pn, pokud f(an)f(pn) < 0 ( f(pn)f(bn) > 0) () an+1 := pn, bn+1 := bn, pokud f(an)f(pn) > 0 ( f(pn)f(bn) < 0) Výše popsaná konstrukce posloupnosti {p}n se nazývá metodou dělení intervalu s dělící funkcí . Věta 2.8. Jestliže bn - an 0, pak |pn - p| bn+1 - an+1 a tedy pn p. Speciální případy (1) Metoda půlení intervalu (bisekce) (x1, x2) := x1+x2 2 , tj. pn = an+bn 2 . Příklady (MATLAB). zktest([ ],'itbisek') Algoritmus (MATLAB). vesely('nummet/iterace'): itbisek.m jmathews('chap 2'): bisect.m (+ demo a2 2.m) (2) Metoda regula-falsi pn určíme jako průsečík spojnice bodů [an, f(an)] a [bn, f(bn)] s osou x: f(bn)-f(an) bn-an = f(bn)-0 bn-pn pn = bn - f(bn)(bn-an) f(bn)-f(an) . Tedy dělící funkce je v tomto případě definována předpisem: (x1, x2) := x2 - f(x2)(x2-x1) f(x2)-f(x1) . Příklady (MATLAB). zktest([ ],'itregfal') Algoritmus (MATLAB). vesely('nummet/iterace'): regfalsi.m, itregfal.m jmathews('chap 2'): regula.m (+ demo a2 3.m) (3) Metoda zlatého řezu y := (x1, x2) je určeno poměrem y-x1 x2-y = x2-x1 y-x1 , což je ekvivalentní s řešením kvadratické rovnice (y - x1)2 - (x2 - y)(x2 - x1) = 0 vzhledem k y. (4) Metoda Monte Carlo (náhodné střílení) (x1, x2) je náhodné číslo mezi x1 a x2. Algoritmus (MATLAB). x1 + (x2-x1)*rand(1) Věta 2.9 (Konvergence metody půlení intervalu). Necht' f C[a, b], f(a)f(b) < 0 a p (a, b) je jediný kořen f(p) = 0. Je-li {pn}n posloupnost metody půlení intervalu, pak |p - pn| b-a 2n+1 pro n = 0, 1, . . . a tedy limn pn = p. Důsledek 2.10. N = int(ln(b-a)-ln ln 2 ) je počet iterací metody půlení intervalu, po jehož dosažení je chyba aproximace kořene p menší než předepsaná tolerance > 0, tj. |p - pN | < . Důsledek 2.11. limn|p-pn+1 p-pn | limn b-a 2n+2 2n+1 b-a = 0, 5 (lineární konvergence). Poznámka . Metoda půlení intervalu konverguje pomalu -- podle 2.11 se po každé iteraci chyba redukuje přibližně na jednu polovinu, neboli přesnost se zvýší o jedno binární místo. Protože 10-1 (1 2 )3,3 je ke zvýšení přesnosti o jedno desetinné místo potřeba nejméně tří iterací. 18 VÍTĚZSLAV VESELÝ Věta 2.12 (Konvergence metody regula falsi). Necht' f C[a, b], f(a)f(b) < 0 a p (a, b) je jediný kořen f(p) = 0. Pak posloupnost iterací {pn}n metody regula falsi konverguje k p lineárně. Věta 2.13. Je-li f konvexní (resp. konkávní) na intervalu [an, bn] a xn {an, bn} je ten z krajních bodů, pro nějž f(xn) > 0 (resp. f(xn) < 0), pak xn zůstává v dalších iteracích metody regula falsi pevný, tj. xn = am nebo xn = bm pro každé m n. Poznámka 2.14. Dá se ukázat, že v konkávním, resp. konvexním případě je |p - pn+1| = An|p - pn|, kde An < 1 2 pro velká n. Konvergence je lineární jako u metody půlení intervalu, avšak s asymp- totickou chybovou konstantou limn An = A < 1 2 . Lze tedy očekávat, že metoda regula falsi bude mírně rychlejší než metoda půlení intervalu. Je rovněž vidět, že na rozdíl od metody půlení intervalu bm - am 0, nebot' bm - am p - an nebo bm - am bn - p pro dostatečně velká m n. 2.3. Metoda pevného bodu. Rovnice: x = g(x), x [a, b]. Definice 2.15. Bod p [a, b], pro nějž p = g(p) nazýváme pevným bodem funkce g(x). Definice 2.16 (Metoda pevného bodu). Necht' p0 R libovolné a pn+1 = g(pn) pro n = 0, 1, . . . . Pak posloupnost {pn} n=0 nazýváme posloup- ností iterací metody pevného bodu (PB-iterací) funkce g. Bod p0 nazýváme startovacím bodem této posloupnosti. Věta 2.17. Necht' g(x) je spojitá na R a {pn} n=0 je posloupností jejích PB-iterací, která je konver- gentní, p = limn pn. Pak p je pevným bodem funkce g. Věta 2.18. Necht' g C[a, b]. Potom (1) g([a, b]) [a, b] g má pevný bod v [a, b]. (2) Jestliže navíc g splňuje Lipschitzovu podmínku s konstantou 0 K < 1, tj. |g(x) - g(x )| K|x - x | x, x [a, b] (nebo existuje g na (a, b), |g (x)| K < 1 x (a, b)), pak g má jediný pevný bod v [a, b]. Věta 2.19 (Věta o pevném bodu). Necht' p je pevný bod funkce g(x) a I = [p - , p + ] [a, b], jeho -okolí ( > 0). Necht' g(x) je na I definována (resp. zde navíc má derivaci g (x)) a {pn} n=0 je její posloupnost iterací se startovacím bodem p0 I. Pak (1) |g(x) - g(p)| K|x - p|, 0 K < 1 pro x I (resp. |g (x)| K < 1, pro x I) pn p a pn I pro každé n = 0, 1, . . . . Takový pevný bod p nazýváme přitahujícím pevným bodem. (2) |g(x) - g(p)| |x - p| pro x I (resp. |g (x)| 1 pro x I) pn p, pokud p0 = p a současně (g(x) = p x = p na [a, b]). Jestliže |g(x) - g(p)| K|x - p|, K > 1, pro x I (resp. |g (x)| K > 1 pro x I), pak navíc pn I m > n : pm / I. Takový pevný bod p nazýváme odpuzujícím pevným bodem. Důsledek 2.20. Za předpokladů věty 2.19(1) platí následující odhady chyby aproximace: (1) |p - pn| Kn |p - p0| pro n 1. (2) |p - pn| Kn |p1-p0| 1-K pro n 1. (3) |p - pn| K 1-K |pn - pn-1| pro n 1. Jestliže |g (x)| K < 1 na I, pak |p - pn| |g (n-1)| 1-|g (n-1)| |pn - pn-1| K 1-K |pn - pn-1| pro vhodné n-1 I[pn-1, p]. (4) K 1 2 |p - pn| |pn - pn-1|. Poznámka 2.21. Pokud g (x) C(I), pak tvrzení (1) a (2) věty 2.19 lze vzhledem ke spojitosti g (x) nahradit po řadě jednoduššími kritérii: (1') |g (p)| < 1, (2') |g (p)| > 1. Totiž z nich pak plynou původní (1) a (2) takto: Pro vhodné > 0 lze kolem p nalézt dostatečně malé okolí O(p) I, v němž platí |g (x)| |g (p)| + < 1 v případě (1), resp. |g (x)| |g (p)| - > 1 v případě (2). Příklad 2.22. (1) Řešte x2 = c, c > 0 (i) x2 = c x = c x g(x) = c x , p1 = c p0 , p2 = c c/p0 = p0, . . . a pro p0 = c dostáváme alternující divergentní posloupnost. Totiž |g (x)| = |- c x2 | |g ( c)| = 1, |g (x)| < 1 pro x > c a |g (x)| > 1 pro x < c. NUMERICKÉ METODY 19 (ii) x2 = c 2x2 = x2 +c x = 1 2 x+ c x g(x) = 1 2 x+ c x , g (x) = 1 2 1- c x2 , g ( c) = 0. Podle 2.19(1) pn c konverguje velmi rychle s uvážením 2.20(1), nebot' konstanta K = 0. Vidíme, že volbou g(x) můžeme výrazně ovlivnit konvergenční vlastnosti příslušné posloup- nosti PB-iterací. (2) Rovnice x = 2 x - 1 má jediný pevný bod p = 2 pro x 1, přičemž g (x) = 1/ x - 1 a tedy g (2) = 1, g (x) > 1 pro 1 x < 2 a g (x) < 1 pro x > 2. Není tedy splněna žádná z podmínek (1') a (2') z poznámky 2.21, abychom mohli jednoznačně rozhodnout o konvergenci či di- vergenci. Numerický experiment (viz zktest(-1,'itpevbod')) ukazuje divergenci zvolíme-li počáteční aproximaci p0 = 1, 5 < 2 v levém okolí a naopak velmi pomalou konvergenci pro p0 = 2.5 > 2 z pravého okolí. Z (1)(i) a z (2) můžeme tedy usuzovat, že v případě g (p) = 1 posloupnost může i nemusí konvergovat. Algoritmus (MATLAB). vesely('nummet/iterace'): fixpoint.m jmathews('chap 2'): fixpt.m (+ demo a2 1.m) 2.4. Newton-Raphsonova metoda. Předpoklady: f(x) C2 [a, b], f(p) = 0, f (p) = 0, p [a, b]. Konstrukce posloupnosti iterací: ˇ p0 . . . počáteční aproximace ˇ Je-li pn již zkonstruováno, pak pn+1 nalezneme jako průsečík tečny v bodě (pn, f(pn)) s osou x: f (pn) = 0 - f(pn) pn+1 - pn . Odtud dostáváme iterační schéma Newton-Raphsonovy metody: pn+1 = pn - f(pn) f (pn) pro n = 0, 1, 2, . . . (2.4) Poznámka 2.23. Iterační schéma (2.4) lze považovat za iterace pevného bodu funkce g(x) = x - f(x) f (x) , nebot' pro x [p - , p + ] při dostatečně malém > 0 je f (x) = 0 a tedy x = x - f(x) f (x) 0 = - f(x) f (x) f(x) = 0. Spočtěme g (x): g (x) = x - f(x) f (x) = 1 - f (x)f (x) - f(x)f (x) f (x)2 = f(x)f (x) f (x)2 = f(x) f (x) f (x) f (x) . (2.5) Protože f(p) = 0 a f (p) = 0, je g (p) = 0 a posloupnost (2.4) Newton-Raphsonových iterací tedy splňuje předpoklad 2.21(1'). Dokázali jsme tak následující tvrzení. Věta 2.24 (Newton-Raphsonova věta). Necht' f C2 [a, b] má kořen p [a, b], tj. f(p) = 0. Jestliže f (p) = 0, pak existuje > 0 tak, že posloupnost (2.4) Newton-Raphsonových iterací {pn} n=0 konverguje k p pro libovolné p0 [p - , p + ], přičemž pn [p - , p + ] pro každé n = 0, 1, 2, . . . . Důsledek 2.25 (Newton-Raphsonovy iterace pro hledání druhé odmocniny). Necht' c R, c > 0 a p0 > 0 je libovolný počáteční odhad c. Definujme posloupnost {pn} n=0 rekurzivně takto: pn+1 = 1 2 pn + c pn pro n = 0, 1, 2, . . . , pak pn c při n (viz též příklad 2.22(1)(ii)). Věta 2.26 (Řád konvergence Newton-Raphsonovy metody). Necht' posloupnost Newton-Raphsonových iterací (2.4) konverguje ke kořenu p funkce f(x). Označíme-li en = p - pn pro n = 0, 1, . . . , pak za předpokladů věty 2.24 tato posloupnost konverguje kvadraticky s asymptotickou chybovou konstantou A := lim n |en+1| |en|2 = |f (p)| 2|f (p)| . 20 VÍTĚZSLAV VESELÝ Důkaz. Vezmeme prvé dva členy Taylorova rozvoje funkce f(x) v bodě pn pro x = p (viz 1.13): 0 = f(p) = f(pn) + f (pn)(p - pn) + f (n) 2! (p - pn)2 , n I(p, pn). (2.6) Po vydělení rovnice (2.6) hodnotou f (pn) a úpravě dostaneme po limitním přechodu n + požadovaný vztah. Věta 2.27 (Zrychlená konvergence Newton-Raphsonovy metody). Necht' posloupnost Newton-Raphsonových iterací (2.4) konverguje lineárně ke kořenu p násobnosti M > 1 funkce f(x). Potom modifikovaná posloupnost pn+1 = pn - M f(pn) f (pn) , n = 0, 1, . . . (2.7) konverguje kvadraticky k p. Poznámka 2.28. Z poznámky 2.23 vidíme, že Newton-Raphsonova metoda (2.4) představuje optimální postup jak řešení rovnice f(x) = 0 převést na řešení rovnice x = g(x) metodou pevného bodu, nebot' g (p) = 0 vede k rychlé konvergenci s nejmenší možnou konstantou K = 0 (viz 2.20(1)). Dobře je toto demonstrováno na případu rovnice x2 = c z příkladu 2.22(1) , kdy nevhodný postup (i) nevede ke konvergentní posloupnosti. Naopak postup (ii) zdůvodněný ve 2.25 dává rychle konvergentní posloupnost. Například pro c = 5 (hledáme 5) a p0 = 2 dostáváme: p1 = 2 + 5/2 2 = 2, 25, p2 = 2, 25 + 5/2, 25 2 = 2, 23611111, p3 = 2, 23611111 + 5/2, 23611111 2 = 2, 236067978, . . . , pk = 2, 236067978 pro k 4. Tedy již po čtyřech iteracích dostáváme 5 s přesností na 9 desetinných míst. Věta 2.29 (Odhad chyby Newton-Raphsonových iterací). Za předpokladů věty 2.24 necht' |f (x) f (x) | L pro x I[pn-1, p]. Pak platí (1) L|pn - pn-1| K < 1 |p - pn| K 1-K |pn - pn-1|, (2) L|pn - pn-1| 1 2 |p - pn| |pn - pn-1|. Uvedeme ještě dvě tvrzení, pomocí nichž lze snadno ověřit, zda počáteční aproximace p0 byla zvolena dostatečně blízko p tak, aby byla zaručena konvergence podle 2.24. Věta 2.30. Necht' {pn} n=0 je posloupnost Newton-Raphsonových iterací (2.4). Označme I0 := I[p0, p0 + 2(p1 - p0)]. Jestliže 2|p1 - p0|M |f (p0)|, kde M = maxxI0 |f (x)|, potom pn I0 pro n = 0, 1, . . . a limn pn = p, kde p je jediný kořen f(x) na intervalu I0. Věta 2.31. Necht' f (x) = 0, f (x) nemění znaménko na intervalu [a, b] (tj. f(x) je konvexní nebo konkávní na [a, b]) a f(a)f(b) < 0. Jestliže | f(a) f (a) | < b - a a současně | f(b) f (b) | < b - a (tj. iterace odstartovaná v krajním bodě intervalu [a, b] nepadne mimo tento interval), potom posloupnost New- ton-Raphsonových iterací (2.4) konverguje pro libovolnou počáteční aproximaci p0 [a, b] k některému kořenu p [a, b] funkce f(x). Příklady (MATLAB). zktest([ ],'itnewton') Algoritmus (MATLAB). vesely('nummet/iterace'): newton.m jmathews('chap 2'): newton.m (+ demo a2 5.m) 2.32 (Úskalí Newton-Raphsonovy metody). Věta 2.24 nedává konstruktivní návod k nalezení -okolí kořene p pro volbu počáteční aproximace p0. Z jejího důkazu vyplývá pouze možnost odhadnout toto okolí z grafu funkce g (x) = f(x)f (x) f (x)2 tak, aby |g (x)| K < 1. Jestliže neprověříme dostatečně předpoklady této věty, může se stát, že zkonstruovaná posloupnost {pn} n=0 bude vykazovat neočekávané chování. Typické jsou následující případy: (1) f(x) > 0 nebo f(x) < 0 pro každé x R. Tedy f(x) nemá reálný kořen a {pn} n=0 nekonverguje, například f(x) = x2 - 4x + 5 = (x - 2)2 + 1 > 0. (2) p0 je příliš daleko od kořene p. NUMERICKÉ METODY 21 a) {pn} n=0 může konvergovat k jinému kořenu (mimo uvažovaný interval). Například u funkce f(x) = cos(x) hledáme p = 2 a zvolíme p0 = 3. Potom p1 -4, 015, p2 -4, 853, . . . , pn -3 2 = -4, 71238898 (viz též zktest(-12,'itnewton')). b) {pn} n=0 je cyklicky divergentní. V posloupnosti se objevují tzv. cykly, kdy po konečném počtu kroků se členy posloup- nosti opakují bud' přesně nebo přibližně. Například f(x) = x3 - x - 3 při p0 = 0 dává p1 = -3, p2 = -1, 961538, p3 = -1, 147176, p4 = -0, 006579 0, takže uvízneme v přibližném cyklu řádu 4, tj. pn+4 pn (viz též zktest(10,'itnewton')). (3) Osa x je asymptotou funkce f(x) pro x + nebo pro x -. V tomto případě může pn + nebo pn - při f(pn) 0 a algoritmus se může zastavit s ohlášením falešného kořene. Jako příklad uvažme funkci f(x) = xe-x a p0 = 2, kdy f(x) 0 rychle pro x +, například f(p15) = 5, 36 × 10-8 (viz též zktest(13,'itnewton')). Používáme-li pro zastavení odhad relativní chyby KRE(b): 2|pn - pn-1|/(|pn| + |pn-1|), může tento vykazovat falešně malou hodnotu. Z těchto důvodů se někdy toto kritérium modifikuje tak, že v něm hodnotu |pn-1| nahradíme dostatečně malým číslem, například 10-6 a tím odhad relativní chyby u falešného kořene poněkud zvětšíme. (4) |g (x)| 1 na intervalu obsahujícím kořen. V tomto případě můžeme dostat tzv. oscilatoricky divergentní posloupnost {pn} n=0, kdy hodnoty pn jsou střídavě vlevo a vpravo od kořene p a přitom se od něj vzdalují. Například v případě funkce f(x) = arctan(x) s kořenem p = 0 je g(x) = x - (1 + x2 ) arctan(x) a g (x) = -2x arctan(x). Pro p0 = 1, 45 je |g (p0)| = 2.8044 > 1 a dostáváme postupně p1 -1, 55, p2 1, 85, p3 -2, 89, atd. (viz též zktest(14,'itnewton')). 2.5. Steffensenova metoda. Konstrukce posloupnosti Steffensenových iterací {sn} n=0: Posloupnost {sn} n=0 vznikne z posloupnosti iterací pevného bodu {pn} n=0 (definice 2.16), resp. New- ton-Raphsonových iterací (2.4) jakožto speciálního případu (viz poznámku 2.23), jejím zrychlením Aitkenovou extrapolací (věta 2.5 a poznámka 2.6). Zpravidla se používá následující schéma využívající (2.3b): s0 := p0 := p0 p1 = g(p0), p2 = g(p1), s1 := p2 p3 = g(p2), p4 = g(p3), s2 := p4 (2.8) p5 = g(p4), p6 = g(p5), s3 := p6 ... p2n-1 = g(p2n-2), p2n = g(p2n-1), sn := p2n (2.3b) = p2n - (p2n - p2n-1)2 p2n - 2p2n-1 + p2n-2 , n = 1, 2, . . . Věta 2.33. Necht' {pn} n=0 je posloupnost PB-iterací funkce g(x) a {sn} n=0 příslušná posloupnost Steffensenových iterací (2.8). Potom {sn} n=0 je posloupností PB-iterací funkce G(x) = x - (g(x) - x)2 g(g(x)) - 2g(x) + x = x g(g(x)) - g2 (x) g(g(x)) - 2g(x) + x . (2.9) Důsledek 2.34 (Steffensenovy iterace pro rovnici f(x) = 0). Bud' dána rovnice f(x) = 0. Necht' {s+ n } n=0 a {s- n } n=0 jsou posloupnosti Steffensenových iterací příslušné PB-iteracím po řadě funkcí g+ (x) := x + f(x) a g- (x) := x - f(x). Pak pro n = 0, 1, . . . platí s+ n+1 = s+ n - f2 (s+ n ) f(s+ n + f(s+ n )) - f(s+ n ) (2.10a) s- n+1 = s- n - f2 (s- n ) f(s- n ) - f(s- n - f(s- n )) (2.10b) Věta 2.35 (Řád konvergence Steffensenovy metody). Pro posloupnost Steffensenových iterací {pn} n=0, kde pn := s+ n nebo pn := s- n , platí tvrzení věty 2.24 22 VÍTĚZSLAV VESELÝ za týchž předpokladů o funkci f(x), přičemž konvergence je opět kvadratická s asymptotickou chybovou konstantou A poněkud větší než ve větě 2.26: A := lim n |en+1| |en|2 = |f (p)| 2|f (p)| |1 + f (p)|, kde en := p - pn. Příklady (MATLAB). zktest([ ],'itsteff') Algoritmus (MATLAB). vesely('nummet/iterace'): steffplu.m, steffmin.m jmathews('chap 2'): steff.m (+ demo a2 7.m) 2.6. Kvazinewtonovy metody. Obecný princip konstrukce posloupnosti iterací: Necht' f(x) je opět definována na intervalu [a, b]. V Newtonově metodě pn+1 = pn - f(pn) f (pn) aprox- imujeme f (pn) (směrnice tečny v bodě (pn, f(pn))) směrnicí přímky procházející body (pn, f(pn)) a (pn, f(pn)), kde pn zvolíme vhodně v okolí pn, tj. f (pn) f(pn)-f(pn) pn-pn . Iterace kvazinewtonových metod mají tedy obecně tvar: pn+1 = pn - f(pn) pn - pn f(pn) - f(pn) pro n = 0, 1, . . . . (2.11) 2.36 (Speciální případy). (A) Metoda regula-falsi (viz odst. 2.2) p0 := a, p1 := b (počáteční iterace) pn := ps, pro n = 1, 2, . . . , kde s, 0 s < n, je největší index takový, že f(pn)f(ps) < 0. (B) Steffensenova metoda (viz tvrzení 2.34) Iterace (2.10a) a (2.10b) lze přepsat na tvar: s+ n+1 = s+ n - f(s+ n ) -f(s+ n ) f(s+ n ) - f(s+ n + f(s+ n )) = s+ n - f(s+ n ) s+ n - (s+ n + f(s+ n )) f(s+ n ) - f(s+ n + f(s+ n )) , což odpovídá při pn := s+ n kvazinewtonově metodě s volbou pn := pn + f(pn) 0 . Podobně: s- n+1 = s- n - f(s- n ) f(s- n ) f(s- n ) - f(s- n - f(s- n )) = s- n - f(s- n ) s- n - (s- n - f(s- n )) f(s- n ) - f(s- n - f(s- n )) , což odpovídá při pn := s- n kvazinewtonově metodě s volbou pn := pn - f(pn) 0 . (C) Metoda sečen p0 := a, p1 := b (počáteční iterace) pn := pn-1, pro n = 1, 2, . . . , tj. obdržíme pn+1 = pn - f(pn) pn-pn-1 f(pn)-f(pn-1) pro n = 1, 2, . . . . Příklady (MATLAB). zktest([ ],'itsecant') Algoritmus (MATLAB). vesely('nummet/iterace'): secant.m jmathews('chap 2'): secant.m (+ demo a2 6.m) Věta 2.37 (Konvergence metody sečen). Necht' f C2 [a, b] má kořen p [a, b], tj. f(p) = 0. Jestliže f (p) = 0, pak existuje > 0 tak, že posloupnost iterací metody sečen {pn} n=0 konverguje k p při libovolné počáteční aproximaci p0, p1 [p - , p + ] s řádem konvergence R = 1 2 (1 + 5) 1, 618 a asymptotickou chybovou konstantou: A := lim n |en+1| |en|R = 1 2 f (p) f (p) 1 R , kde en := p - pn. NUMERICKÉ METODY 23 Poznámka 2.38 (Přehled konvergenčních vlastností (kvazi)newtonových metod). Metoda Výpočetní náročnost jedné iterace Řád konvergence Newton-Raphsonova f(pn), f (pn) 2 Steffensenova f(pn), f(pn f(pn)) 2 sečen f(pn) 1,618 regula-falsi f(pn) 1 Dá se ukázat, že neexistuje kvadraticky konvergentní iterační metoda vyžadující pouze jedno funkční vyhodnocení f(pn) v každé iteraci. 24 VÍTĚZSLAV VESELÝ 2.7. Systémy nelineárních rovnic. f1(x1, x2, . . . , xm) = 0 x1 = g1(x1, x2, . . . , xm) f2(x1, x2, . . . , xm) = 0 x2 = g2(x1, x2, . . . , xm) ... nebo ... fm(x1, x2, . . . , xm) = 0 xm = gm(x1, x2, . . . , xm) vektorový zápis f(x) = 0 x = g(x) Příklady (MATLAB). zktest([ ],'itnelrov') 2.7.1. Metoda pevného bodu. Hledáme p := [p1, p2, . . . , pm] : p = g(p). m = 1 (viz poznámku 2.21(1')) m > 1 a) g spojitá v okolí p a) Spojité parciální derivace gi xj v okolí p pro i, j = 1, 2, . . . , m b) |g (p)| < 1 b) gi x1 (p) + + gi xm (p) < 1 pro i = 1, 2, . . . , m Pro p0 dostatečně blízko k p plati: Pro p(0) = [p (0) 1 , p (0) 2 , . . . , p (0) m ] dostatečně blízko k p platí: pn p pro n , kde pn+1 = g(pn), n = 0, 1, 2, . . . p(n) p pro n , kde p(n+1) = g(p(n) ), n = 0, 1, 2, . . . p (n+1) 1 = g1(p (n) 1 , p (n) 2 , . . . , p (n) m ) p (n+1) 2 = g2(p (n) 1 , p (n) 2 , . . . , p (n) m ) ... ... p (n+1) m = gm(p (n) 1 , p (n) 2 , . . . , p (n) m ) 2.7.2. Seidelova metoda (vylepšená metoda pevného bodu). p (n+1) 1 = g1(p (n) 1 , p (n) 2 , . . . , p (n) m ) p (n+1) 2 = g2(p (n+1) 1 , p (n) 2 , . . . , p (n) m ) ... ... p (n+1) m = gm(p (n+1) 1 , p (n+1) 2 , . . . , p (n) m ) p (n+1) i = g1(p (n+1) 1 , . . . , p (n+1) i-1 , p (n) i , . . . , p (n) m ) pro i = 1, 2, . . . , m. Algoritmus (MATLAB). jmathews('chap 2'): fix2dim.m (+ demo a2 9.m) NUMERICKÉ METODY 25 2.7.3. Newton-Raphsonova metoda. Hledáme p := [p1, p2, . . . , pm] : f(p) = 0. m = 1 (viz větu 2.24) m > 1 a) f spojitá v okolí p a) Spojité druhé parciální derivace 2 gi xj xk v okolí p pro i, j, k = 1, 2, . . . , m b) f (p) = 0 b) Jakobián (matice m × m) J(p) := f1 x1 (p) . . . f1 xm (p) ... . . . ... fm x1 (p) . . . fm xm (p) je regulární matice. Pro p0 dostatečně blízko k p plati: Pro p(0) = [p (0) 1 , p (0) 2 , . . . , p (0) m ] dostatečně blízko k p platí: pn p pro n , kde pn+1 = pn - f(pn) f (pn) f (pn)-1f(pn) , n = 0, 1, 2, . . . p(n) p pro n , kde p(n+1) = p(n) - J(p(n) )-1 f(p(n) ), n = 0, 1, 2, . . . p (n+1) 1 = p (n) 1 - m j=1 a1,jfj(p (n) 1 , . . . , p (n) m ) p (n+1) 2 = p (n) 2 - m j=1 a2,jfj(p (n) 1 , . . . , p (n) m ) ... ... p (n+1) m = p (n) m - m j=1 am,jfj(p (n) 1 ,. . . , p (n) m ) , kde J(p(n) )-1 =: [ai,j]. Algoritmus (MATLAB). jmathews('chap 2'): new2dim.m (+ demo a2 10.m) Příklad 2.39 (m = 2, x := x1, y := x2). x2 - 2x - y + 0, 5 = 0 . . . parabola x2 + 4y2 - 4 = 0 . . . elipsa Průsečíky obou křivek určují dvě řešení: Řešení (1): p = [-0, 2; 1, 0] Řešení (2): p = [1, 9; 0, 3]. Demonstrace nalezení numerického řešení Seidelovou a Newtonovou metodou v MATLABu: viz jmathews('chap 2'): a2 9.m, a2 10.m. Řešení (1) (Seidelovou) metodou pevného bodu pro p(0) = [0; 1] Od 2. rovnice odečteme 8y, abychom zajistili splnění podmínky 2.7.1 b) pro i = 2. Obdržíme ekvivalentní systém pro iterace pevného bodu ve tvaru: x = x2 - y + 0, 5 2 y = -x2 - 4y2 + 8y + 4 8 . Řešení (2) (Seidelovou) metodou pevného bodu pro p(0) = [2; 0] Od 1. rovnice odečteme 2x a od 2. rovnice odečteme 11y, abychom zajistili splnění podmínky 2.7.1 b) pro i = 1, 2 -- předchozí volba totiž v tomto kořenu tuto podmínku nesplňuje. Obdržíme 26 VÍTĚZSLAV VESELÝ ekvivalentní systém pro iterace pevného bodu ve tvaru: x = -x2 + 4x + y - 0, 5 2 y = -x2 - 4y2 + 11y + 4 11 . Řešení (2) Newtonovou metodou, např. pro [x0, y0] := p(0) = [2; 0, 25] J(x, y) = 2x - 2 -1 2x 8y , xn+1 yn+1 = xn yn - J(xn, yn)-1 x2 n - 2xn - yn + 0, 5 x2 n + 4y2 n - 4 . NUMERICKÉ METODY 27 3. Řešení systému lineárních rovnic (Algoritmy v MATLABu: vesely('nummet/linrov'), jmathews('chap 3')) Formulace problému: Je dán systém m lineárních rovnic (SLR) o n neznámých x1, x2, . . . , xn ve tvaru: a1,1x1 + a1,2x2 + . . . + a1,nxn = b1 . . . 1. rovnice a2,1x1 + a2,2x2 + . . . + a2,nxn = b2 . . . 2. rovnice ... ... . . . ... ... . . . ... am,1x1 + am,2x2 + . . . + am,nxn = bm . . . m-tá rovnice . (3.1a) Ekvivalentní maticový zápis systému (3.1a): Ax = b, (3.1b) kde je A := [ai,j]1im 1jn . . . matice rozměru m × n nazývaná maticí SLR (3.1a), ai,j R b := b1 b2 ... bm . . . pravá strana SLR (3.1a), bi R x := x1 x2 ... xn . . . vektor neznámých SLR (3.1a), xj R. Každý vektor x R vyhovující všem rovnicím SLR (3.1a) nazýváme jeho řešením. Poznámka 3.1. (1) Je-li A matice rozměru m × n, pak značíme: a:,j . . . j-tý sloupec matice A, ai,: . . . i-tý řádek matice A a píšeme A = [a:,1, a:,2, . . . , a:,n] = [a1,:; a2,:; . . . ; am,:] = a1,: a2,: ... am,: . Je Ax = a:,1x1 + a:,2x2 + + a:,nxn R(A), kde R(A) značí vektorový podprostor v Rm generovaný všemi sloupci matice A (jeho prvky jsou všechny možné lineární kombinace těchto sloupců). Je-li m = n, pak matice A se nazývá čtvercovou maticí řádu n. Je-li ai,j = 0 pro i < j (resp. i > j), pak matice A se nazývá dolní (resp. horní) trojúhelníkovou maticí. (2) r(A) značí hodnost matice A. (3) Im := [i,j] . . . jednotková matice m × m. (4) p = [p(1), p(2), . . . , p(m)]T určuje permutaci p na vzestupně uspořádané množině {1, 2, . . . , m}, tj. i udává pozici, na kterou je permutací p přemístěn prvek j := p(i), i = 1, 2, . . . , m. Permutaci lze ekvivalentně popsat permutační maticí P := [pi,j], kde pi,j := p(i),j, 1 i, j m. Zřejmě platí A[1, 2, . . . , m]T = p, takže pro libovolnou matici A o m řádcích a matici B o m sloupcích dostáváme: P A = [ap(1),:; ap(2),:; . . . ; ap(m),:] . . . p(i)-tý řádek A je přemístěn na i-tou pozici, BP T = [b:,p(1), b:,p(2), . . . , b:,p(m)] . . . p(j)-tý sloupec B je přemístěn na j-tou pozici, nebot' BP T = (P BT )T . Odtud zejména dostáváme P Im = P , takže permutační matice P se dostane permutací řádků jednotkové matice. Je-li D := diag(d1, d2, . . . , dm) diagonální matice, pak snadno nahlédneme, že P DP T = diag(dp(1), dp(2), . . . , dp(m)) je rovněž diagonální maticí, jejíž prvky na diagonále jsou permutací p diagonálních prvků původní matice. Jsou-li všechny diagonální prvky v D 28 VÍTĚZSLAV VESELÝ stejné, pak zřejmě P DP T = D. Speciálně tedy P P T = P ImP T = Im, takže P je orto- gonální matice a tedy P -1 = P T je permutační matice příslušná k inverzní permutaci p-1 . Definice 3.2. Necht' A je matice SLR (3.1b), pak matici A := [A, b] rozměru m × (n + 1) nazýváme rozšířenou maticí tohoto systému. Poznámka 3.3. Každý SLR je jednoznačně popsán svou rozšířenou maticí soustavy A. Uvážíme-li poznámku 3.1(1), můžeme vyslovit následující tvrzení o existenci řešení SLR: SLR daný maticí A má řešení b R(A) b je lineární kombinací sloupců matice A matice A a A mají stejnou hodnost. SLR má přitom nejvýše jedno řešení sloupce matice A jsou lineárně nezávislé (tvoří bázi). Obdrželi jsme tak jednoduchou úvahou známé tvrzení z lineární algebry. Definice 3.4. Řekneme, že dva SLR Ax = b a A x = b jsou ekvivalentní, jestliže oba systémy mají tytéž množiny řešení. V takovém případě píšeme A A , kde A := [A, b] a A := [A , b] značí příslušné rozšířené matice těchto systémů (s týmž počtem sloupců). Lemma 3.5. Necht' A vznikne z A provedením konečného počtu elementárních operací jednoho z následujících tří typů: ˇ záměna dvou řádků, ˇ vynásobení řádku nenulovou konstantou, ˇ přičtením lineární kombinace některých řádků k jinému řádku. Potom A A . 3.1. Přímé metody pro řešení systému lineárních rovnic. Příklady (MATLAB). zktest([ ],'linrov') Věta 3.6 (Algoritmus dopředného a zpětného dosazování). Necht' Ax = b je SLR, kde A je čtvercová dolní nebo horní trojúhelníková matice řádu n s nenulovými prvky na diagonále (ak,k = 0 pro k = 1, 2, . . . , n). Pak existuje právě jedno řešení x tohoto systému, které nalezneme pomocí algoritmu dopředného nebo zpětného dosazování. Algoritmus dopředného dosazování pro dolní trojúhelníkovou matici A: x1 = b1/a1,1 xk = bk - k-1 j=1 ak,jxj /ak,k pro k = 2, 3, . . . , n. (3.2a) Algoritmus zpětného dosazování pro horní trojúhelníkovou matici A: xn = bn/an,n xk = bk - n j=k+1 ak,jxj /ak,k pro k = n - 1, n - 2, . . . , 1. (3.2b) 3.7. Výpočetní složitost algoritmů dopředného a zpětného dosazování. Označíme-li O() počet operací sečítání a odčítání (aditivní složitost) a O(, /) počet operací násobení a dělení (multiplikativní složitost), pak užitím známého vztahu pro součet členů aritmet- ické posloupnosti dostáváme pro algoritmy dopředného a zpětného dosazování kvadratickou aditivní i multiplikativní složitost: Aditivní složitost: O() = 0 + 1 + . . . n - 1 = n(n - 1) 2 = n2 - n 2 . Multiplikativní složitost: O(, /) = 1 + 2 + . . . n = n(n + 1) 2 = n2 + n 2 . Pokud matice A má jednotkovou diagonálu (ak,k = 1 pro k = 1, 2, . . . , n), pak v každém kroku ušetříme jedno dělení a multiplikativní i aditivní složitost bude stejná: O(, /) = 0 + 1 + . . . n - 1 = O() = n2 - n 2 . NUMERICKÉ METODY 29 Algoritmus (MATLAB). vesely('nummet/linrov'): zpetdos.m jmathews('chap 3'): backsub.m (+ demo a3 1.m) Věta 3.8 (Gaussova eliminace s výběrem pivotního prvku). Necht' Ax = b je SLR se čtvercovou regulární maticí soustavy A řádu n. Pak existuje ekvivalentní SLR Ux = y, [U, y] [A, b] se čtvercovou horní trojúhelníkovou maticí soustavy U řádu n a nenulovou diagonálou (uk,k = 0 pro k = 1, 2, . . . , n). Důkaz: Algoritmus Gaussovy eliminace s výběrem pivotního prvku. Postupně pro k = 1, 2, . . . , n konstruujeme posloupnost ekvivalentních systémů: [A, b] = U1 Un = [U, y], kde U1 := [A, b] =: [u (1) j,i ] Uk := a (1) 1,1 a (1) 1,2 . . . a (1) 1,k-1 a (1) 1,k . . . a (1) 1,n a (1) 1,n+1 0 a (2) 2,2 . . . a (2) 2,k-1 a (2) 2,k . . . a (2) 2,n a (2) 2,n+1 0 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 . . . 0 u (k) k,k . . . u (k) k,n u (k) k,n+1 0 0 . . . 0 u (k) k+1,k . . . u (k) k+1,n u (k) k+1,n+1 ... ... . . . ... ... . . . ... ... 0 0 . . . 0 u (k) n,k . . . u (k) n,n u (k) n,n+1 =: [Ak, bk] [A, b], k = 2, . . . , n a(n) n,n := u(n) n,n, a (n) n,n+1 := u (n) n,n+1 U := An = [a (j) j,i ], y := bn = [a (j) j,n+1]. k-tý krok algoritmu Gaussovy eliminace, k = 1, 2, . . . , n - 1: Konstrukce Uk+1 Uk, je-li Uk již zkonstruováno. Provedeme vhodné elementární transformace dle lemmatu 3.5: Výběr pivotního prvku: j k : u (k) j,k = 0, nebot' v opačném případě by determinant det(Ak) = a (1) 1,1a (2) 2,2 . . . a (k-1) k-1,k-1det([u (k) j,i ]) = 0, což je v rozporu s [Ak, bk] [A, b]. Necht' p k je nejmenší řádkový index v Uk, kde |u (k) p,k| je maximální, tj. |u (k) p,k| = max(|u (k) k,k|, |u (k) k+1,k|, . . . |u (k) n,k|). Prvek u (k) p,k je tedy nenulový a nazývá se pivotní prvek. Záměna p-tého a k-tého řádku: Necht' P k značí permutaci, která zamění p-tý a k-tý řádek v Uk. Pivotní prvek se tak stane diagonálním, přičemž pro j k necht' a (k) j,: značí j-tý řádek v takto permutované matici P kUk. Eliminace prvků v k-tém sloupci pod pivotním prvkem: Pro j = k + 1, . . . , n přičteme k j-tému řádku vhodný násobek k-ho řádku tak, aby prvek na (j, k)-té pozici se vynuloval: u (k+1) j,: := a (k) j,: - mj,ka (k) k,: , kde mj,k := a (k) j,k a (k) k,k . (3.3) Skutečně dostáváme u (k+1) j,k = 0, přičemž pro i < k zůstává u (k+1) j,i = 0, nebot' a (k) j,i = 0 i a (k) k,i = 0. Eliminaci (3.3) můžeme užitím syntaxe MATLABu algoritmizovat takto: for j = k + 1 : n mj,k = a (k) j,k /a (k) k,k; a (k+1) j,k = 0; for i = k + 1 : n + 1 u (k+1) j,i = a (k) j,i - mj,k a (k) k,i ; end end (3.4) Skutečnost, že při výpočtu mj,k se dělí právě pivotním prvkem a (k) k,k, který byl zvolen v absolutní hodnotě největší možný, příznivě ovlivňuje numerickou stabilitu (snižuje se riziko dělení čísly blízkými k nule, což by mohlo vést k hodnotám mj,k blízkým k ). Důsledek 3.9 (LU-rozklad). P A = LU, kde matice P , L a U jsou určeny algoritmem Gaussovy eliminace: 30 VÍTĚZSLAV VESELÝ P = P n-1 . . . P 2P 1 je permutační matice popisující výslednou permutaci řádků po skončení Gaussovy eliminace, L =: [lj,k] je dolní trojúhelníková matice s jednotkovou diagonálou lj,j = 1 a lj,k = mj,k pro j > k, U =: [uj,k] je horní trojúhelníková matice po Gaussově eliminaci. Příklad 3.10 (UKÁZKA KROKŮ GAUSSOVY ELIMINACE A LU-ROZKLADU). Řešený systém lineárních rovnic: 1 2 1 4 2 0 4 3 4 2 2 1 -3 1 3 2 x = 13 28 20 6 Gaussova eliminace s pivotováním ve sloupcích pod diagonálou: krok 1: U1 = 1 2 1 4 13 2 0 4 3 28 4 2 2 1 20 -3 1 3 2 6 , p = 1 2 3 4 P 1U1 = 4 2 2 1 20 2 0 4 3 28 1 2 1 4 13 -3 1 3 2 6 , p = 3 2 1 4 , m = 1/2 1/4 -3/4 krok 2: U2 = 4 2 2 1 20 0 -1 3 5/2 18 0 3/2 1/2 15/4 8 0 5/2 9/2 11/4 21 , p = 3 2 1 4 P 2U2 = 4 2 2 1 20 0 5/2 9/2 11/4 21 0 3/2 1/2 15/4 8 0 -1 3 5/2 18 , p = 3 4 1 2 , m = 3/5 -2/5 krok 3: U3 = 4 2 2 1 20 0 5/2 9/2 11/4 21 0 0 -11/5 21/10 -23/5 0 0 24/5 18/5 132/5 , p = 3 4 1 2 P 3U3 = 4 2 2 1 20 0 5/2 9/2 11/4 21 0 0 24/5 18/5 132/5 0 0 -11/5 21/10 -23/5 , p = 3 4 2 1 , m = -11/24 NUMERICKÉ METODY 31 krok 4: U4 = 4 2 2 1 20 0 5/2 9/2 11/4 21 0 0 24/5 18/5 132/5 0 0 0 15/4 15/2 , p = 3 4 2 1 Výsledek: x = 3 -1 4 2 , m(2, 1) = 1/2 patří k 2. řádku v A - 3. řádek v P A m(3, 1) = 1/4 patří k 1. řádku v A - 4. řádek v P A m(4, 1) = -3/4 patří k 4. řádku v A - 2. řádek v P A m(3, 2) = 3/5 patří k 1. řádku v A - 4. řádek v P A m(4, 2) = -2/5 patří k 2. řádku v A - 3. řádek v P A m(4, 3) = -11/24 patří k 1. řádku v A - 4. řádek v P A Získaný LU-rozklad: P = 0 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 , L = 1 0 0 0 -3/4 1 0 0 1/2 -2/5 1 0 1/4 3/5 -11/24 1 , U = 4 2 2 1 0 5/2 9/2 11/4 0 0 24/5 18/5 0 0 0 15/4 3.11. Konstrukce matic L, U, P pomocí modifikované Gaussovy eliminace Algoritmus Gaussovy eliminace z důkazu věty 3.8 modifikujeme takto: ˇ b := [1, 2, . . . , n]T , ˇ (3.3) provádíme jen s prvými n sloupci matice Uk, přičemž místo nul klademe u (k+1) j,k = mj,k pro j = k + 1, . . . , n; neboli na místo eliminovaných prvků ukládáme multiplikativní faktory mj,k. ˇ Po poslednim kroku obdržíme matici U =: [a (j) j,i ], kde a (j) j,i = lj,i pro j > i, a (j) j,i = uj,i pro j i n a a (j) j,n+1 = pj pro j = 1, 2, . . . , n. Dostáváme dva výsledné postupy pro řešení SLR na bázi Gaussovy eliminace: 3.12. Přímý algoritmus Sestává ze dvou kroků: (1) Pomocí algoritmu Gaussovy eliminace popsaném v důkazu věty 3.8 převedeme Ax = b na ekvivalentní systém Ux = y s horní trojúhelníkovou maticí soustavy U. (2) Soustavu Ux = y řešíme algoritmem zpětného dosazování (3.2b). Algoritmus (MATLAB). vesely('nummet/linrov'): gelim.m, jelim.m (jelim.m je implementace tzv. Jordanovy metody, která převádí SLR na systém s diagonální maticí soustavy, kdy eliminace probíhá i nad diagonálou) jmathews('chap 3'): uptrbk.m (+ demo a3 2.m) 3.13. Algoritmus nepřímé faktorizace Sestává ze čtyř kroků: (1) Konstrukce matic L, U a P podle 3.11. (2) Provedení permutace pravé strany P b. (3) Řešíme soustavu Ly = P b algoritmem dopředného dosazování (3.2a). (4) Řešíme soustavu Ux = y algoritmem zpětného dosazování (3.2b). Kroky (3) a (4) jsme obdrželi takto: Ax = b P Ax = P b LUx = P b, kde položíme Ux =: y. Algoritmus (MATLAB). vesely('nummet/linrov'): lufakt.m, lureseni.m jmathews('chap 3'): lufact.m, lusolve.m (+ demo a3 3.m) 32 VÍTĚZSLAV VESELÝ 3.14. Výpočetní složitost algoritmů 3.12 a 3.13. Oba algoritmy mají stejnou kubickou aditivní i multiplikativní složitost: Aditivní složitost: O() = n3 - n 3 + n2 - n 2 Multiplikativní složitost: O(, /) = n3 - n 3 + n2 Z hlediska výpočtové složitosti jsou oba algoritmy ekvivalentní. Vzhledem k tomu, že v obou případech má nejvyšší výpočetní nároky triangularizace (převod na horní trojúhelníkový tvar) dle věty 3.8 (faktor n3 -n 3 ), bude algoritmus 3.13 výhodnější v případě, že řešíme systém opakovaně s různými pravými stranami b a stejnou maticí soustavy A. V takovém případě totiž stačí provést triangularizaci P A = LU jen jednou, nebot' nezávisí na pravé straně, a pak opakovaně provádět kroky (2) až (4) pro jednotlivé pravé strany. Tohoto postupu lze například využít při výpočtu inverze regulární čtvercové matice A řádu n, kdy AA-1 = In je ekvivalentní s řešením n systémů lineárních rovnic Ax1 = 1 0 ... 0 , Ax2 = 0 1 ... 0 , . . . , Axn = 0 0 ... 1 , kde xi reprezentuje i-tý sloupec hledané inverzní matice A-1 . 3.15. Problémy s numerickou nestabilitou při řešení SLR. Matice A se nazývá špatně podmíněná, jestliže malé změny pravé strany b mají za následek velké změny řešení x = A-1 b. Tato situace nastává v případech, kdy A je "skoro singulární", tj. některý sloupcový (řádkový) vektor v A svírá malý úhel s nadrovinou generovanou ostatními sloupci (řádky). Jako příklad uvažme následující systém dvou lineárních rovnic: x + 2y = 2 2x + 3y = 3, 4 Přesné řešení x = 0, 8, y = 0, 6 aproximujme dosti nepřesně hodnotami x = 1, y = 0, 48. Po jejich dosazení do obou rovnic dostáváme po řadě 1, 96 místo 2 a 3, 44 místo 3, 4, což jsou malé odchylky řádově v setinách ve srovnání s chybou řešení, která je řádově v desetinách. Příčina tkví v tom, že hledáme průsečík dvou přímek, které svírají velmi malý úhel o směrnicích -1 2 a -2 3 . Řádky [1, 2] a [2, 3] matice soustavy v tomto případě totiž reprezentují souřadnice kolmic k oběma přímkám a svírají tedy stejný úhel jako přímky samotné. Se zmenšujícím se úhlem se blížíme ke stavu lineární závislosti obou vektorů, tj. ke stavu, kdy matice soustavy je singulární. Podobně lze zkonstruovat situace, kdy při libovolně velké odchylce řešení lze dosáhnout libovolně malé odchylky od pravé strany. ZÁVĚR: Nelze tedy hodnotit přesnost nalezeného řešení podle odchylek od pravé strany po jeho dosazení do rovnic! Lze pouze hodnotit podmíněnost matice. Ta se měří číslem podmíněnosti matice, které se zpravidla počítá jako poměr největšího ku nejmenšímu singulárnímu číslu matice (singulární čísla matice A jsou druhé odmocniny vlastních čísel symetrické matice AT A -- tedy v případě symetrické matice splývají s jejími vlastními čísly). Toto číslo podmíněnosti v systému MATLAB počítá funkce cond. Čím větší číslo podmíněnosti má daná matice, tím je hůře podmíněná a více se blíží k singulární matici. Jinak počítané reciproké číslo podmíněnosti v intervalu [0, 1] vrací funkce rcond. V tomto případě špatnou podmíněnost indikují hodnoty blízké k nule. Poznámka 3.16. Gaussovu eliminaci z věty 3.8 lze modifikovat pro případ jakékoliv (i nečtvercové) mat- ice soustavy A. Výsledkem je úprava na schodovitý horní trojúhelníkový tvar. Zatímco v případě, kdy sloupce A jsou lineárně nezávislé, dostáváme "šířku" schodů rovnou jedné jako v případě věty 3.8, obecně se může stát, že některé "schody" budou širší (při eliminaci nelze vybrat nenulový pivotní prvek, nebot' příslušná poddiagonálová část sloupce je nulová -- v takovém případě jej ponecháme a přejdeme na další sloupec, atd.). NUMERICKÉ METODY 33 Dostaneme-li po skončení schodovitou matici Un = [An, y], kde je "schod" mezi An a y, pak systém nemá řešení, nebot' poslední nenulový prvek vektoru y nemůžeme dostat jako lineární kombinaci nul z posledního řádku matice An. To znamená, že hodnost r(Un) = r(An) + 1. Pokud tato situace nenastane, pak xi korespondující s nadbytečnými sloupci v širších "schodech" odpovídají volným proměnným. Lineární kombinace těchto sloupců s těmito volnými proměnnými převedeme na pravou stranu. Tím se U stane opět standardní horní trojúhelníkovou maticí a systém dořešíme zpětným dosazováním jako v regulárním případě. Věta 3.17. Necht' Ax = b, kde A je regulární pásová tridiagonální matice (ai,j = 0 pro |i - j| > 1) rozměru n × n: A = a1 c1 0 . . . 0 0 0 b2 a2 c2 . . . 0 0 0 ... ... ... ... ... ... ... 0 0 0 . . . bn-1 an-1 cn-1 0 0 0 . . . 0 bn an . Potom její LU-rozklad A = LU je tvaru L = 1 0 . . . 0 0 2 1 . . . 0 0 ... ... ... ... ... 0 0 . . . 1 0 0 0 . . . n 1 , U = 1 c1 0 . . . 0 0 0 2 c2 . . . 0 0 ... ... ... ... ... ... 0 0 0 . . . n-1 cn-1 0 0 0 . . . 0 n , kde 1 = a1 i = bi i-1 i = ai - ici-1 pro i = 2, 3, . . . , n. (3.5a) Poznámka 3.18. (1) Jinou variantu algoritmu (3.5a) lze odvodit pro L = 1 0 . . . 0 0 b2 2 . . . 0 0 ... ... ... ... ... 0 0 . . . n-1 0 0 0 . . . bn n , U = 1 1 0 . . . 0 0 0 1 2 . . . 0 0 ... ... ... ... ... ... 0 0 0 . . . 1 n-1 0 0 0 . . . 0 1 , kde 1 = a1, 1 = c1 1 i = ai - bii-1 i = ci i pro i = 2, 3, . . . , n. (3.5b) (2) Analogické algoritmy obdržíme pro blokově tridiagonální matice, kde ai, bi, ci, i, i (resp. i) jsou submatice kompatibilních rozměrů. Pak místo 1 i-1 násobíme inverzní maticí -1 i-1 zprava v (3.5a) (resp. místo 1 i maticí -1 i zleva v (3.5b)). (3) V případě obecné regulární čtvercové matice A lze výše popsaný postup zobecnit pro přímý výpočet prvků matic L a U tak, že postupně sestavujeme rovnice pro prvky 1. řádku U, 1. sloupce L, 2. řádku U, 2. sloupce L, atd. 34 VÍTĚZSLAV VESELÝ 3.2. Iterační metody pro řešení systému lineárních rovnic. Příklady (MATLAB). zktest([ ],'itlinrov') Princip: Rovnici Ax = b, kde A je čtvercová regulární matice řádu n, převedeme na rovnici tvaru x = x + a provádíme iterace metodou pevného bodu 2.7.1 nebo Seidelovou metodou 2.7.2 pro vektorovou lineární funkci g(x) := x + o n proměnných. Obecný postup hledání a : Matici A vhodně rozložíme do tvaru A = N - M, kde N je regulární matice. Pak platí Ax = b (N - M)x = b Nx = Mx + b x = N-1 M x + N-1 b . Iterační proces: x(0) . . . počáteční iterace x(k+1) = x(k) + pro k = 0, 1, . . . . (3.6) Definice 3.19 (Jacobiho prosté iterace). Položíme N := diag(A) = a1,1 0 . . . 0 0 a2,2 . . . 0 ... ... ... ... 0 0 . . . an,n , M := N - A = 0 -a1,2 . . . -a1,n -a2,1 0 . . . -a2,n ... ... ... ... -an,1 -an,2 . . . 0 , pak = 0 - a1,2 a1,1 . . . - a1,n a1,1 - a2,1 a2,2 0 . . . - a2,n a2,2 ... ... ... ... - an,1 an,n - an,2 an,n . . . 0 , = b1 a1,1 b2 a2,2 ... bn an,n . Rozepsáním (3.6) do jednotlivých rovnic dostáváme Jacobiho iterační proces: x (0) 1 , . . . , x (0) n . . . počáteční iterace x (k+1) j = 1 aj,j bj - n i=1 i=j aj,i x (k) i , j = 1, . . . , n; k = 0, 1, . . . . (3.7) Vidíme, že vztah pro x (k+1) j odpovídá formálnímu vyřešení j-té rovnice vzhledem k proměnné xj. Algoritmus (MATLAB). jmathews('chap 3'): jacobi.m (+ demo a3 4.m) Jestliže nyní v (3.7) použijeme pro i < j již spočtenou iteraci x (k+1) i místo x (k) i , obdržíme následující zřejmě poněkud rychleji konvergující iterační proces, který je obdobou Seidelovy metody 2.7.2. Definice 3.20 (Gauss-Seidelovy iterace). Položíme N := a1,1 0 . . . 0 a2,1 a2,2 . . . 0 ... ... ... ... an,1 an,2 . . . an,n , M := N - A = 0 -a1,2 . . . -a1,n 0 0 . . . -a2,n ... ... ... ... 0 0 . . . 0 , pak z Nx(k+1) = Mx(k) + b dostáváme ihned rovnice pro NUMERICKÉ METODY 35 Gauss-Seidelův iterační proces: x (0) 1 , . . . , x (0) n . . . počáteční iterace x (k+1) j = 1 aj,j bj - j-1 i=1 aj,i x (k+1) i - n i=j+1 aj,i x (k) i , j = 1, . . . , n; k = 0, 1, . . . . (3.8) Algoritmus (MATLAB). jmathews('chap 3'): gseid.m (+ demo a3 5.m) Definice 3.21. Čtvercová matice A řádu n se nazývá striktně diagonálně dominantní, jestliže |aj,j| > n i=1 i=j |aj,i| pro j = 1, 2, . . . , n. Věta 3.22 (Konvergence Jacobiho a Gauss-Seidelovy metody). Necht' A je striktně diagonálně dominantní matice. Potom SLR Ax = b má jediné řešení a Jacobiho i Gauss-Seidelův iterační proces konverguje k tomuto řešení při libovolně zvolené počáteční aproximaci x(0) . Poznámka 3.23. Obecně z konvergence jedné z obou metod neplyne konvergence druhé. Pokud kon- vergují obě, pak Gauss-Seidelova metoda konverguje rychleji. Příklad 3.24. Jak ukazuje následující příklad, může pouhou změnou pořadí rovnic dojít ke ztrátě konvergence. Uvažme SLR 4x - y + z = 7 4x - 8y + z = -21 -2x + y + 5z = 15 Jacobiho iterace tohoto SLR odstartované v počáteční aproximaci x(0) = [1, 2, 2] konvergují k přesnému řešení x = [2, 4, 3], nebot' matice soustavy je striktně diagonálně dominantní. Jestliže zaměníme první a poslední rovnici, tak se tato vlastnost poruší a SLR -2x + y + 5z = 15 4x - 8y + z = -21 4x - y + z = 7 při téže počáteční iteraci diverguje od přesného řešení. Poznámka 3.25. Iterační metody jsou vhodné zejména pro rozsáhlé systémy lineárních rovnic s řídkou maticí soustavy (malým počtem nenulových prvků), kde u Gaussovy eliminace bývají problémy s nu- merickou stabilitou a větší výpočetní nároky vzhledem k nemožnosti účelně využít případné řídkosti matice soustavy. Snažíme se vhodnou záměnou pořadí rovnic (řádkovou permutací rozšířené matice soustavy) a změnou pořadí nezávisle proměnných (sloupcová permutace matice soustavy) soustředit největší nenulové prvky kolem hlavní diagonály s cílem dosáhnout (pokud možno) striktně diagonální domi- nantnosti. To je zpravidla snazší u řídké matice, která se stane pásovou (nenulové prvky soustředěny kolem hlavní diagonály). 36 VÍTĚZSLAV VESELÝ 4. Interpolace a aproximace funkcí pomocí polynomů (Algoritmy v MATLABu: vesely('nummet/interp'), jmathews('chap 4')) Formulace aproximační úlohy: Snaha nahradit danou funkci f(x) vhodnou funkcí z předem dané třídy funkcí , které mají jisté příznivější výpočetní vlastnosti (např. jednodušší výpočet derivace, integrálu apod.). Funkci (x) volíme tak, aby se v nějakém smyslu co nejméně lišila od dané funkce f(x). Velikost odchylky f(x) - (x) měříme zpravidla vhodnou normou na nějakém intervalu [a, b]. Nejčastěji se používá tzv. p-norma: f(x) - (x) p := b a |f(x) - (x)|p dx 1/p pro 1 p < maxx[a,b] |f(x) - (x)| pro p = . Je-li funkce f(x) zadána pouze tabulkou svých hodnot f(x0), f(x1), . . . , f(xn) na konečné diskrétní síti uzlových bodů (uzlů) x0, x1, . . . , xn, n N0, xi = xj pro i = j, pak místo integrální normy používáme vhodnou diskrétní (vektorovou) normu. Diskrétní analogie integrální p-normy měří velikost odchylky mezi vektory hodnot dané a aproximované funkce na síti uzlů x := [x0, x1, . . . , xn]T : f(x) - (x) p := n i=0|f(xi) - (xi)|p 1/p pro 1 p < maxi=0,1,...,n |f(xi) - (xi)| pro p = , kde f(x) := [f(x0), f(x1), . . . , f(xn)]T a (x) := [(x0), (x1), . . . , (xn)]T . Formulace interpolační úlohy: Interpolační úloha je speciálním případem diskrétní aproximační úlohy, kde f(x) - (x) = 0, tj. aproximační funkce (x) nabývá ve všech uzlech xi stejných hodnot jako funkce f(x): (xi) = f(xi) pro i = 0, 1, . . . , n. Pak říkáme, že funkce (x) interpoluje funkci f(x) na síti uzlů x v třídě funkcí . ˇ Je-li x I[x0, x1, . . . , xn], pak (x) nazýváme interpolovanou hodnotou funkce f v bodě x. ˇ Je-li x / I[x0, x1, . . . , xn], pak (x) nazýváme extrapolovanou hodnotou funkce f v bodě x. Nadále se budeme zabývat pouze polynomiální interpolací na síti uzlů x0, x1, . . . , xn v třídě všech polynomů (s reálnými koeficienty) := Pn := {Pn(x) | Pn(x) polynom stupně nejvýše n}. Poly- nom Qn(x) Pn, Q(x0) = f(x0), Q(x1) = f(x1), . . . , Q(xn) = f(xn) nazveme interpolačním poly- nomem funkce f(x) na síti uzlových bodů x0, x1, . . . , xn. Při řešení praktických aproximačních, resp. interpolačních úloh se často pracuje s bohatší třídou funkcí := Sn tvořenou po částech polynomickými funkcemi, tzv. splajny, stupně nejvýše n (viz Spline Toolbox v MATLABu). Báze v Pn: Každý polynom Pn(x) = p0 + p1x + + pnxn Pn, pi R, lze z hlediska sečítání a násobení skalárem ztotožnit s jeho vektorem koeficientů p := [p0, p1, . . . , pn], takže Pn tvoří vektorový prostor dimenze n + 1 nad tělesem reálných čísel. Zvolíme-li v Pn vhodný systém bázových polynomů Bn := {Bn,0, Bn,1, . . . , Bn,n}, pak každý polynom Pn Pn lze jediným způsobem vyjádřit jako lineární kombinaci bázových polynomů Bn,i. Zejména pro interpolační polynom Qn(x) =: a0 +a1x+ +anxn existují jednoznačně určené souřadnice qi tohoto polynomu v bázi Bn tak, že Qn(x) = n i=0 qiBn,i(x). Dále ukážeme, že vždy existuje jediný interpolační polynom a budeme se zabývat jeho vyjádřením v různých bázích. 4.1. Základní interpolační úloha v přirozené bázi Bn = {x0 , x1 , . . . , xn }. Polynomy xi skutečně tvoří bázi, nebot' korespondují s n + 1 lineárně nezávislými vektory [0, . . . , 0, 1 i+1 , 0, . . . , 0]T . Zřejmě ai = qi a tedy koeficienty interpolačního polynomu jsou přímo jeho souřadnicemi v této bázi. Úloha nalezení Qn(x) je pak ekvivalentní s řešením SLR: a0 + a1x0 + . . . + anxn 0 = f(x0) a0 + a1x1 + . . . + anxn 1 = f(x1) ... ... ... ... a0 + a1xn + . . . + anxn n = f(xn) V n+1a = f(x) , NUMERICKÉ METODY 37 kde V n+1 = 1 x0 . . . xn 0 1 x1 . . . xn 1 ... ... ... 1 xn . . . xn n = [x0 , x1 , . . . , xn ] je tzv. Vandermondova matice řádu n + 1. Pro determinant této matice platí známý vztah |V n+1| = i>k (xi - xk) = 0, nebot' xi = xk pro i = k. Matice V n+1 je tedy regulární a SLR má právě jedno řešení. Dokázali jsme tak následující tvrzení: Věta 4.1. Existuje právě jeden interpolační polynom Qn(x) stupně nejvýše n procházející body (x0, f(x0)), (x1, f(x1)), . . . , (xn, f(xn)). Věta 4.2 (Chyba interpolace). Necht' f Cn+1 [a, b] a Qn(x) je její interpolační polynom v n + 1 uzlových bodech xi [a, b], i = 0, 1, . . . , n. Pak pro každé x [a, b] existuje := (x) I[x0, x1, . . . , xn, x] tak, že En(x) := f(x) - Qn(x) = n+1(x) f(n+1) () (n + 1)! , (4.1a) |En(x)| max x[a,b] |n+1(x)| Mn+1 (n + 1)! , x [a, b], (4.1b) kde n+1(x) := (x - x0)(x - x1) . . . (x - xn) a |f(n+1) (x)| Mn+1 na [a, b] (f(n+1) C, takže dle věty 1.3 je f(n+1) (x) ohraničená na [a, b]). Důsledek 4.3. Necht' za předpokladů věty 4.2 je a = x0, b = xn a xi = x0 + ih, i = 0, 1, . . . , n, ekvidistantní sít' s krokem h > 0. Pak pro každé x [x0, xn] existuje := (x) [x0, xn] tak, že En(x) = En(x0 + th) = n+1(t) f(n+1) () (n + 1)! hn+1 , (4.2a) |En(x)| max t[0,n] |n+1(t)| Mn+1 (n + 1)! hn+1 , tj. f(x) = Qn(x) + O(hn+1 ) pro x [x0, xn], (4.2b) kde n+1(t) := t(t - 1) . . . (t - n), t [0, n]. Speciálně pro n = 1, 2, 3 dostáváme odhady |E1(x)| M2 8 h2 pro x [x0, x1], |E2(x)| M3 9 3 h3 pro x [x0, x2], |E3(x)| M4 24 h4 pro x [x0, x3]. (4.2c) Poznámka 4.4 (Rungeho jev). Otázka: Platí maxx[a,b]|En(x)| 0 pro n ? Odpověd': a) ano pro funkci f, jejíž všechny derivace jsou ohraničeny společnou konstantou M (Mn M pro n = 1, 2, . . . ). Pak totiž maxx[a,b]|n+1(x)| roste pomaleji, než (n + 1)! pro n a tudíž celý výraz v (4.1b) konverguje k nule. Tato podmínka je splněna například pro funkce sin(x) a ex . b) obecně nikoliv: například v případě funkce f(x) = 1 1+12x2 je maxx[-1,1]|En(x)| rostoucí nade všechny meze (tzv. Rungeho jev), nebot' Mn zesiluje maxima oscilující funkce |n+1(x)| tak, že se zvětšujícím se počtem uzlů roste amplituda oscilací rychleji, než (n + 1)!. Rungeho efekt je pro tuto funkci výrazný již na 11-ti bodové ekvidistantní síti (n = 10). V MATLABu koeficienty a interpolačního polynomu pro y = f(x) snadno spočteme pomocí funkce a=polyfit(x,y,n), která je vrací v obráceném pořadí. Nakonec zvolíme jemnou sít' X=linspace(-1,1).' pro kreslení grafů. Rungeho efekt pak znázorníme vykreslením průběhů Y = f(X) a P = Qn(X) (P vypočteme pomocí P=polyval(a,X)) do jednoho grafu příkazem plot(X,[Y,P]). Pro dané n následující věta udává optimální (neekvidistantní) volbu uzlů, pro niž maxima oscilací |n+1(x)| jsou minimalizována a je tak odstraněn Rungeho efekt. 38 VÍTĚZSLAV VESELÝ Věta 4.5 (Čebyševovy uzly). Položme xk = tk b-a 2 + a+b 2 , kde tk = cos (2n + 1 - 2k)/(2n + 2) pro k = 0, 1, . . . , n. Potom pro n+1(x) = (x - x0)(x - x1) . . . (x - xn) platí maxx[a,b]|n+1(x)| 2 b-a 4 n+1 , přičemž pro každou jinou sít' n+1 uzlů x0, x1, . . . , xn je maxx[a,b]|n+1(x)| maxx[a,b]|n+1(x)|. Jestliže f Cn+1 [a, b], pak interpolační polynom v uzlech xk (tzv. Čebyševovy uzly) dává nejmenší možný odhad chyby |En(x)| 2 b - a 4 n+1 maxx[a,b]|f(n+1) (x)| (n + 1)! , x [a, b]. (4.2d) Poznámka 4.6. Použijeme-li pro funkci f(x) = 1 12+x2 z poznámky 4.4 Čebyševovy uzly, pak maxx[a,b]|En(x)| 0 pro n a Rungeho efekt je tak odstraněn. Toto platí obecně pro každou funkci i při snížených požadavcích na její hladkost. Dá se totiž ukázat, že pro každou funkci f C1 [a, b] její interpolační polynom Qn(x) v Čebyševových uzlech vždy konverguje, a to dokonce stejnoměrně, k funkci f na intervalu [a, b] a je tedy odstraněn Rungeho efekt. 4.2. Lagrangeův interpolační polynom (Bn = {Ln,0(x), Ln,1(x), . . . , Ln,n(x)}). Ln,i(x) := (x - x0) . . . (x - xi-1)(x - xi+1)(x - xn) (xi - x0) . . . (xi - xi-1)(xi - xi+1)(xi - xn) = n j=0 j=i (x - xj) n j=0 j=i (xi - xj) = = n+1(x) (x - xi)n+1(xi) . . . tzv. Lagrangeův interpolant. (4.3a) Zřejmě Ln,i(xj) = i,j (Kroneckerův symbol), takže ihned snadno ověříme lineární nezávislost poly- nomů Ln,i: a(x) := n j=0 ajLn,j(x) 0 0 = a(xi) = ai pro i = 0, 1, . . . , n. Tedy Ln,0, Ln,1, . . . , Ln,n skutečně tvoří tzv. Lagrangeovu bázi v Pn. V této bázi dostáváme následující vyjádření interpolačního polynomu: Qn(x) = n i=0 f(xi)Ln,i(x) = n i=0 f(xi) n+1(x) (x - xi)n+1(xi) (4.3b) Tedy qi = f(xi) a funkční hodnoty f(xi) tak jsou přímo souřadnicemi interpolačního polynomu v Lagrangeově bázi. Poznamenejme ještě, že v případě n = 0 klademe n j=0 j=i (x - xj) = 1 ve (4.3a), a tedy L0,0 1. Obecně Ln,i(xi) = 1 není maximální hodnotou polynomu Ln,i(x) na intervalu I[x0, x1, . . . , xn]. Věta 4.7 (Lagrangeův interpolační polynom pro ekvidistantní uzly). Je-li xi = x0 +ih, h > 0, i = 0, 1, . . . , n, sít' ekvidistantních uzlů, pak (4.3a) a (4.3b) nabudou po řadě tvaru: Ln,i(x0 + th) = (-1)n-i i! (n - i)! n j=0 j=i (t - j) = (-1)n-i n! n i n+1(t) t - i (4.4a) a Qn(x0 + th) = n i=0 f(xi)Ln,i(x0 + th) = n+1(t) n! n i=0 (-1)n-i n i f(xi) t - i . (4.4b) Algoritmus (MATLAB). vesely('nummet/interp'): lagr.m jmathews('chap 4'): lagran.m (+ demo a4 3.m) NUMERICKÉ METODY 39 4.3. Newtonův interpolační polynom (Bn = {0(x), 1(x), . . . , n(x)}). Nevýhodou vyjádření interpolačního polynomu v přirozené či Lagrangeově bázi je nutnost sestrojit interpolační polynom zcela znovu, kdykoliv potřebujeme ke stávající síti přidat další uzel. Tento problém odstraňuje použití Newtonovy báze Bn = {0(x), 1(x), . . . , n(x)}, kde 0(x) 1 a pro i > 0 je i(x) = (x - x0)(x - x1) . . . (x - xi-1) polynom stupně i jako ve větě 4.2. Snadno nahlédneme, že tento systém je lineárně nezávislý. Totiž matice, jejímiž sloupci jsou po řadě vektory koeficientů i(x) pro i = 0, 1, . . . , n, tvoří čtvercovou horní trojúhelníkovou matici řádu n + 1 s jed- notkovou diagonálou (i(x) má jednotkový koeficient u nejvyšší mocniny xi ). Tato matice je regulární s determinantem rovným jedné, takže její sloupce musí být lineárně nezávislé. Věta 4.8. Necht' Qk(x) = q0+q11(x)+ +qkk(x) je vyjádření interpolačního polynomu funkce f(x) v uzlových bodech x0, x1, . . . , xk (k = 0, 1, . . . , n). Pro tento tzv. Newtonův tvar interpolačního polynomu platí následující rekurentní vztah pro výpočet souřadnic qi: Q0(x) = q0 = f(x0) Qk(x) = Qk-1(x) + qkk(x) pro k = 1, . . . , n, kde qk = f(xk) - Qk-1(xk) k(xk) = f(xk) - Qk-1(xk) k-1 i=0 (xk - xi) . (4.5) Poznamenejme, že pro vyhodnocení Newtonova interpolačního polynomu Qn(x) v proměnné x = z lze použít zobecnění Hornerova schématu z odstavce 1.2: (. . . (( n-1 qn(z - xn-1) + qn-1)(z - xn-2) + qn-2)(z - xn-3) + + q1)(z - x0) + q0. Definice 4.9. Koeficient qk v (4.5) se nazývá diferenční kvocient k-ho řádu a značí se qk =: f[x0, x1, . . . , xk] (alternativně též qk =: f0,1,...,k nebo qk =: [x0, x1, . . . , xk]f ). Newtonův interpolační polynom pak můžeme psát ve tvaru: Qn(x) = f[x0] + f[x0, x1](x - x0) + + f[x0, x1, . . . , xn](x - x0)(x - x1) . . . (x - xn-1), (4.6a) kde f[x0] = f(x0). Věta 4.10. Diferenční kvocienty lze vyjádřit v explicitním tvaru qk = f[x0, x1, . . . , xk] = k i=0 f(xi) k j=0 j=i (xi - xj) = k i=0 f(xi) k+1(xi) . (4.6b) Důsledek 4.11. V případě ekvidistantní sítě xi = x0 +ih, h > 0, i = 0, 1, . . . , k, nabude (4.6b) tvaru qk = f[x0, x1, . . . , xk] = 1 k! hk k i=0 (-1)k-i k i f(xi). (4.6c) Důkaz. V důkazu věty 4.7 jsme ukázali k j=0 j=i (xi - xj) = hk i! (k - i)! (-1)k-i = hk k! (-1)k-i k i . Po dosazení do (4.6b) dostaneme ihned (4.6c). Důsledek 4.12. Je-li p libovolná permutace indexů {0, 1, . . . , k}, pak qk = f[x0, x1, . . . , xk] = f[xp(0), xp(1), . . . , xp(k)]. Hodnota diferenčního kvocientu f[x0, x1, . . . , xk] tedy nezávisí na pořadí uzlových bodů x0, x1, . . . , xk. Důkaz. Tvrzení je bezprostředním důsledkem (4.6b), nebot' sumace ani součin zde nezávisí na pořadí. Věta 4.13 (Rekurentní vztah pro diferenční kvocienty). f[x0, x1, . . . , xk] = f[x1, x2, . . . , xk] - f[x0, x1, . . . , xk-1] xk - x0 (4.7a) 40 VÍTĚZSLAV VESELÝ Důsledek 4.14. V případě ekvidistantní sítě xi = x0 + ih, h > 0, i = 0, 1, . . . , k, nabude (4.7a) tvaru f[x0, x1, . . . , xk] = k f(x0) k! hk , (4.7b) kde pro k-tou diferenci k f(x0) (viz definici 2.4) platí vztah k f(x0) = k i=0 (-1)k-i k i f(xi). (4.7c) Důsledek 4.15 (1. Newtonova interpolační formule pro ekvidistantní uzly). Je-li xi = x0 + ih, h > 0, i = 0, 1, . . . , n, ekvidistantní sít' uzlů a položíme-li x = x0 + th, t [0, n], pak (4.6a) lze psát ve tvaru Qn(x0 + th) = f(x0) + t 1! f(x0) + t(t - 1) 2! 2 f(x0) + + t(t - 1) . . . (t - (n - 1)) n! n f(x0) = (4.8a) = n k=0 k(t) k! k f(x0). Důsledek 4.16 (2. Newtonova interpolační formule pro ekvidistantní uzly). Je-li xi = x0 + ih, h > 0, i = 0, 1, . . . , n, ekvidistantní sít' uzlů a položíme-li x = xn + th, t [-n, 0], pak (4.6a) lze psát ve tvaru Qn(xn + th) = f(xn) + t 1! f(xn-1) + t(t + 1) 2! 2 f(xn-2) + + t(t + 1) . . . (t + n - 1) n! n f(x0). (4.8b) Poznámka 4.17. Pokud x0 < x1 < < xn, pak vzorec (4.8a) je z hlediska numerické stability vhodný pro interpolaci na začátku intervalu, kde t = x-x0 h nabývá malých hodnot, zatímco (4.8b) se naopak hodí pro interpolaci na konci intervalu. Algoritmus 4.18 (Diferenční schéma). xi f[] f[] f[ ] f[ ] . . . x0 f0 x1 f1 f1-f0 x1-x0 = f0,1 x2 f2 f2-f1 x2-x1 = f1,2 f1,2-f0,1 x2-x0 = f0,1,2 x3 f3 f3-f2 x3-x2 = f2,3 f2,3-f1,2 x3-x1 = f1,2,3 f1,2,3-f0,1,2 x3-x0 = f0,1,2,3 ... ... ... ... ... ... V ekvidistantním případě stačí dle (4.7b) počítat místo diferenčních kvocientů pouze diference: xi fi fi 2 fi 3 fi . . . x0 f0 x1 f1 f1 - f0 = f0 x2 f2 f2 - f1 = f1 f1 - f0 = 2 f0 x3 f3 f3 - f2 = f2 f2 - f1 = 2 f1 2 f1 - 2 f0 = 3 f0 ... ... ... ... ... ... Algoritmus (MATLAB). jmathews('chap 4'): newpoly.m (+ demo a4 5.m) Věta 4.19 (Jiné vyjádření chyby interpolace). Necht' Qn(x) je interpolační polynom funkce f(x) v uzlových bodech x0, x1, . . . , xn, pak pro libovolné x / {x0, x1, . . . , xn} lze chybu interpolace vyjádřit ve tvaru En(x) = f(x) - Qn(x) = n+1(x)f[x0, x1, . . . , xn, x]. (4.9) Poznámka 4.20. Porovnejme za dodatečného předpokladu f Cn+1 [a, b] oba vztahy pro chybu (4.1a) a (4.9). Vidíme, že musí existovat = (x) I[x0, x1, . . . , xn, x] takové, že f[x0, x1, . . . , xn, x] = f(n+1) ((x)) (n + 1)! . NUMERICKÉ METODY 41 Dá se dokonce ukázat, že pro f Ck [a, b], lze f[x0, x1, . . . , xk] spojitě prodloužit na funkci k + 1 proměnných x0, x1, . . . , xk [a, b] (viz následující dvě tvrzení), kdy již nemusíme předpokládat, že xi jsou navzájem různé. Věta 4.21. Necht' f Ck [a, b], m f(k) (x) M pro každé x I[x0, x1, . . . , xk], kde x0, x1, . . . , xk [a, b] jsou libovolné. Pak platí f[x0, x1, . . . , xk] je spojitá funkce k + 1 proměnných na [a, b]k+1 , jejíž hodnota nezávisí na pořadí x0, x1, . . . , xk ve smyslu důsledku 4.12 (4.10) m k! f[x0, x1, . . . , xk] M k! (4.11) f[x0, x1, . . . , xk] = f(k) () k! pro vhodné I[x0, x1, . . . , xk] (4.12) f[x, x, . . . , x (k+1)× ] = f(k) (x) k! (4.13) Důsledek 4.22. Vztah (4.9) platí pro každé x [a, b]. Poznámka 4.23. Jestliže totožné uzly uspořádáme vedle sebe, lze pro výpočet f[x0, x1, . . . , xk] mod- ifikovat (4.7a) a výsledné diferenční schéma 4.18 tak, že za výrazy typu f[xi, xi, . . . , xi (j+1)× ] dosazujeme f(j) (xi) j! podle (4.13). Tento algoritmus je základem tzv. Hermitovy interpolace, kde v (j + 1)-násobném uzlu xi předepisujeme kromě funkční hodnoty f(xi) ještě dalších j derivací f (xi), f (xi), . . . , f(j) (xi). Poznámka 4.24. Pro funkci f Cn+1 [a, b] lze podle důsledku 4.22 a s uvážením tvrzení (4.10) a (4.12) vybrat (x) I[x0, x1, . . . , xn, x] tak, že (n + 1)! f[x0, x1, . . . , xn, x] = f(n+1) ((x)) je spojitá funkce v proměnné x na [a, b] při pevně zvolených uzlových bodech x0, x1, . . . , xn. Následující věta říká, že zvýšením hladkosti funkce f(x) o r řádů se obdobně zvýší hladkost funkce f(n+1) ((x)) až do řádu r. Věta 4.25. Necht' f Cn+r+1 [a, b] (r N0) je dána v uzlech x0, x1, . . . , xn [a, b] a (x) I[x0, x1, . . . , xn, x] je prvek z chyby interpolace (4.1a) vybraný podle poznámky 4.24 stejně jako v (4.12). Potom f(n+1) ((x)) Cr [a, b] a platí dr dxr f(n+1) ((x)) (n + 1)! = r! f(n+r+1) (r(x)) (n + r + 1)! , (4.14) kde r := r(x) I[x0, x1, . . . , xn, x], 0(x) := (x). 42 VÍTĚZSLAV VESELÝ 5. Numerické derivování (Algoritmy v MATLABu: jmathews('chap 6')) Princip numerického derivování: Numerické derivování používáme pro přibližný výpočet prvé nebo vyšší derivace funkce f(x) zadané tabulkou hodnot v uzlových bodech x0, x1, . . . , xn. Této metody můžeme s výhodou použít i v případě, kdy výpočet derivací f(x) je obtížný pro složité analytické vyjádření. Nejčastěji klademe f (x) (x), kde (x) je vhodná funkce aproximující dobře f(x) na zvoleném intervalu [a, b]. Nejjednodušší formule pro numerické derivování funkce f(x) lze získat derivováním jejího interpolačního polynomu: f (x) Qn(x). V případě derivací vyšších řádů postupujeme obdobně. V každém případě je nutno počítat s tím, že numerické derivování je méně přesné než pouhá interpolace, přičemž chyby vzrůstají progresivně s řádem derivace. Je to dáno především tím, že z blízkosti funkčních hodnot v sousedních uzlových bodech nemusí obecně plynout blízkost derivací, pokud funkce není dostatečně hladká. Vysoká hladkost funkce bude tedy nutným požadavkem pro získání hodnot derivací s vyhovující přesností. Necht' Qn(x) interpoluje funkci f(x)v uzlových bodech x0, x1, . . . , xn. Pak z vyjádření f(x) = Qn(x) + En(x) dostáváme f (x) = Qn(x) + En(x) f(r) (x) = Q(r) n (x) + E(r) n (x), r N. (5.1) Zatímco vyjádření derivace Q (r) n = n i=0 f(xi)L (r) n,i(x) dostáváme přímo z Lagrangeova tvaru (4.3b), formální vyjádření chyby derivace jako r-té derivace chybového členu En(x) z (4.1a) je podstatně pracnější a vyžaduje užití věty 4.25. Omezme se nejprve na případ 1. derivace. Věta 5.1 (Chyba 1. derivace). Necht' f Cn+2 [a, b]) a Qn(x) je její interpolační polynom v n + 1 uzlových bodech x0, x1, . . . , xn [a, b]. Pak pro každé x [a, b] platí En(x) = E(0)(x) + E(1)(x), kde (5.2a) E(0)(x) = n+1(x) f(n+1) (0(x)) (n + 1)! , E(1)(x) = n+1(x) f(n+2) (1(x)) (n + 2)! , kde 0, 1 I[x0, x1, . . . , xn, x] a E(1)(x) = 0 pro x {x0, x1, . . . , xn}, kdy stačí předpokládat pouze f Cn+1 [a, b]) . (5.2b) Podobně pro r-tou derivaci dostaneme po úpravě. Věta 5.2 (Chyba r-té derivace). Necht' f Cn+r+1 [a, b]) , r N, a Qn(x) je její interpolační polynom v n + 1 uzlových bodech x0, x1, . . . , xn [a, b]. Pak pro každé x [a, b] platí E(r) n (x) = E(0)(x) + E(1)(x) + + E(r)(x), kde (5.3a) E(j)(x) = r! (r - j)! (r-j) n+1 (x) f(n+j+1) (j(x)) (n + j + 1)! , j I[x0, x1, . . . , xn, x], j = 0, 1, . . . , r, a E(r)(x) = 0 pro x {x0, x1, . . . , xn}, kdy stačí předpokládat pouze f Cn+r [a, b]) . (5.3b) V případě ekvidistantní sítě xi = x0 + ih, h > 0, t = x-x0 h , je E (r) n (x) = O(hn+1-r ) a vidíme jak přesnost klesá s řádem derivace. Vztah (5.3b) lze v tomto případě psát ve tvaru E(j)(t) = r!hn+1-r (r-j) n+1 (t) (r - j)! f(n+j+1) (j(x)) (n + j + 1)! , j I[x0, x1, . . . , xn, x], j = 0, 1, . . . , r, a E(r)(t) = 0 pro t {0, 1, . . . , n}, kdy stačí předpokládat pouze f Cn+r [a, b]) . (5.3c) Užijeme-li pro Qn(x) Lagrangeova tvaru interpolačního polynomu pro obecné uzly (viz (4.3a) a (4.3b)), resp. pro ekvidistantní uzly (viz (4.4a) a (4.4b)), obdržíme následující výsledné tvary formulí pro výpočet r-té derivace (r 1) funkce f. NUMERICKÉ METODY 43 Věta 5.3. Necht' x0, x1, . . . , xn [a, b], n 0, jsou uzlové body funkce f a x [a, b] je libo- volné. Jestliže označíme r = r - 1 pro x {x0, x1, . . . , xn} r pro x / {x0, x1, . . . , xn} , kde r 1 udává řád derivace a f Cn+r +1 [a, b], pak platí a) pro obecné uzly xi: f(r) (x) = n i=0 f(xi)L (r) n,i(x) + r! r j=0 (r-j) n+1 (x) (r - j)! f(n+j+1) (j) (n + j + 1)! , (5.4a) j I[x0, x1, . . . , xn, x]. b) pro ekvidistantní uzly xi = x0 + ih, h > 0: f(r) (x) = 1 hr n i=0 f(xi)D (r) i (t) + r! hn+1-r r j=0 (r-j) n+1 (t) (r - j)! f(n+j+1) (j) (n + j + 1)! , (5.4b) j I[x0, x1, . . . , xn, x], kde Di(t) = Ln,i(x0 + th) a tedy D (r) i (t) = (-1)n-i i! (n - i)! dr dtr n j=0 j=i (t - j) , t = x - x0 h , tj. zejména x = xk t = k. (5.4c) Důsledek 5.4. Speciálně pro r = 1 a f Cn+1 [a, b] pro x {x0, x1, . . . , xn} Cn+2 [a, b] pro x / {x0, x1, . . . , xn} dostáváme a) pro obecné uzly xi: f (x) = n i=0 f(xi)Ln,i(x) + n+1(x) f(n+1) (0) (n + 1)! + n+1(x) f(n+2) (1) (n + 2)! E(1)(x) , (5.5a) 0, 1 I[x0, x1, . . . , xn, x] a E(1)(x) = 0 pro x {x0, . . . , xn}. b) pro ekvidistantní uzly xi = x0 + ih, h > 0: f (x) = 1 h n i=0 f(xi)Di(t) + hn n+1(t) f(n+1) (0) (n + 1)! + n+1(t) f(n+2) (1) (n + 2)! E(1)(t) , (5.5b) 0, 1 I[x0, x1, . . . , xn, x] a E(1)(t) = 0 pro t {0, . . . , n}, kde Di(t) = (-1)n-i i! (n - i)! n j=0 j=i (t - j) , t = x - x0 h , tj. zejména x = xk t = k. (5.5c) Poznámka 5.5 (Speciální případy). V praxi se často setkáváme s ekvidistantní sítí uzlů a úlohou počítat derivace v těchto uzlech. Z hlediska přesnosti je ovšem výhodnější volba Čebyševových uzlů (Věta 4.5), nebot' výskyt Rungeho jevu může zcela znehodnotit odhady derivací v důsledku oscilačního průběhu Qn(x). Pro ekvidistantní uzly můžeme ze vztahů (5.5b) a (5.5c) odvodit následující nejčastěji používané formule pro 1. derivaci (značíme fi = f(xi), fi = f (xi)): (1) 2-bodová formule (n = 1, f C2 [a, b], E1 = O(h)) D0(t) = (-1)1 0! 1! (t - 1) = -1 D1(t) = (-1)0 1! 0! (t - 0) = 1 2(t) = t(t - 1) = 2t - 1 2(0) = -1, 2(1) = 1. Odtud 44 VÍTĚZSLAV VESELÝ f0 = 1 h (f1 - f0) - h 2 f (0), f1 = 1 h (f1 - f0) + h 2 f (1), 0, 1 I[x0, x1]. (5.6) (2) 3-bodová formule (n = 2, f C3 [a, b], E2 = O(h2 )) D0(t) = (-1)2 0! 2! (t - 1)(t - 2) = 1 2 [t2 - 3t + 2] = 1 2 (2t - 3). D1(t) = (-1)1 1! 1! t(t - 2) = -2(t - 1). D2(t) = (-1)0 2! 0! = t(t - 1) = 1 2 (2t - 1). 3(t) = t(t - 1)(t - 2) = t(t2 - 3t + 2) = 3t2 - 6t + 2 3(0) = 2, 3(1) = -1, 3(2) = 2. Odtud f0 = 1 h - 3 2 f0 + 2f1 - 1 2 f2 + h2 3 f (0), f1 = 1 h - 1 2 f0 + 1 2 f2 - h2 6 f (1), f2 = 1 h 1 2 f0 - 2f1 + 3 2 f2 + h2 3 f (2), 0, 1, 2 I[x0, x1, x2]. (5.7) Podobným způsobem lze konstruovat (n + 1)-bodové formule pro libovolné n > 2. Jsou-li dány hodnoty (xi, fi) pro i = 0, 1, . . . , n, kde n je velké, není nutné (a ani vhodné) používat (n+1)-bodovou formuli pro výpočet f0, f1, . . . , fn. Stačí totiž provádět postupně pouze lokální interpolaci dat v okolí xi pro i = 0, 1, . . . , n polynomem Qm nižšího stupně m n, takže například pro m = 1 můžeme (5.6) přepsat do přeindexovaného tvaru: fi = 1 h (fi+1 - fi) - h 2 f (i), i [xi, xi+1]. (5.6.1) fi = 1 h (fi - fi-1) + h 2 f (i), i [xi-1, xi]. (5.6.2) Podobně pro m = 2 přepíšeme (5.7) do tvaru: fi = 1 2h (-3fi + 4fi+1 - fi+2) + h2 3 f (i), i [xi, xi+2] (5.7.1) fi = 1 2h (fi+1 - fi-1) - h2 6 f (i), i [xi-1, xi+1] (5.7.2) fi = 1 2h (fi-2 - 4fi-1 + 3fi) + h2 3 f (i), i [xi-2, xi] (5.7.3) Způsob použití formulí: (5.7.1) pro i = 0 . . . tj. na začátku intervalu (tzv. levostranná formule) (5.7.2) pro 0 < i < n . . . tj. uvnitř intervalu (tzv. centrální formule) (5.7.3) pro i = n . . . tj. na konci intervalu (tzv. pravostranná formule) Analogicky postupujeme v případě formulí s lichým počtem bodů (tj. m = 2q > 2 sudé), kdy dostaneme jednu centrální formuli, q levostranných a q pravostranných formulí. Všimněte si, že velikost chyby centrální formule (5.7.2) je poloviční ve srovnání s (5.7.1) a (5.7.3). Je logické, že centrální formule dávají nejmenší chybu, nebot' využívají stejný počet hodnot z obou stran k xi, které jsou nejblíž k fi a tudíž nejpřesněji popisují chování f(x) v okolí xi. Všimněte si dále rovněž symetričnosti formulí (5.7.1) a (5.7.3). Formule (5.7.3) se dostane z (5.7.1) inverzním uspořádáním uzlů: f(xi-j) = f(xi + j(-h)) a záměnou h za -h, která odpovídá obrácení znaménka derivace při záměně nezávisle proměnné t za -t. NUMERICKÉ METODY 45 Podobné vlastnosti vykazují i formule s vyšším lichým počtem bodů (5-bodové, 7-bodové atd.), které se nejčastěji používají vzhledem k menší chybě centrální formule. Necentrální formule se aplikují pouze na okrajích intervalu. Poznámka 5.6 (Problémy s přesností u formulí pro numerické derivování). Jestliže jsou hodnoty fi dány s chybou fi = ~fi + Ei, |Ei| i, může tato okolnost podstatně ovlivnit kvalitu výsledného odhadu fi . Ukážeme to na příkladě formule (5.7.2). Jelikož chyba i a chyba formule se sčítají (viz 1.33), lze celkovou chybu T shora odhadnout takto: |T| = |fi - 1 2h ( ~fi+1 - ~fi-1)| = |fi - 1 2h (fi+1 - fi-1) + 1 2h (Ei+1 - Ei-1)| h2 6 |f (i)| + 1 2h (i+1 + i-1). Položíme-li = max(i-1, i+1) a M3 = maxx[xi-1,xi+1]|f (x)|, pak |T| h + h2 6 M3 První člen chyby h (např. může být chyba způsobená zaokrouhlováním, měřením apod.) závisí nepřímo úměrně na h, zatímco druhý člen (chyba metody) závisí přímo úměrně na h. Vzniká otázka, jak optimálně volit h, aby celková chyba byla minimální. 1. Strategie: hledáme minimum funkce g(h) = h + h2 6 M3. Ze vztahu g (h) = - h2 + h 3 M3 = 0 dostáváme hopt = 3 3 M3 . 2. Strategie: hledáme h tak, aby příspěvek obou členů chyby byl stejný. Ze vztahu h = h2 6 M3 dostáváme hopt = 3 6 M3 . V případě, že hodnoty fi jsou dány s malou přesností (např. byly získány empiricky), není vhodné použít formule pro numerické derivování přímo, nebot' hrozí riziko výrazného zkreslení výsledků. V takových případech je lépe nejdříve naměřené hodnoty "vyhladit" například metodou nejmenších čtverců a pak teprve použít formule pro numerické derivování. Algoritmus (MATLAB). jmathews('chap 6'): difflim.m (+ demo a6 1.m) V tomto algoritmu je snaha nalézt "nejlepší" aproximaci f (x) zvolením optimálního kroku h pro formuli (5.7.2) jeho postupným zjemňováním metodou půlení: f (x) Dk(x) := f x + h 2k - f x - h 2k h/2k-1 pro k = 0, 1, . . . až do dosažení ukončující podmínky: |Dk+1 - Dk| |Dk - Dk-1| . . . odhad chyby začíná vzrůstat, nebo |Dk - Dk-1| < Tolerance . . . odhad chyby dosáhl požadované přesnosti. 46 VÍTĚZSLAV VESELÝ 6. Numerické integrování (Algoritmy v MATLABu: jmathews('chap 7')) Definice 6.1. Necht' - a x0 < x1 < < xn b + a W(x)f(x) je funkce inte- grovatelná na intervalu [a, b], kde W(x) je předem daná, tzv. váhová funkce, a f(x) je funkce zadaná v uzlových bodech x0, x1, . . . , xn tabulkou funkčních hodnot (případně i derivacemi některých řádů). Kvadraturní formulí nazýváme formuli pro přibližný výpočet integrálu b a W(x)f(x) dx, závisející pouze na W(x), uzlech x0, x1, . . . , xn a zadaných tabulkových hodnotách. Obecný tvar kvadraturní formule závisející lineárně na tabulkových hodnotách f(xi), i = 0, 1, . . . , n: b a W(x)f(x) dx = n i=0 Aif(xi) + Ef , (6.1) kde je Ef . . . chyba kvadraturní formule, a pro i = 0, 1, . . . , n jsou xi . . . kvadraturní uzly Ai . . . váhové koeficienty (váhy) závislé pouze na W(x) a xi parametry formule. Definice 6.2. Kvadraturní formule (6.1) se nazývá ˇ uzavřená, jestliže a = x0 a xn = b, ˇ otevřená, jestliže a < x0 a xn < b, tj. formule nezávisí na hodnotách funkce v krajních bodech intervalu. ˇ polouzavřená, resp. polootevřená, jestliže (1) a = x0 a xn < b, tj. je uzavřená zleva a otevřená zprava, nebo (2) a < x0 a xn = b, tj. je otevřená zleva a uzavřená zprava. Definice 6.3. Celé číslo 0 nazveme stupněm přesnosti kvadraturní formule (6.1), jestliže EP = 0 pro všechny polynomy P(x) stupně nejvýše (tj. P P) a EP = 0 pro některý polynom P(x) stupně + 1. Poznámka 6.4. Zřejmě pro libovolná , R a funkce f(x), g(x) splňující předpoklady 6.1 je Ef+g = Ef + Eg, nebot' dle (6.1) je Ef = b a W(x)f(x) dx - n i=0 Aif(xi), kde výraz vpravo závisí lineárně na f. Odtud ihned vidíme, že v případě formule (6.1) se stupněm přesnosti je Ep = 0 pro každý polynom P(x) =: a0 + a1x + . . . ax + a+1x+1 stupně + 1 (a+1 = 0), nebot' EP = a0 E1 =0 +a1 Ex =0 + + a Ex =0 + a+1 =0 Ex+1 = 0 Ex+1 = 0. Protože P(+1) () ( + 1)! a+1, dostáváme odtud ihned vyjádření pro chybu ve tvaru: EP = Ex+1 ( + 1)! =:e P(+1) () =: eP(+1) (), R. Dá se tedy očekávat, že Ef bude mít výše uvedený tvar pro každou funkci f C+1 [a, b] a vhodné [a, b] (viz dále). V dalším se omezíme na ohraničený interval [a, b], tj. - < a < b < +. Věta 6.5. Je-li W(x) integrovatelná na [a, b] a f Cn+1 [a, b], potom pro formuli (6.1) se stupněm přesnosti n platí Ai = b a W(x)Ln,i(x) dx pro i = 0, 1, . . . , n, (6.1a) Ef = 1 (n + 1)! b a W(x)n+1(x)f(n+1) ((x)) dx, (x) [a, b]. (6.1b) NUMERICKÉ METODY 47 Poznámka 6.6. (1) Pro n je tedy (6.1) jednoznačně určena vztahy (6.1a) a (6.1b) a dostane se formálním zintegrováním Lagrangeova interpolačního vztahu (4.3b) f(x) = n i=0 Ln,i(x)f(xi) + n+1(x) f(n+1) ((x)) (n + 1)! En(x) , s využitím vyjádření (4.1a) pro chybu En(x). (2) Metoda neurčitých koeficientů pro určení Ai při daných uzlech xi, i = 0, 1, . . . , n: Necht' 0(x), 1(x), . . . , n(x) je vhodná polynomická báze v Pn (například j(x) = xj ). Z podmínky, aby (6.1) byla přesná pro všechny polynomy j(x), 0 j n, dostáváme systém n + 1 lineárních rovnic pro n + 1 neznámých Ai: n i=0 Aij(xi) = b a W(x)j(x) dx pro j = 0, 1, . . . , n, (6.1c) Protože j(x) jsou bázové polynomy stupně nejvýše n bude vzhledem k linearitě chyby 6.4 formule přesná nejen pro ně, ale pro všechny polynomy stupně nejvýše n a bude tedy formulí dosahující přesnosti alespoň n. V případě ekvidistantní sítě obdržíme tzv. Newton-Cotesovy formule, které budou dále podrobně popsány v odstavci 6.1. (3) Ve (2) byly uzly xi předepsány a volnými parametry byly pouze koeficienty Ai, i = 0, 1, . . . , n. Všech parametrů je tedy celkem 2n + 2. Zvolíme-li obecně k volných parametrů (1 k 2n + 2), je snaha sestavit k rovnic typu (6.1c). Pokud tento systém je řešitelný, obdržíme analogicky formuli přesnosti alespoň k - 1. Poznamenejme ovšem, že tento systém je lineární pouze v koeficientech Ai, nikoliv ale v uzlech xi, což komplikuje nalezení řešení. Ponecháme-li všechny parametry volné (k = 2n + 2), pak lze pro několik typů váhových funkcí W(x) tímto postupem zkonstruovat tzv. Gaussovy kvadraturní formule, které dosahují maximálního stupně přesnosti = 2n + 1. (4) Obecný postup konstrukce kvadraturních formulí lze vzhledem ke (3) formulovat takto: ˇ Zvolíme vhodnou váhovou funkci W(x) a zadáme r (0 r < 2n+2) vhodných omezujících podmínek pro parametry (například ve (2) bylo předepsáno r = n + 1 uzlů). ˇ Sestrojíme 2n + 2 rovnic tvořených r předepsanými podmínkami doplněnými o k = 2n + 2 - r podmínek typu (6.1c). ˇ Vyřešením tohoto systému nalezneme všech 2n + 2 parametrů (váhy + uzly) a obdržíme formuli se stupněm přesnosti alespoň k - 1 ( 2n + 1 - r). 6.1. Newton-Cotesovy kvadraturní formule. Definice 6.7. (1) Řekneme, že uzly x0, x1, . . . , xn jsou symetrické na intervalu [a, b], jestliže jsou rozmístěny symetricky vzhledem ke středu s = a+b 2 intervalu, tj. s-xi = xn-i -s =: hi pro i = 0, 1, . . . , n. (2) Je-li n celé číslo, pak par(n) := 1 pro n sudé -1 pro n liché udává jeho paritu. (3) Je-li funkce g(x) definována na [a, b], pak řekneme, že je sudá, resp. lichá na [a, b], jestliže g(s + t), t [-(b - a)/2, (b - a)/2], je sudá, resp. lichá; neboli platí g(s - t) = g(s + t), resp. g(s - t) = -g(s + t). V takovém případě píšeme para,b(g) = 1, resp. = -1. Pokud g není ani sudá, ani lichá, klademe para,b(g) = 0. Lemma 6.8. Necht' x0, x1, . . . , xn jsou symetrické uzly na [a, b]. Potom para,b(n+1) = par(n+1) = -par(n), neboli n+1(x) je sudá, resp. lichá funkce na [a, b] v závislosti na tom, zda počet uzlů je sudý, resp. lichý. Věta 6.9. Necht' x0, x1, . . . , xn jsou symetrické uzly na [a, b] a W(x) je integrovatelná sudá nebo lichá funkce na [a, b], potom pro koeficienty Ai určené vztahem (6.1a) platí symetrie An-i = para,b(W) Ai pro i = 0, 1, . . . , n. (6.2) 48 VÍTĚZSLAV VESELÝ Věta 6.10. Necht' x0, x1, . . . , xn jsou symetrické uzly na [a, b], W(x) spojitá na [a, b], para,b(W) = -par(n + 1) a f Cn+2 [a, b]. Potom chybu formule určenou vztahem (6.1b) lze vyjádřit ve tvaru Ef = - b a q(x) f(n+1) ((x)) (n + 1)! dx = - b a q(x) f(n+2) ((x)) (n + 2)! dx, (6.3) kde q(x) := x a W(u)n+1(u) du a (x), (x) [a, b], přičemž f(n+2) ((x)) C[a, b]. V dalším se již pro jednoduchost omezíme pouze na případ jednotkové váhové funkce W(x) 1. Lemma 6.11. Necht' a = x-1 < x0 < < xn+1 = b a xi = x0 + ih, h > 0 (i = 0, 1, . . . , n) jsou ekvidistantní uzly na intervalu [a, b]. Označíme-li Ij := xj+1 xj n+1(x) dx pro j = -1, 0, 1, . . . , n, pak platí (1) Ij = par(n + 1)In-1-j pro j = -1, 0, 1, . . . , n. (2) Ij-1 = (j -xn+1) (j -x0) Ij pro vhodné j, xj < j < xj+1, 0 j < n 2 . (3) |Ij-1| > |Ij| pro 0 j < n 2 . (4) qj(x) := x xj n+1(u) du nemění znaménko na [xj, xn-j], pro n sudé, -1 j < n 2 . Je-li n liché, -1 j < n-1 2 , pak qj(x) nemění znaménko na [xj, xn-1-j]. Důkaz. Viz A. Ralston: Základy numerické matematiky, Academia Praha, 1973, str. 178, cv. 49. Věta 6.12 (NEWTON-COTESOVY KVADRATURNÍ FORMULE). Necht' x0, x1, . . . , xn jsou ekvidistantní uzly, xi = x0 + ih, h > 0 a a = x0 - h, b = xn + h, kde = 0 (v případě (n + 1)-bodové uzavřené formule) nebo = 1 (v případě (n + 1)-bodové otevřené formule). Potom existuje [a, b] takové, že platí b a f(x) dx = n i=0 Aif(xi) + hn+3 f(n+2) () (n + 2)! n+ - t n+1(t) dt Ef (6.4a) pro n sudé a f Cn+2 [a, b] a b a f(x) dx = n i=0 Aif(xi) + hn+2 f(n+1) () (n + 1)! n+ - n+1(t) dt Ef (6.4b) pro n liché a f Cn+1 [a, b], kde Ai = h (-1)n-i i!(n - i)! n+ - n+1(t) t - i dt (6.4c) jsou tzv. Cotesova čísla, pro něž platí An-i = Ai, i = 0, 1, . . . , n. Poznámka 6.13. ˇ n sudé (6.4a) = = n + 1 . . . stupeň přesnosti, Ef = h+2 enf(+1) () = O(h+2 ) = O(hn+3 ) pro h 0, en = 1 (n+2)! n+ - t2 (t - 1) . . . (t - n) dt. ˇ n liché (6.4b) = = n . . . stupeň přesnosti, Ef = h+2 enf(+1) () = O(h+2 ) = O(hn+2 ) pro h 0, en = 1 (n+1)! n+ - t(t - 1) . . . (t - n) dt. Vidíme, že pro lichý počet uzlů (n sudé) je řád konvergence i stupeň přesnosti Newton-Cotesových formulí o 1 vyšší, než pro sudý počet uzlů (n liché). Poznámka 6.14 (Speciální případy Newton-Cotesových formulí pro W(x) 1). Označíme-li Ai = hAai, kde ai = an-i, dostaneme následující tabulky popisující NUMERICKÉ METODY 49 a) Uzavřené ( = 0) Newton-Cotesovy formule pro 1 n 8 Název formule n A a0 a1 a2 a3 a4 en Lichoběžníkové pravidlo 1 1/2 1 1 1 -1/12 Simpsonovo pravidlo 2 1/3 1 4 1 3 -1/90 Simpsonovo 3/8 pravidlo 3 3/8 1 3 3 1 3 -3/80 Booleovo pravidlo 4 2/45 7 32 12 32 7 5 -8/945 5 5/288 19 75 50 50 75 5 -275/12096 6 1/140 41 216 27 272 27 7 -9/1400 7 7/17280 751 3577 1323 2989 2989 7 -8183/518400 8 4/14175 989 5888 -928 10946 -4540 9 -2368/467775 Pomocí obecných vzorců z věty 6.12, resp. z poznámky 6.13, odvodíme jako ukázku koefi- cienty pro prvé dvě uzavřené kvadraturní formule z předchozí tabulky. (1) Lichoběžníkové pravidlo (2-bodová uzavřená formule): = 0, n = 1, f C2 [a, b] Podle 6.13 je = n = 1 a tedy Ef = O(h3 ). Odtud a z (6.4c) dostáváme: A0 = h (-1)1-0 0! (1 - 0)! 1 0 (t - 1) dt = -h (t - 1)2 2 1 0 = -h 0 - 1 2 = h 1 2 . A1 = A1-1 = A0 = h 1 2 A = 1 2 , a0 = a1 = 1. e1 = 1 (1 + 1)! 1 0 t(t - 1) dt = 1 2 t3 3 - t2 2 1 0 = 1 2 1 3 - 1 2 = - 1 12 . (2) Simpsonovo pravidlo (3-bodová uzavřená formule): = 0, n = 2, f C4 [a, b] Podle 6.13 je = n + 1 = 3 a tedy Ef = O(h5 ). Odtud a z (6.4c) dostáváme: A0 = h (-1)2-0 0! (2 - 0)! 2 0 (t - 1)(t - 2) t3-3t+2 dt = h 1 2 t3 3 - 3t2 2 + 2t 2 0 = h 1 2 8 3 - 12 2 + 4 = h 1 3 . A1 = h (-1)2-1 1! (2 - 1)! 2 0 t(t - 2) dt = -h t3 3 - t2 2 0 = -h 8 3 - 4 = h 4 3 . A2 = A0 = h 1 3 , A1 = h 4 3 A = 1 3 , a0 = a2 = 1, a1 = 4. e2 = 1 (2 + 2)! 2 0 t2 (t - 1)(t - 2) dt = 1 24 t5 5 - 3t4 4 + 2t3 3 2 0 = 1 24 32 5 - 3 16 4 + 2 8 3 = = 4 15 - 1 2 + 2 9 = 24 - 45 + 20 90 = - 1 90 . Analogicky bychom postupovali pro n > 2. Zejména tedy dostáváme následující tvary nejčastěji používaných uzavřených Newton-Cotesových kvadraturních formulí. Lichoběžníkové pravidlo (n = 1): b a f(x) dx = h 2 (f0 + f1) - h3 12 f (), (6.5a) Simpsonovo pravidlo (n = 2): b a f(x) dx = h 3 (f0 + 4f1 + f2) - h5 90 f(4) (), (6.5b) 50 VÍTĚZSLAV VESELÝ Simpsonovo 3/8 pravidlo (n = 3): b a f(x) dx = 3h 8 (f0 + 3f1 + 3f2 + f3) - 3h5 80 f(4) (), (6.5c) Booleovo pravidlo (n = 4): b a f(x) dx = 2h 45 (7f0 + 32f1 + 12f2 + 32f3 + 7f4) - 8h7 945 f(6) () (6.5d) pro vhodné [a, b]. b) Otevřené ( = 1) Newton-Cotesovy formule pro 0 n 4 Název formule n A a0 a1 a2 en Obdélníkové pravidlo 0 2 1 1 1/3 1 3/2 1 1 1 3/4 2 4/3 2 -1 2 3 14/45 3 5/24 11 1 1 3 95/144 4 3/10 11 -14 26 5 41/140 Pomocí obecných vzorců z věty 6.12, resp. z poznámky 6.13, opět odvodíme jako ukázku koeficienty pro prvé dvě otevřené kvadraturní formule z předchozí tabulky. (1) Obdélníkové pravidlo (1-bodová otevřená formule): = 1, n = 0, f C2 [a, b] Podle 6.13 je = n + 1 = 1 a tedy Ef = O(h3 ). Odtud a z (6.4c) dostáváme: A0 = h (-1)0-0 0! (0 - 0)! 1 -1 1 dt = h(1 - (-1)) = h2 A = 2, a0 = 1. e0 = 1 (0 + 2)! 1 -1 t2 dt = 1 2 t3 3 1 -1 = 1 2 1 3 + 1 3 = 1 3 . (2) 2-bodová otevřená formule: = 1, n = 1, f C2 [a, b] Podle 6.13 je = n = 1 a tedy Ef = O(h3 ). Odtud a z (6.4c) dostáváme: A0 = h (-1)1-0 0! (1 - 0)! 2 -1 (t - 1) dt = -h (t - 1)2 2 2 -1 = -h 1 2 - 4 2 = h 3 2 . A1 = A1-1 = A0 = h 3 2 A = 3 2 , a0 = a1 = 1. e1 = 1 (1 + 1)! 2 -1 t(t - 1) dt = 1 2 t3 3 - t2 2 2 -1 = 1 2 8 3 - 4 2 + 1 3 + 1 2 ) = 1 2 3 - 3 2 = 3 4 . Analogicky bychom postupovali pro n > 1. Zejména tedy dostáváme následující tvary nejčastěji používaných otevřených Newton-Cotesových kvadraturních formulí. Obdélníkové pravidlo (n = 0): b a f(x) dx = 2hf0 + h3 3 f (), (6.6a) n = 1: b a f(x) dx = 3h 2 (f0 + f1) + 3h3 4 f (), (6.6b) NUMERICKÉ METODY 51 n = 2: b a f(x) dx = 4h 3 (2f0 - f1 + 2f2) + 14h5 45 f(4) (), (6.6c) n = 3: b a f(x) dx = 5h 24 (11f0 + f1 + f2 + 11f3) + 95h5 144 f(5) () (6.6d) pro vhodné [a, b]. Věta 6.15 (SLOŽENÉ NEWTON-COTESOVY KVADRATURNÍ FORMULE). Necht' x0, x1, . . . , xN , N = (n+2)M (M 1, n 0 celá) jsou ekvidistantní uzly, xj = x0+jh, h > 0 a a = x0-h, b = xn+h, kde = 0 (v případě uzavřené formule) nebo = 1 (v případě otevřené formule). Potom existuje [a, b] takové, že platí b a f(x) dx = N j=0 Ajf(xj) NC(f,h) + hn+2 b - a n + 2 enf(n+2) () ENC (f,h)=O(hn+2) (6.7a) pro n sudé, f Cn+2 [a, b] a en z poznámky 6.13 a b a f(x) dx = N j=0 Ajf(xj) NC(f,h) + hn+1 b - a n + 2 enf(n+1) () ENC (f,h)=O(hn+1) (6.7b) pro n liché, f Cn+1 [a, b] a en z poznámky 6.13, kde Aj = 2(1 - )A0 pro j = m(n + 2), 0 < m < Mcelé A(j mod(n+2))- pro ostatní j, 0 j N, (6.7c) přičemž A-1 = 0 a Ai jsou určeny vztahem (6.4c) pro i = 0, 1, . . . , n. Poznámka 6.16 (Speciální případy složených Newton-Cotesových formulí). ˇ Uzavřené složené formule ( = 0): (1) Lichoběžníkové složené pravidlo (NC= T . . . Trapezoidal, n = 1) 1 2 1 2 1 2 1 2 . . . 1 2 ... 1 2 1 2 Aj : 1 2 1 1 . . . 1 1 2 × h viz (6.5a) a (6.7c) b a f(x) dx = T(f, h) + ET (f, h), kde užitím (6.5a) a (6.7b) T(f, h) = h 2 (f0 + 2f1 + + 2fN-1 + fN ) = h f(a) + f(b) 2 + N-1 j=1 f(xj) , N = M, ET (f, h) = -h2 b - a 12 f () = O(h2 ), [a, b]. (6.8a) 52 VÍTĚZSLAV VESELÝ Algoritmus (MATLAB). jmathews('chap 7'): traprl.m (+ demo a7 1.m) (2) Simpsonovo složené pravidlo (NC= S, n = 2) 1 3 4 3 1 3 1 3 4 3 1 3 . . . 1 3 ... 1 3 4 3 1 3 Aj : 1 3 4 3 2 3 4 3 2 3 . . . 2 3 4 3 1 3 × h viz (6.5b) a (6.7c) b a f(x) dx = S(f, h) + ES(f, h), kde užitím (6.5b) a (6.7a) S(f, h) = h 3 (f0 + 4f1 + 2f2 + 4f3 + + 2fN-2 + 4fN-1 + fN ) = = h 3 f(a) + f(b) + 2 M-1 m=1 f(x2m) + 4 M m=1 f(x2m-1) , N = 2M, ES(f, h) = h4 b - a 2 - 1 90 e2 f(4) () = -h4 b - a 180 f(4) () = O(h4 ), [a, b]. (6.8b) Algoritmus (MATLAB). jmathews('chap 7'): simprl.m (+ demo a7 2.m) ˇ Otevřené složené formule ( = 1): Obdélníkové složené pravidlo (NC= R . . . Rectangular, n = 0) 0 2 0 0 2 0 . . . 0 ... 0 2 0 Aj : 0 2 0 2 0 . . . 0 2 0 × h viz (6.6a) a (6.7c) b a f(x) dx = R(f, h) + ER(f, h), kde užitím (6.6a) a (6.7a) R(f, h) = 2h(f1 + f3 + + fN-1) = 2h M m=1 f(x2m-1), N = 2M, ER(f, h) = h0+2 b - a 0 + 2 1 1 3 e0 f () = h2 b - a 6 f () = O(h2 ), [a, b]. (6.9a) NUMERICKÉ METODY 53 Katedra aplikované matematiky, Masarykova univerzita v Brně, Janáčkovo nám. 2a, 662 95 Brno E-mail address: vesely@@math.muni.cz