Numerická derivace a integrace
1 Úvod
Při řešení praktických příkladů často narazíme na problém. V mnoha případech známe funkční hodnoty nějaké funkce, ale neznáme její analytické zadání. Někdy se může stát, že funkce má tak složitý tvar, že nejsme schopni vyjádřit analyticky její derivaci nebo ji zintegrovat. V takových případech nám nezbude nic jiného, než použít numerické metody pro výpočet derivací neb integrálů. Jsou to tedy metody odhadu derivace nebo integrálu funkce. Pro řešení některých praktických problémů se nám nabízí i metody Monte Carlo. My si ale ukážeme další metody.
2 Numerická derivace
2.0.1   První derivace
f(x+h) f(x)
x x+h
Obrázek 1: Grafické znázornění derivace funkce f(x) v bodě x.
Vrátíme se na okamžik k definici derivace funkce jedné proměnné. Pomůže nám obrázek 1. Zde máme graf funkce fix) a vybereme si dva body x a x + h, pro které známe funkční hodnoty. Přímka, která prochází body fix) a f (x + h) svírá s osou x úhel, jehož tangens je definován jako ^x^]tI^ tedy ÍÍ2±ÍLl£M_ pf[ vhodně malém h přejde přímka do tečny ke
1
grafu v bodě x a my už víme, že směrnice tečny v bodě je derivací funkce v daném bodě. Připomeneme si nám již velmi dobře známou rovnici:
f (x) = hm---
Protože někdy nemáme k dispozici analytické vyjádření funkce, ale známe její funkční hodnoty v určených bodech, musíme použít numerické metody. V takovém případě proložíme známými body polynom, který následně derivujeme a vypočtená derivace polynomu v bodě odpovídá hodnotě derivace původní funkce v daném bodě. Stupeň polynomu musí být větší než stupeň počítané derivace. Cím vyšší je stupeň polynomu , tím přesnější je odhad funkce a tím je také přesnější výpočet derivace. Záleží na tvaru funkce, ale v některých případech dosahujeme uspokojivých výsledků i interpolací polynomem prvního stupně - tedy přímkou.
Pro odvození použijeme interpolační polynom v Lagrangeově tvaru 1 kdy budeme nahrazovat funkci f(x) v intervalu (xj; xl + h). Máme zadané dva body o souřadnicích [xf, f(xl)] a [xí + h; f(xl + h)], potom Lagrangeův polynom prvního stupně má tvar:
L1(x) = f(xl)X~^ + kKf(xl + h)- X~Xl
+ h) xl + h — xl
sečtením jmenovatelů získáme zjednodušený tvar:
Ll{x) = f(Xlf~iX] + h) + f(Xi + hf~ Xl
-h J v 1      ' h
Grafem získané funkce Li(x) je přímka. Derivací funkce Li(x) podle proměnné x do-
staneme:
L[(x) = ^.(f(xt + h)- fix,)) odtud plynou vzorce pro derivaci funkce fix) v bodech
Příklad 1 Vypočítejte hodnotu první derivace funkce f(x) = ^ v bodě x0 = 2 a pro h = 0, 2. Tuto jednoduchou funkci dokážeme snadno zderivovat a vypočítat přesnou hodnotu derivace v bodě x = 2. Derivace funkce f'(x) = —^2, tedy pro x = 2 je /'(2) = —0.25. Představme si ale, že funkci neumíme zderivovat nebo ji vůbec neznáme a známe jen její funkční body. Při použití hodnot x0 = 2 a pro x0 + h = 2,2 dostaneme pro derivaci hodnotu:
/'(„> = =-0,2273
lagrangeův polynom je jedním ze známějších a také snadných způsobů interpolace funkce zadané pouze v diskrétních bodech. Teorie k tomuto tématu přesahuje rámec naší výuky. Pro zájemce doporučuji vysokoškolské učebnice matematické analýzy.
2
nebo při použití hodnot xq = 2 a pro xq — h = 1, 8 dostaneme hodnotu:
/'(,,) = /P>-/(M)= _0| 277g
0 ^ 2
Vidíme, že výsledky se liší už na druhém desetinném místě. Je to ale důsledek použití toho nejjednoduššího vzorce, jaký jsme mohli použít.
Pokud budeme mít k dispozici tři body se souřadnicemi: [xl — h; fixl — h)], [xt; f(xt)] a [xí + h; f(xl + h)], můžeme sestrojit Lagrangeův polynom druhého stupně:
yX     XijyX     Xj ilj
L2{x) = f{x - h)
+f(x
+ f(Xi + h)
[X, —
(x — Xj + /i)(x — Xj — h) {x% - x% + -Xi-h)
(x — Xj + /l)(x — Xj)
(xj + /i — Xj + h) (xj + /i — Xj úpravou tohoto výrazu do mocnin (x — x i) dostaneme:
L2(x) = (Z~Ji)2(/(^ " >0 " 2/(xť) +       + h))
X"Xl(/(xl + /l)-/(xl-/l)) + /(xl)
2/i
Po derivaci funkce podle x dostaneme výraz:
Ľ2(x) = ^ [/(x, - h)(2(x - x,) - h) - 2/(xl).2(x - x,) + f(Xi - h)(2(x - x i) + h)]
Nyní do vyjádřené derivace L'2(x) dosadíme po řadě body xl5 (xj — h) a (xj + h) a dostaneme vztahy:
L'2(x, - h) = -^(-3/(zí - /i) + 4/(x,) - / (x, + /*)),
£'2(^) = ^(ffa + /i) - /(a;i - h)),
L'2(xí + h) = ^(f(xi -h)- Af(Xi) + 3/(x, + fc)) Po dosazení xo = Xj — h, x\ = Xj a x2 = Xj + h dostaneme následující vzorce:
f(xo) = ^(-3/(^o) + 4/(xx) - /(x2)), (2)
3
/'(*i) = ^(/(*2)-/(*<>)), (3)
nx2) = if(x0)-Af(x1) + 3f(x2)). (4)
Pro běžné výpočty je důležitý vzorec 3. Další dva vzorce 2 a 4 se používají tehdy, pokud nejsme schopni vypočítat funkční hodnoty nalevo od klíčového bodu, ve kterém nás derivace zajímá (uvedené vzorce pro derivaci v konkrétním bodě totiž užívají funkční hodnoty ve třech různých bodech).
Příklad 2   Nyní zkusíme výpočet jako v předchozím příkladu za použití aproximace Lagrangeovým polynomem druhého stupně. Hodnoty xq — 2 či h = 0,2. Po dosazení do vzorce 3 dostaneme:
/'(2) = äÍ^(/(2, 2) - /(l, 8)) = -0, 2525
Vidíme, že při stejných vstupních datech výsledek s přesností na dvě desetinná místa. Pokud bychom zmenšili krok h z hodnoty 0,2 na hodnotu 0,1 dostaneme pro polynom prvního stupně hodnotu:
f (2) = /(2' ^ ~ /(2) = -0, 2381 Pro polynom druhého stupně dostaneme:
/'(2) = 2^(/(2,1) " /(l, 9)) = -0, 2506
Jak jsme předpokládali zmenšením kroku pro funkční hodnoty obdržíme přesnější výsledek. To samé platí pro interpolaci polynomem vyššího řádu. Je ale nutné zopakovat, že pro interpolační polynom ž-tého stupně potřebujeme i + 1 bodů.
Příklad 3 Vypočtěte přibližnou hodnotu první derivace funkce fix) = ex(l — x) v bodě xq = 1 pro h = 0,1.
Nejprve použijeme vzorce la??.
x     el'Hl - 1,1) -e^l - 1) n
f(l) =-i--i-'- = -e1'1 w -3,0041
J ( ; 0,1
^     6 (1-1 -e ' (1-0,9) no
f(l) = —i-'--i-^ = -e0'9 w -2,4596
J y ' 0,1
4
Při použití vzorce 3 obdržíme:
/'(l) = -±-(^(1 - 1,1) - e°'9(l - 0, 9)) = el'1(-°'1} " e°'9(°'V>
0,2V    v       ' '        v       ' " 0,2
ei,i + eo
-2,7318
0,2 2
Protože tuto funkci umíme zderivovat, můžeme vypočítat přesnou hodnotu derivace funkce v bodě 1.
f{x) = ex{\ -x) + ex{-\) = -xex
/'(l) = -le1 = -e w -2,7182
Jak jsme předpokládali, nej lepšího výsledku dosáhneme při použití vzorce 3. První dva vzorce dosahují chyby v řádu velikosti h, tedy desetin. Poslední vzorec dosahuje chyby v řádu h2, tedy v našem případě v řádu setin.
2.1   Druhá derivace
Analogicky odvodíme druhou derivaci. Výraz L'2(xí) zderivujeme ještě jednou a obdržíme výraz:
Ľfai) = 7^(/(x* ~h)- 2f(xi) + f(xi + h))
Analogickým dosazením jako v případě první derivace xq — Xi — hj Xi — Xi ěl x2 — Xi -\- h dostaneme vzorec pro výpočet druhé derivace v bodě x\\
f"(x1) = ^(f(x0)-2f(x1) + f(x2)) (5)
Příklad 4   Vypočítejte druhou derivaci funkce f(x) = ^ pro x = 3 a pro h = 0, 2 a h = 0,1 podle vzorce 5. Pro h = 0, 2 dostáváme:
/"(3) =       (/(2, 8) - 2./(3, 0) + /(3, 2)) = 0, 0744
Pro h = 0,1 dostáváme:
/"(3) =       (/(2, 9) - 2./(3, 0) + /(3,1)) = 0, 0742
Protože opět danou funkci umíme zderivovat, můžeme vypočítat zadanou derivaci analyticky. Druhá derivace funkce f(x) = - je rovna f"(x) = \ a po dosazení hodnoty x = 3 dostáváme /"(3) = 0, 07407. Pro oba dva kroky jsme obdrželi shodu na tři desetinná místa.
5
3   Numerická integrace
Existuje několik důvodů, proč určitý integrál nepočítáme přesně, například
• integrál funkce f(x) neumíme spočítat analytickými metodami;
• analytický výpočet je příliš pracný;
• funkce fix) je dána jen tabulkou.
Z předchozích přednášek víme, že určitý integrál (Riemannův integrál) funkce fix) na intervalu (a; b) je obsah plochy ohraničené osou x, grafem funkce f(x) a kolmicemi na osu x v bodech a a b, viz Obrázek 2.
a b x
Obrázek 2: Grafické znázornění významu určitého integrálu funkce f(x) na intervalu (a; b)
Připomeňme si odvození Riemannova integrálu (Obrázek 3). Interval si rozdělíme na stejné úseky a v nich vytvoříme obdélníky s výškou nejnižší funkční hodnoty v daném intervalu a ve druhém případě nejvyšší funkční hodnotou. Součet obsahů obdélníků vytvořených s nej nižšími funkčními hodnotami se nazývá dolní součet, součet obsahů obdélníků s nejvyššími funkčními hodnotami se nazývá horní součet. V případě rozdělení integračního intervalu na limitně vysoký počet úseků, dojde ke sblížení velikosti dolních a horních součtů, ty jsou pak rovny hledanému určitému integrálu.
3.1    Obdélníková metoda
Tato metoda je nejpodobnější Riemannovu odvození. Nepoužíváme ale horní a dolní součty, místo nich vybereme funkční hodnotu uprostřed intervalu. Na obrázku 4 máme znázorněný celý interval (a; b). V praxi jej ale rozdělíme na stejné části. Čím menší úseky zvolíme, tím přesnější hodnotu získáme. Počet podintervalů na které rozdělujeme zadaný interval se nazývá dělení.
6
Obsah vyznačeného obdélníku vypočítáme S = (b — a)./(^-). V praxi interval {a; b) rozdělíme na m stejných úseků o šířce h. Funkční hodnotu v bodě a označíme f(x0), funkční hodnotu v bodě b označíme f(xm). Body dělení budou: xq = a, x±, x2, £m-i, xm = b. Pro výpočet integrálu potřebujeme znát funkční hodnoty středů jednotlivých podintervalů. Potom odhad určitého integrálu vypočítáme podle následujícího vztahu:
f(x)dx = h.(f(---) + /(---) + ... + /(---)) = /i- -2-^
Příklad 5 Za použití obdélníkové metody vypočítejte integrál J^_ľ exdx pro dělení m = 4 a m = 8. Výsledky porovnejte. Analyticky vypočítaná hodnota je 2,3504024.
Řešení:
Pro m = 4 bude šířka intervalu /i = (1 — (—1)) : 4 = 0,5. Interval (—1; 1) rozdělíme na body: xq = —1, x\ = —0,5, x2 = 0, x% = 0,5 a x4 = 1. Pro výpočet integrálu potřebujeme znát funkční hodnoty v bodech xt = —0, 75; —0, 25; 0, 25; 0, 75. Hledaný integrál tedy vypočítáme následovně:
I' exdx = 0, 5 * (e"0'75 + e"0'25 + e0'25 + e0'75) = = 0, 5 * (0,472367 + 0, 778801 + 1, 284025 + 2,117000) = 0, 5 * 4, 652193 = 2, 326097
Pro m = 8 bude šířka intervalu h = 0,25. Uzlové body intervalu budou: xq = —1, x\ = —0, 75, x2 = —0, 50, x3 = —0, 25, x4 = 0, x5 = 0, 25, x6 = 0, 50, x 7 = 0, 75 a x§ = 1. Pro výpočet integrálu budeme potřebovat funkční hodnoty: x, = -0, 875; -0, 625; -0, 375; -0,125; 0,125; 0, 375; 0, 625; 0, 875
7
[a+b)/2 b
Obrázek 4: Grafické znázornění obdélníkové metody řešení integrálu funkce fix) na intervalu (a;b)
Integrál bude mít tvar:
I' ďdx = 0, 25 * (e"0'875 + e"0'625 + e"0'375 + e"0'125 + e0'125 + e0'375 + e0'625 + e0'875) =
= 0, 25 * (0, 416862 + 0, 535261 + 0, 687289 + 0, 882497 + 1,133148 + 1,454991+
+1, 868246 + 2, 398875) = 0, 25 * 9, 377169 = 2, 344292
Při porovnání s analyticky vypočítanou hodnotou vidíme, že získané výsledky liší už na místě setin. Jak je vidět už z obrázku 4, tato metoda dosahuje dobré přesnosti až při velkém dělení. Proto si ukážeme další přesnější metodu.
3.2   Lichoběžníková metoda
U této metody budeme plochu pod křivkou nahrazovat stejně širokými lichoběžníky. Metoda je znázorněna pro jeden interval na obrázku 5. Obsah lichoběžníku vypočítáme podle nám dobře známého vzorce S = \iz\ + z2), kde v je výška lichoběžníku a z\ a z2 jsou jeho základny. V případě našeho integrálu odpovídají základny funkčním hodnotám v krajních bodech a výška lichoběžníka je šířkou intervalu, v případě dělení je tedy výška rovna hodnotě h.
V našem případě pro aproximaci jedním lichoběžníkem dostáváme výraz: Sj{x)dx = b-^(f(a)+f(b)) = ^(f(a) +/(&))
8
a b
Obrázek 5: Grafické znázornění lichoběžníkové metody řešení integrálu funkce f(x) na intervalu (a;b)
Budeme postupovat jako v předchozím případě. Funkční hodnotu v bodě a označíme f(xo), funkční hodnotu v bodě b označíme fixm). Body dělení budou: xq = a, x±, x2, xm-i, xm = b. Pro výpočet potřebujeme znát pouze šířku intervalů a funkční hodnoty v jejich krajních bodech. Potom pro aproximaci určitého integrálu dostáváme výraz:
^6/W^ = |[/(a;o)+/(a;1)]+|[/(^i)+/(^2)]+. - ■+^[/(^m-2)+/(a:m-i)]+|[/(a:m-i)+/(a:rra)]
Ve výrazu můžeme vytknout | a také vidíme, že každý bod uvnitř intervalu se ve výrazu vyskytuje dvakrát. Proto získáme zjednodušený tvar:
ľ f(x)dx = ^(f(xQ)) + 2 * [f(Xl) + f(x2) + ... + /(zm_i)] + f(xm))
Příklad 6 Nyní vyzkoušíme vypočítat určitý integrál /_1 exdx lichoběžníkovou metodou pro dělení m = 4 a m = 8. Výsledky porovnáme s předchozím příkladem.
Pro dělení m = 4 je šířka podintervalů h = 0,5. Uzlovými body jsou xq = —1, x\ = —0,5, x2 = 0, x3 = 0, 5 a x4 = 1. Pro výpočet integrálu budeme potřebovat funkční hodnoty v těchto uzlových bodech. Hledaná aproximace určitého integrálu má hodnotu:
9
/  e*dx = M(e-i + 2 * e"0-5 + 2 * e° + 2 * e°<5 + e1) = J-i 2
= —(O, 367879 + 2*0, 606531 + 2* 1 + 2 * 1, 648721 + 2, 718282) = = 0,25 * 9, 596665 = 2, 399166
Pro m = 8 bude šířka intervalu h = 0,25. Uzlové body intervalu budou: x0 = —1, x\ = —0, 75, x2 = —0, 50, x% = —0, 25, x4 = 0, x5 = 0, 25, xq = 0, 50, x 7 = 0, 75 a x§ = 1. Integrál má potom tvar:
ľ ďdx = ^ * (e-i + 2 * (e-°>75 + e"0-5 + e-°>25 + e° + e°'25 + e°'5 + e°>75) + e1) = J-i 2
= 0,125*(0, 367879+2*(0,472367+0, 606531+0, 778801+1+1, 284025+1, 648721+2,117000) +
+2, 718282) = 0,125 * 18, 901051 = 2, 362631
Vidíme, že i při použití této metody jsou rozdíly již na místě setin. Poslední metodou, kterou si představíme bude Simpsonova metoda.
3.3   Sipsonova metoda
V lichoběžníkové metodě jsme nahrazovali graf funkce f(x) přímkou - polynomem prvního stupně. Nyní ji nahradíme polynomem stupně druhého - kvadratickou funkcí. Pro tuto metodu je ale třeba mít interval rozdělený na sudý počet podintervalů. To je z důvodu, že budeme vytvářet dvojice těchto podintervalů. Pro každou sousední dvojici máme k dispozici 3 body a k nim příslušné 3 funkční hodnoty a těmi proložíme polynom druhého stupně - parabolu. Tento polynom zintegrujeme na obou intervalech součtem integrálů přes dvojice podintervalů získáme přibližnou hodnotu hledaného integrálu. Grafické znázornění této metody najdeme na obrázku 6. Zde tedy nebudeme řešit problém výpočtu obsahu rovinného obrazce jako v předchozích případech.
Interpolační polynom pro dané tři body xt — h, xt a xt + h má tvar:
L2(x) = {X [f(Xl -h)- 2f(Xi) + f(Xi + h)] +
+ X-^1[f^ + h)-f{xl-h)} + f{xl)
10
a (a+b)/2 b
Obrázek 6: Grafické znázornění Simpsonovy metody řešení integrálu funkce f{x) na intervalu (a;b)
Tuto funkci zintegrujeme na intervalu xl+i) a protože jsou podintervaly stejně
široké, můžeme při úpravě použít: h = — xl = xl — xl_i. Po dosazení mezí do zinte-grované funkce a dosazením dříve vyplývající podmínky 2h = xl+i —       obdržíme tvar:
/      f(x)dx = /      L2ix)dx = + 4/(xj) +
Sečtením výsledků přes všechny páry při dělení m dostáváme výraz:
řf(x)dx = ^(/(x0)+4/(x1) + 2/(x2) + 4/(x3) + ... + 2/(xm_2) + 4/(xm_1 + /(xm)) (6)
Pro snazší zapamatování výrazu (6) je tu následující pomůcka: poslední a první koeficienty jsou rovny 1, druhý a předposlední jsou rovny 4 a od třetího se střídají koeficienty 2 a 4.
Příklad 7 Náš již dobře známý příklad s integrálem funkce j\ exdx s dělením m = 4 a m = 8 za použití Simpsonovy metody.
Řešení: Pro dělení m = 4 je šířka podintervalů h = 0,5. Uzlovými body jsou x0 = —1, xi = —0, 5, x2 = 0, x3 = 0, 5 a x4 = 1. Pro výpočet integrálu budeme potřebovat funkční hodnoty v těchto uzlových bodech. Hledaná aproximace určitého integrálu má hodnotu:
ľ exdx = ^(e"1 + 4 * e"0'5 + 2 * e° + 4 * e0'5 + e1) = J-i 3
= —(0, 367879 + 4*0, 606531 + 2* 1 +4* 1, 648721 + 2, 718282) =
11
= 0,16667 * 14,107169 = 2, 351199
Pro m = 8 bude šířka intervalu h = 0,25. Uzlové body intervalu budou: xq = —1, x\ = —0, 75, x2 = —0, 50, x% = —0, 25, x4 = 0, x5 = 0, 25, xq = 0, 50, x 7 = 0, 75 a x§ = 1. Integrál má potom tvar:
ľ ďdx= 2l^*(e-l+4*e-0J5 + 2#e-0,5+4#e-0,25 + 2#e0+4#e0,25+2jKe0,5+4j|ee0,75+el) =
J-i 3
= 0, 0833333*(0, 367879+4*0, 472367+2*0, 606531+4*0, 778801+2*1+4*1, 284025+2*1, 648721+
+4 * 2,117000 + 2, 718282) = 0, 0833333 * 28, 205437 = 2, 350452
Vidíme, že jsme dosáhli mnohem lepších výsledků při stejném dělení při srovnání s obdélníkovou a lichoběžníkovou metodou. Při dělení m = 8 máme shodu s analytickým integrálem na čtvrtém desetinném místě.
Pro lepší procvičení Vám nabízím další příklady, které si můžete vypočítat sami.
1. Vypočtěte integrál Jq1 j^dx při dělení m = 10 za použití všech uvedených metod.
2. Vypočtěte integrál Jq ex2dx při dělení m = 10 za použití všech uvedených metod.
3. Za pomocí lichoběžníkové metody vypočtěte integrál Jq1 x2+\x+8dx pro dělení m = 4 a m = 8.
12