Fakulta elektrotechniky a informatiky Katedra aplikované matematiky. P. Beremlijski, M. Sadowská. Matematika s radostí

Vysoká ˇskola bánˇská – Technická univerzita Ostrava Fakulta elektrotechniky a informatiky Katedra aplikované matematiky

ˇn´ı u ´ lohy pro radost Aplikac P. Beremlijski, M. Sadowsk´ a

Matematika s radost´ı

Petr Beremlijski, Marie Sadowská Aplikaˇcn´ı u ´ lohy pro radost

c Petr Beremlijski, Marie Sadowská, 2015

´ Uvod Tato kniha je urˇcena pro kaˇzdého, koho zaj´ımaj´ı aplikace matematiky a kdo se rád matematikou bav´ı. Kniha ˇctenáˇre provede ˇradou matematických aplikac´ı. Jej´ım obsahem je také ˇreˇsen´ı tˇechto aplikac´ı. Text je urˇcen zejména pro studenty posledn´ıch roˇcn´ık˚ u stˇredn´ıch ˇskol, ale m˚ uˇze být zaj´ımavou inspirac´ı i pro jejich uˇcitele a dalˇs´ı zájemce o uˇziteˇcné a zaj´ımavé aplikace matematiky. Tento text vznikl v letech 2013–2015 a je u ´ zce spjatý s u ´ lohami, které ˇ ˇreˇs´ı studenti v pr˚ ubˇehu semináˇre Skola matematického modelován´ı. Tento ˇ semináˇr, pro který pouˇz´ıváme zkratku SKOMAM, organizuje Katedra aplikované matematiky Fakulty elektrotechniky a informatiky Vysoké ˇskoly bán ˇ ské – Technické univerzity Ostrava jednou roˇcnˇe jiˇz od roku 2005. Podrobnosti o semináˇri naleznete na webové adrese http://skomam.vsb.cz. Text byl napsán d´ıky finanˇcn´ı podpoˇre projektu Matematika s radost´ı. D´ıky realizaci tohoto projektu vznikla ˇrada interaktivn´ıch test˚ u, párovac´ıch her, her Neriskuj, AZ kv´ız˚ u a Poznej obrázek. Vˇsechny tyto materiály jsou na webových stránkách projektu http://msr.vsb.cz. Naˇse kniha je rozdˇelena do ˇsesti kapitol. V prvn´ı kapitole se vˇenujeme u ´ lohám z matematické analýzy. Na u ´ vod této kapitoly se seznám´ıme s nˇekterými základn´ımi pojmy a ty pak vyuˇzijeme napˇr´ıklad pro pˇribliˇzné urˇcen´ı ˇc´ısla π, výpoˇcet pˇribliˇzného obsahu zadané plochy nebo nalezen´ı pˇribliˇzného ˇreˇsen´ı zadané rovnice. Druhá kapitola je zamˇeˇrena na aplikace lineárn´ı algebry. Po u ´ vodu do této oblasti si ukáˇzeme, jak m˚ uˇzeme vyuˇz´ıt základn´ı algebraické struktury jako jsou vektory a matice pro vyhledáván´ı v databáz´ıch pomoc´ı kl´ıˇcových slov. Dalˇs´ı aplikac´ı, kterou se v této kapitole zabýváme, je realizace geometrických zobrazen´ı (stejnolehlost, rotace, posunut´ı) pomoc´ı nástroj˚ u lineárn´ı algebry. T´ımto zp˚ usobem m˚ uˇzeme geometrická zobrazen´ı efektivnˇe implementovat do poˇc´ıtaˇcové grafiky. V tˇret´ı kapitole se seznám´ıme se speciáln´ım typem rovnic, tzv. diferenciáln´ımi rovnicemi. A pomoc´ı nich si zavedeme tzv. poˇcáteˇcn´ı u ´ lohy. V u ´ vodu této kapitoly si nav´ıc ukáˇzeme, jak tyto u ´ lohy ˇreˇsit. Poˇcáteˇcn´ı u ´ lohy maj´ı ˇsiroké uplatnˇen´ı a popisuje se pomoc´ı nich ˇrada u ´ loh z oblasti biologie, fyziky iii

´ Uvod ˇci spoleˇcenských vˇed. My si ukáˇzeme, jak popsat napˇr´ıklad koncentraci léku v krvi pacienta v závislosti na ˇcase, jak modelovat pohyb kyvadla, ˇci jak odhadovat vývoj populace daného státu. Ve ˇctvrté kapitole se opˇet potkáme s diferenciáln´ımi rovnicemi, ale budeme zde ˇreˇsit tzv. okrajové u ´ lohy. Pomoc´ı takového typu u ´ loh si m˚ uˇzeme modelovat pr˚ uhyb struny pevnˇe upevnˇené na obou jejich konc´ıch, na kterou p˚ usob´ı zadaná s´ıla. V páté kapitole se pod´ıváme na optimalizaˇcn´ı u ´ lohy a ukáˇzeme si, jakými metodami lze tyto u ´ lohy ˇreˇsit. Pomoc´ı tohoto aparátu si pak pop´ıˇseme kuliˇcku ˇ sen´ım optimalizaˇcn´ı u dané hmotnosti zavˇeˇsenou na pruˇzinˇe. Reˇ ´ lohy z´ıskáme polohu, ve které se kuliˇcka ustál´ı. Nejrozsáhlejˇs´ı ˇsestá kapitola obsahuje podrobná a komentovaná ˇreˇsen´ı u ´ loh popsaných v prvn´ıch pˇeti kapitolách. Souˇcást´ı ˇreˇsen´ı jsou i kódy pro poˇc´ıtaˇcové ˇreˇsen´ı tˇechto u ´ loh. Pro poˇc´ıtaˇcové ˇreˇsen´ı jsme zvolili komerˇcn´ı systém Matlab, proto jsou vˇsechny kódy v této kapitole zapsány s pouˇzit´ım syntaxe tohoto systému. Matlab se skládá z nˇekolika souˇcást´ı. Jde zejména o programovac´ı jazyk zamˇeˇrený na numerické výpoˇcty a vývoj numerických algoritm˚ u. Matlab obsahuje také knihovnu funkc´ı pro ˇreˇsen´ı ˇrady numerických u ´ loh. A d˚ uleˇzitou sloˇzkou Matlabu je grafické prostˇred´ı pro interaktivn´ı zadáván´ı pˇr´ıkaz˚ u. Podrobný popis tohoto jazyka je k dispozici v [7]. Pokud ˇctenáˇr naˇsi knihy nemá Matlab k dispozici, m˚ uˇze pouˇz´ıt systém Octave, který se znaˇcnˇe podobá Matlabu a je volnˇe k dispozici na webové adrese ˇ sen´ı jsou nav´ıc doplnˇena tam, kde https://www.gnu.org/software/octave. Reˇ je to vhodné, o analýzu chyby aproximace ˇreˇsen´ı. Na závˇer je kniha doplnˇena nˇekolika pˇr´ılohami.

Ostrava, u ´ nor 2015

Petr Beremlijski Marie Sadowsk´ a

Text byl vytvoˇren v rámci projektu Matematika s radost´ı (http://msr.vsb.cz).

iv

Podˇ ekov´ an´ı Rádi bychom na tomto m´ıstˇe podˇekovali naˇs´ım koleg˚ um z Katedry aplikované matematiky Fakulty elektrotechniky a informatiky Vysoké ˇskoly bán ˇ ské – Technické univerzity Ostrava za pozorné pˇreˇcten´ı celého textu. Velký d´ık patˇr´ı zejména J. Bouchalovi a B. Krajcovi za mnohé cenné pˇripom´ınky ke ´ covi a A. Markopoukapitolám 1, 3 a 6. Dále dˇekujeme L. Malému, R. Cosi´ losovi za podnˇety ke kapitolám 2 a 4–6. M. Theuer nám pomohl vylepˇsit kapitolu o geometrických transformac´ıch a M. Pˇetroˇs z Fakultn´ı nemocnice Ostrava nám ochotnˇe poskytl rady ohlednˇe modelován´ı zmˇeny koncentrace lék˚ u v krvi pacient˚ u, moc jim za to dˇekujeme. P. B., M. S.

v

vi

Obsah ´ Uvod

iii

Podˇ ekov´ an´ı

v

Znaˇ cen´ı

ix

´ 1 Ulohy z matematick´ e anal´ yzy 1.1 Nˇekolik definic na u ´ vod: limita a derivace 1.2 Numerický výpoˇcet urˇcitého integrálu . . 1.3 Sˇc´ıtáme ˇrady - proˇc? . . . . . . . . . . . 1.4 Monte Carlo . . . . . . . . . . . . . . . . 1.5 Numerické ˇreˇsen´ı rovnic . . . . . . . . . 1.6 Jak pracuje kalkulaˇcka – k ˇcemu je dobrý

funkce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Taylor˚ uv polynom?

1 1 4 9 12 17 22

´ 2 Ulohy z line´ arn´ı algebry 2.1 Nˇekolik definic na u ´ vod: matice, souˇcin matic a inverze matice 2.2 Jednoduchý model vyhledáván´ı v databáz´ıch . . . . . . . . . . 2.3 Geometrické transformace . . . . . . . . . . . . . . . . . . . .

27 27 29 33

3 Poˇ c´ ateˇ cn´ı u ´ lohy 3.1 Struˇcný u ´ vod do poˇcáteˇcn´ıch u ´ loh a jak je ˇreˇsit 3.2 Modelován´ı zmˇeny koncentrace lék˚ u v krvi . . . 3.3 Tak je to padˇelek nebo to nen´ı padˇelek aneb jak nˇekterých Vermeerových“ obraz˚ u? . . . . . . . ” 3.4 Populaˇcn´ı modely . . . . . . . . . . . . . . . . . 3.5 Zmˇena teploty tˇelesa . . . . . . . . . . . . . . . 3.6 RC obvod . . . . . . . . . . . . . . . . . . . . . 3.7 Matematické kyvadlo . . . . . . . . . . . . . . .

37 . . . . . . . . 37 . . . . . . . . 39 poznat stáˇr´ı . . . . . . . . 42 . . . . . . . . 46 . . . . . . . . 47 . . . . . . . . 47 . . . . . . . . 49

4 Okrajov´ eu ´ lohy 53 4.1 Jednorozmˇerné (1D) okrajové u ´ lohy a 1D metoda s´ıt´ı . . . . . 53 vii

OBSAH 5 Optimalizaˇ cn´ı u ´ lohy 57 5.1 Nˇekolik definic na u ´ vod: gradient, Hessián . . . . . . . . . . . 57 5.2 Rovnováˇzná poloha tˇelesa . . . . . . . . . . . . . . . . . . . . 58 ˇ sen´ı u 6 Reˇ ´ loh ´ 6.1 Ulohy kapitoly ´ 6.2 Ulohy kapitoly ´ 6.3 Ulohy kapitoly ´ 6.4 Ulohy kapitoly ´ 6.5 Ulohy kapitoly

viii

1 2 3 4 5

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

65 65 85 92 112 115

Znaˇ cen´ı A\B N Z R R+ R∗ Rn Df Hf lim an lim f (x)

... ... ... ... ... ... ... ... ... ... ...

rozd´ıl mnoˇzin A a B mnoˇzina vˇsech pˇrirozených ˇc´ısel mnoˇzina vˇsech celých ˇc´ısel mnoˇzina vˇsech reálných ˇc´ısel mnoˇzina vˇsech kladných reálných ˇc´ısel rozˇs´ıˇrená ˇc´ıselná osa, tj. R∗ = R ∪ {+∞, −∞} mnoˇzina vˇsech uspoˇrádaných n-tic reálných ˇc´ısel definiˇcn´ı obor funkce f obor hodnot funkce f limita posloupnosti (an ) limita funkce f v bodˇe x0

f (x) f ′′ (x) f (n) (x) Rb f (x) dx a ∞ P an

... ... ... ...

(prvn´ı) derivace funkce f v bodˇe x druhá derivace funkce f v bodˇe x derivace n-tého ˇrádu funkce f v bodˇe x urˇcitý integrál funkce f od a do b

x→x0 ′

. . . ˇrada reálných ˇc´ısel

n=1

Uδ (x0 ) kxk

. . . okol´ı bodu x0 o polomˇeru δ . . . velikost vektoru x

ix

x

Kapitola 1 ´ Ulohy z matematick´ e anal´ yzy 1.1

Nˇ ekolik definic na u ´ vod: limita a derivace funkce

V této kapitole si velmi struˇcnˇe pˇripomeneme základn´ı pojmy matematické analýzy, které budeme v dalˇs´ıch kapitolách pouˇz´ıvat. N´ıˇze uvedené pojmy si lze mnohem podrobnˇeji prostudovat napˇr. v [1]. ˇ Nejprve si pˇripomeneme definici limity posloupnosti reálných ˇc´ısel. Rekneme, ˇze posloupnost (an ) m´ a limitu a ∈ R, a p´ıˇseme lim an = a nebo an → a, jestliˇze1 (∀ε ∈ R+ )(∃n0 ∈ N)(∀n ∈ N, n > n0 ) : |an − a| < ε. Má-li posloupnost (koneˇcnou) limitu, ˇr´ıkáme, ˇze je konvergentn´ ı. V opaˇcném ˇ pˇr´ıpadˇe se jedná o divergentn´ ı posloupnost. Rekneme, ˇze posloupnost (an ) m´ a limitu +∞, a p´ıˇseme lim an = +∞ nebo an → +∞, jestliˇze (∀k ∈ R)(∃n0 ∈ N)(∀n ∈ N, n > n0 ) : an > k.

ˇ Rekneme, ˇze posloupnost (an ) m´ a limitu −∞, a p´ıˇseme lim an = −∞ nebo an → −∞, 1

Pˇripomeˇ nme si, ˇze symbolem ∀ oznaˇcujeme vˇ seobecn´ y kvantifik´ ator a ˇcteme jej pro vˇsechna“. Symbol ∃ oznaˇcuje existenˇ cn´ ı kvantifik´ ator a ˇcteme jej existuje“. ” ”

1

´ Ulohy z matematické analýzy jestliˇze (∀ℓ ∈ R)(∃n0 ∈ N)(∀n ∈ N, n > n0 ) : an < ℓ. Pomoc´ı limity posloupnosti definujeme limitu reálné funkce jedné reálné promˇenné. Vztahu x0 6= xn → x0 rozum´ıme tak, ˇze xn → x0 a xn 6= x0 pro vˇsechna dost velká n ∈ N. Obdobnˇe budeme chápat i vztahy x0 < xn → x0 a ˇ x0 > xn → x0 . Rekneme, ˇze funkce f : R 7→ R m´ a v bodˇ e x0 ∈ R∗ limitu ∗ a ∈ R , a p´ıˇseme lim f (x) = a, x→x0

pokud x0 6= xn → x0

⇒

f (xn ) → a.

ˇ Rekneme, ˇze funkce f m´ a v bodˇ e x0 ∈ R limitu zprava a ∈ R∗ , a p´ıˇseme lim f (x) = a,

x→x0 +

pokud x0 < xn → x0

⇒

f (xn ) → a.

ˇ Rekneme, ˇze funkce f m´ a v bodˇ e x0 ∈ R limitu zleva a ∈ R∗ , a p´ıˇseme lim f (x) = a,

x→x0 −

pokud x0 > xn → x0

⇒

f (xn ) → a.

Plat´ı, ˇze funkce f má v bodˇe x0 ∈ R∗ nejvýˇse jednu limitu. ˇ Pojem limita funkce vyuˇzijeme dále pro definici spojitosti funkce. Rekneme, ˇze funkce f je spojit´ a v bodˇ e x0 ∈ R, jestliˇze lim f (x) = f (x0 ).

x→x0

ˇ Rekneme, ˇze funkce f je spojit´ a zprava v bodˇ e x0 ∈ R, jestliˇze lim f (x) = f (x0 ).

x→x0 +

ˇ Rekneme, ˇze funkce f je spojit´ a zleva v bodˇ e x0 ∈ R, jestliˇze lim f (x) = f (x0 ).

x→x0 −

ˇ Rekneme, ˇze funkce f je spojit´ a v intervalu I ⊂ R, plat´ı-li souˇcasnˇe: • f je spojitá v kaˇzdém vnitˇrn´ım bodˇe intervalu I, 2

• patˇr´ı-li poˇcáteˇcn´ı bod intervalu I do I, je v nˇem f spojitá zprava, • patˇr´ı-li koncový bod intervalu I do I, je v nˇem f spojitá zleva. Pomoc´ı limity funkce budeme definovat pojem derivace funkce. Bud’ f : R 7→ R a x ∈ R. Existuje-li f (x + h) − f (x) , h→0 h lim

nazveme ji derivac´ ı funkce f v bodˇ e x a znaˇc´ıme ji f ′ (x). Nen´ı-li ˇreˇceno jinak, budeme derivac´ı rozumˇet vlastn´ı (tj. koneˇcnou) derivaci. Lze ukázat, ˇze má-li funkce f v bodˇe x0 ∈ R (koneˇcnou) derivaci, je v bodˇe x0 spojitá. Derivac´ ı funkce f budeme dále rozumˇet funkci f ′ definovanou pˇredpisem f ′ (x) := f ′ (x). Zadefinujme si jeˇstˇe derivace vyˇsˇs´ıch ˇrád˚ u. Pˇredpokládejme nejprve, ˇze funkce f má derivaci v bodˇe x0 ∈ R. Jestliˇze funkce derivace f ′ má derivaci v bodˇe x0 , definujeme druhou derivaci funkce f v bodˇ e x0 jako f ′′ (x0 ) := (f ′)′ (x0 ). Podobnˇe postupujeme pˇri definic´ıch tˇ ret´ ı, ˇ ctvrt´ e, . . . derivace funkce f v bodˇ e x0 . Indukc´ı tedy definujeme pro n ∈ N derivaci n-t´ eho ˇ r´ adu funkce f v bodˇ e x0 jako f (n) (x0 ) := (f (n−1) )′ (x0 ), pˇriˇcemˇz f (0) := f . Na závˇer této kapitoly si jeˇstˇe uved’me nˇekterá základn´ı pravidla pro derivován´ı funkc´ı. Bud’ f, g : R 7→ R a x ∈ R. Pak plat´ı • (f ± g)′ (x) = f ′ (x) ± g ′ (x), má-li pravá strana rovnosti smysl, • (f g)′(x) = f ′ (x)g(x) + f (x)g ′(x), existuj´ı-li (vlastn´ı) derivace f ′ (x) a g ′ (x), ′ f f ′ (x)g(x) − f (x)g ′ (x) • (x) = , existuj´ı-li (vlastn´ı) derivace f ′ (x) g g 2 (x) a g ′ (x) a je-li g(x) 6= 0. 3

´ Ulohy z matematické analýzy Povˇezme si také, jak derivovat sloˇzenou funkci.2 Necht’ x ∈ R a necht’ existuj´ı koneˇcné derivace f ′ (x) a g ′ (f (x)). Potom je (g ◦ f )′ (x) = g ′ (f (x))f ′ (x). ˇ aˇri se nav´ıc mohou hodit pˇredpisy pro derivován´ı nˇekterých bˇeˇznˇe pouCten´ ˇz´ıvaných funkc´ı: • (c)′ = 0, c ∈ R (konst.), x ∈ R, • (xr )′ = rxr−1 , r ∈ R, x ∈ R+ ,3 • (sin x)′ = cos x, x ∈ R, • (cos x)′ = − sin x, x ∈ R, nπ o 1 , x∈R\ + kπ : k ∈ Z , • (tg x) = cos2 x 2 ′

• (cotg x)′ = −

1 , x ∈ R \ {kπ : k ∈ Z}, sin2 x

• (ex )′ = ex , x ∈ R, • (ln x)′ =

1.2

1 , x ∈ R+ . x

Numerick´ y v´ ypoˇ cet urˇ cit´ eho integr´ alu

Nˇ ekolik pozn´ amek k urˇ cit´ emu integr´ alu Funkci F : R 7→ R nazveme primitivn´ ı funkc´ ı k funkci f : R 7→ R na otevˇreném intervalu I ⊂ R, pokud pro kaˇzdé x ∈ I plat´ı F ′ (x) = f (x). Lze ukázat, ˇze pokud je f spojitá na otevˇreném intervalu I, pak existuje primitivn´ı funkce F k f na tomto intervalu. Pomoc´ı primitivn´ı funkce zadefinujme urˇcitý integrál funkce f od a do b, kde a, b ∈ R a a < b. Jsou-li funkce f a F spojité na uzavˇreném intervalu 2

Pˇripomeˇ nme, ˇze pro funkce f, g : R 7→ R takové, ˇze Hg ⊂ Df , definujeme sloˇ zenou funkci f ◦ g pˇredpisem (f ◦ g)(x) := f (g(x)). 3 Pˇripomeˇ nme si, ˇze x0 := 1.

4

ha, bi a je-li F primitivn´ı k f na (a, b), definujeme urˇ cit´ y integr´ al funkce f od a do b jako Z b

a

f (x) dx := F (b) − F (a).

(1.1)

Vztah (1.1) se ˇcasto nazývá Newton˚ uv-Leibniz˚ uv vzorec.4,5 Pro zjednoduˇsen´ı zápisu se pro rozd´ıl hodnot F (b) a F (a) zavád´ı toto znaˇcen´ı: ozn.

[F (x)]ba = F (b) − F (a). Popiˇsme si nyn´ı geometrický význam urˇcitého integrálu. Nejprve pˇredpokládejme, ˇze je funkce f spojitá a nezáporná na intervalu ha, bi. Pak integrál Z b f (x) dx a

je roven obsahu rovinného obrazce ohraniˇceného osou x, grafem funkce f a pˇr´ımkami x = a a x = b, viz obr. 1.1. Nyn´ı rozˇsiˇrme naˇse u ´ vahy na spojité y

f (x) Rb a

O

a

f (x) dx

b

x

Obrázek 1.1: Význam urˇcitého integrálu funkce, které mohou být v intervalu ha, bi záporné. Necht’ je napˇr. funkce f záporná na intervalu (c, d) ⊂ ha, bi, viz obr. 1.2. Potom také Z d f (x) dx < 0. c

Chceme-li pak pomoc´ı tohoto integrálu vypoˇc´ıst obsah plochy ohraniˇcené osou x, grafem funkce f a pˇr´ımkami x = c a x = d, mus´ıme na této ˇcásti vz´ıt integrál s opaˇcným znaménkem. Bude-li nás dále zaj´ımat obsah plochy, kterou ohraniˇcuje osa x, graf funkce f a pˇr´ımky x = a a x = b, a budeli funkce f prot´ınat v intervalu ha, bi osu x, mus´ıme tyto pr˚ useˇc´ıky naj´ıt 4 5

Isaac Newton (1643–1727) – v´ yznamn´ y anglick´ y fyzik, matematik a astronom Gottfried Wilhelm Leibniz (1646–1716) – v´ yznamn´ y nˇemeck´ y matematik

5

´ Ulohy z matematické analýzy a rozdˇelit interval ha, bi na intervaly, na nichˇz má f totéˇz znaménko. Na tˇechto intervalech spoˇcteme integrály funkce f . Výsledný obsah plochy pak dostaneme jako souˇcet vˇsech vypoˇctených integrál˚ u, pˇriˇcemˇz integrály na tˇech intervalech, kde je f nekladná, mus´ıme uvaˇzovat se záporným znaménkem (viz obr. 1.2). y f (x)

Rc a

O

a

f (x) dx

Rb d

c

Rd − c f (x) dx

d

f (x) dx b

x

Obrázek 1.2: Obsah rovinného obrazce a urˇcitý integrál Pˇ ribliˇ zn´ y v´ ypoˇ cet urˇ cit´ eho integr´ alu V pˇr´ıpadˇe, ˇze primitivn´ı funkci nelze vyjádˇrit elementárn´ımi nebo tabelovanými funkcemi, Newton˚ uv-Leibniz˚ uv vzorec (1.1) nem˚ uˇzeme pouˇz´ıt. Nˇekdy zase m˚ uˇze být hledán´ı primitivn´ı funkce pˇr´ıliˇs sloˇzité ˇci ˇcasovˇe nároˇcné. V tˇechto pˇr´ıpadech ˇcasto pˇristupujeme k tzv. numerickému výpoˇctu daného urˇcitého integrálu, který nám dá pˇribliˇznou hodnotu integrálu s danou pˇresnost´ı. V následuj´ıc´ım textu si pˇribl´ıˇz´ıme dvˇe základn´ı metody pˇribliˇzného výpoˇctu urˇcitého integrálu, a to obdéln´ıkové a lichobˇeˇzn´ıkové pravidlo. Obd´ eln´ ıkov´ e pravidlo Rozdˇelme nejprve interval ha, bi na n stejných d´ılk˚ u o délce

b−a . n Krajn´ı body i-tého d´ılku postupnˇe oznaˇcme xi−1 a xi ; plat´ı tedy h=

x0 = a < x1 < · · · < xn−1 < xn = b. Vypoˇcteme si dále stˇredy jednotlivých d´ılk˚ u: ci = 6

xi−1 + xi , 2

i = 1, 2, . . . , n.

(1.2)

(1.3)

Na i-tém d´ılku pak funkci f nahrad´ıme konstantn´ı funkc´ı o hodnotˇe f (ci ) a hledaný integrál budeme aproximovat takto: Z b n X f (x) dx ≈ hf (c1 ) + hf (c2 ) + · · · + hf (cn ) = h f (ci ). a

i=1

Napˇr´ıklad integrál funkce f od a do b z obr. 1.3 tak nahrazujeme souˇctem obsah˚ u pˇr´ısluˇsných obdéln´ık˚ u. y f (x)

O

a = x0

x1

x2

x3

x4 = b

x

Obrázek 1.3: Aproximace obdéln´ıky (n = 4) Lze ukázat, ˇze pokud existuje spojitá f ′′ na ha, bi, potom pro chybu aproximace plat´ı Z n b (b − a)3 X f (x) dx − h f (ci ) ≤ max |f ′′ (x)|. (1.4) 2 a 24n x∈ha,bi i=1

Lichobˇ eˇ zn´ ıkov´ e pravidlo

Interval ha, bi rozdˇel´ıme stejnˇe jako u obdéln´ıkového pravidla na d´ılky shodné délky, viz (1.2) a (1.3). Hledaný integrál budeme aproximovat takto: Z b f (x1 ) + f (x2 ) f (xn−1 ) + f (xn ) f (x0 ) + f (x1 ) f (x) dx ≈ h +h +· · ·+h = 2 2 2 a h h = (f (x0 )+2f (x1 )+· · ·+2f (xn−1 )+f (xn )) = 2 2

f (x0 ) + f (xn ) + 2

n−1 X

!

f (xi ) .

i=1

Napˇr´ıklad integrál funkce f od a do b z obr. 1.4 tak nahrazujeme souˇctem obsah˚ u pˇr´ısluˇsných lichobˇeˇzn´ık˚ u. Je moˇzné ukázat, ˇze pokud existuje spojitá f ′′ na ha, bi, potom pro chybu aproximace plat´ı Z ! n−1 b (b − a)3 X h f (x) dx − f (x0 ) + f (xn ) + 2 f (xi ) ≤ max |f ′′ (x)|. 2 x∈ha,bi a 2 12n i=1

7

´ Ulohy z matematické analýzy y f (x)

O

a = x0

x1

x2

x3

x4 = b

x

Obrázek 1.4: Aproximace lichobˇeˇzn´ıky (n = 4) Vˇsimnˇeme si, ˇze u lichobˇeˇzn´ıkového pravidla máme odhad chyby aproximace horˇs´ı neˇz u obdéln´ıkového pravidla, pˇrestoˇze u obdéln´ıkového pravidla uˇz´ıváme nahrazen´ı konstantn´ımi funkcemi a u lichobˇeˇzn´ıkového pravidla pouˇz´ıváme nahrazen´ı lineárn´ımi funkcemi.6 Z uvedených vzorc˚ u pro odhad chyby aproximace lze urˇcit, jaký poˇcet d´ılk˚ u zaruˇc´ı poˇzadovanou pˇresnost aproximace. Jelikoˇz ovˇsem odhady obsahuj´ı druhé derivace, jejichˇz hodnoty nen´ı vˇzdy lehké odhadnout, m˚ uˇze být výpoˇcet poˇctu potˇrebných d´ılk˚ u nároˇcný a výsledek pesimistický. Pˇ r´ıklad 1 Pomoc´ı obdéln´ıkového a lichobˇeˇzn´ıkového pravidla spoˇctˇete pˇribliˇznˇe obsah jednotkového kruhu a výsledky porovnejte se skuteˇcným obsahem daného kruhu. y 1

O

√

1 − x2

1

x

Obrázek 1.5: Jednotkový kruh

6

Aproximaˇcn´ı vlastnosti obdéln´ıkového a lichobˇeˇzn´ıkového pravidla si lze detailnˇeji prostudovat napˇr. v [10].

8

2

p

1

1 − (x − 1)2

y

0 −1 −2

arcsin(x − 1) −

−3 −4

−3

−2

−1

0

1

2

π 2

3

x

Obrázek 1.6: Matematické srdce Pˇ r´ıklad 2 Pomoc´ı obdéln´ıkového a lichobˇeˇzn´ıkového pravidla spoˇctˇete pˇribliˇznˇe obsah srdce na obr. 1.6 a výsledky porovnejte se skuteˇcným obsahem srdce. Pˇ r´ıklad 3 Aproximujte hodnotu integrálu Z 1 x2 dx 0

pomoc´ı obdéln´ıkového pravidla tak, aby chyba aproximace byla nejvýˇse 10−4 . Urˇcete poˇcet d´ılk˚ u dˇelen´ı intervalu h0, 1i, který zaruˇc´ı dosaˇzen´ı poˇzadované pˇresnosti.

1.3

Sˇ c´ıt´ ame ˇ rady - proˇ c?

V této kapitole si nejprve pov´ıme, co se skrývá pod pojmy reálná ˇc´ıselná ˇrada a jej´ı konvergence. ˇ Radou re´ aln´ ych ˇ c´ ısel rozum´ıme výraz a1 + a2 + · · · + an + · · · = kde an ∈ R pro kaˇzdé n ∈ N. Napˇr´ıklad: i)

∞ X n=1

7

∞ X

an ,

(1.5)

n=1

(−1)n = −1 + 1 + (−1) + 1 + · · · (pˇr´ıklad tzv. alternuj´ ıc´ ı ˇ rady7 ),

Alternuj´ıc´ı ˇrady jsou takové ˇrady, jejichˇz ˇcleny pravidelnˇe stˇr´ıdaj´ı znaménka.

9

´ Ulohy z matematické analýzy

ii)

∞ X

n = 1 + 2 + 3 + 4 + · · · (pˇr´ıklad tzv. aritmetick´ e ˇ rady),

∞ X

q n−1 = 1 + q + q 2 + q 3 + · · ·, q ∈ R (pˇr´ıklad tzv. geometrick´ e ˇ rady),

n=1

iii)

n=1

∞ X 1 1 1 1 iv) = 1 + + + + · · · (tzv. harmonick´ a ˇ rada). n 2 3 4 n=1

ˇ ıslo an nazýváme n-t´ C´ ym ˇ clenem ˇ rady (1.5), posloupnost (sn ) definovanou pˇredpisem n X sn := a1 + a2 + · · · + an = ak k=1

nazýváme posloupnost´ ı ˇ c´ asteˇ cn´ ych souˇ ct˚ u ˇ rady (1.5). Existuje-li limita lim sn = s ∈ R∗ , nazýváme ji souˇ ctem ˇ rady (1.5) a p´ıˇseme ∞ X

an = s.

n=1

Napˇr´ıklad: i) Souˇcet

∞ X

n

(−1) neexistuje, jelikoˇz sn =

n=1

ii) Souˇcet

∞ X

n = +∞, jelikoˇz sn =

(

0 pro n sudé, 1 pro n liché.

n(n+1) . 2

n=1

iii) Souˇcet

∞ X

q n−1

n=1

jelikoˇz sn =

(

 je-li q ≥ 1,   +∞, 1 , je-li |q| < 1, = 1−q   neexistuje, je-li q ≤ −1,

n

pro q = 1,

1−q n 1−q

pro q 6= 1.

iv) Lze ukázat, ˇze souˇcet 10

∞ X 1 = +∞. n n=1

ˇ ıkáme dále, ˇze ˇrada konverguje, jestliˇze je jej´ı souˇcet roven nˇejakému R´ reálnému ˇc´ıslu. V opaˇcném pˇr´ıpadˇe, tj. pokud má ˇrada souˇcet roven +∞ nebo −∞ nebo pokud ˇrada souˇcet nemá, ˇr´ıkáme, ˇze ˇrada diverguje. Napˇr´ıklad: i)

∞ X

(−1)n diverguje.

n=1

ii)

∞ X

n diverguje.

n=1

iii)

∞ X n=1

q n−1 konverguje, pokud |q| < 1, a diverguje, pokud |q| ≥ 1.

∞ X 1 iv) diverguje. n n=1

Výˇse uvedené pojmy m˚ uˇzeme nyn´ı vyuˇz´ıt k ˇreˇsen´ı následuj´ıc´ıch u ´ loh. Pˇ r´ıklad 1 Prozrad’me, ˇze ˇrada ∞ X n=1

1 n(n + 1)

konverguje a jej´ı posloupnost ˇcásteˇcných souˇct˚ u konverguje rychle“ k souˇctu ” této ˇrady. Zkuste odhadnout souˇcet této ˇrady. Pˇ r´ıklad 2 Mezi dvˇema mˇesty A a B se nacház´ı pˇr´ımá ˇzelezniˇcn´ı trat’, jej´ıˇz délka je 90 km. Z mˇesta A do mˇesta B vyjede vlak rychlost´ı 10 km/h. V tu samou chv´ıli vyjede z mˇesta B vlak do mˇesta A po té samé koleji stejnou rychlost´ı. Ve chv´ıli, kdy se vlaky rozjedou vstˇr´ıc jisté zkáze, z pˇredn´ıho okna lokomotivy vlaku jedouc´ıho z A do B se odraz´ı moucha rychlost´ı 100 km/h a let´ı vstˇr´ıc druhému vlaku. Ve chv´ıli, kdy k nˇemu dolet´ı, dotkne se noˇzkou jeho pˇredn´ıho skla a let´ı zpátky. Takto moucha l´ıtá mezi vlaky, neˇz ji rozmáˇcknou na placku pˇri jejich sráˇzce. Kolik kilometr˚ u moucha celkem nalétala? Jakou vzdálenost moucha urazila mezi 15. a 16. odraˇzen´ım se od oken lokomotiv (jinými slovy: jak dlouhý byl jej´ı 15. let mezi vlaky)? Pˇ r´ıklad 3 Kolik ˇclen˚ u harmonické ˇrady mus´ıte nejménˇe seˇc´ıst, aby tento ˇcásteˇcný souˇcet ˇrady mˇel hodnotu alespoˇ n 10 (15, 20)?

11

´ Ulohy z matematické analýzy Pˇ r´ıklad 4 Prozrad’me, ˇze ˇrada ∞ X (−1)n+1 4 2n − 1 n=1

konverguje a jej´ı posloupnost ˇcásteˇcných souˇct˚ u konverguje rychle“ k souˇctu ” této ˇrady. Zkuste odhadnout souˇcet této ˇrady.

1.4

Monte Carlo

V této kapitole se budeme vˇenovat odhad˚ um Ludolfova8 ˇc´ısla π. Nejprve hledejme obsah kruhu o polomˇeru r (pˇredstavme si na chv´ıli, ˇze neznáme pˇr´ısluˇsný vzorec a nic nev´ıme o ˇc´ıslu π). Prvn´ım nápadem by mohlo být zhotoven´ı válcových nádob s r˚ uznými polomˇery podstav (napˇr. s polomˇery o délce 1 a 2 jednotek) a jednotkovou výˇskou, viz obr. 1.7. Objem vody, který

v =1 r=1

v =1 r=2

Obrázek 1.7: Válcové nádoby se do takových nádob vejde, je roven obsahu podstavy válce, tj. obsahu kruhu s polomˇerem r. Rychle si vˇsimneme, ˇze pokud zvˇetˇs´ıme polomˇer podstavy dvakrát, zvˇetˇs´ı se objem ˇctyˇrikrát, a následnˇe pozorujeme, ˇze obsah kruhu je pˇr´ımo u ´ mˇerný druhé mocninˇe polomˇeru. Také zjist´ıme, ˇze druhou mocninu polomˇeru kruhu mus´ıme vynásobit vhodnou konstantou, abychom dostali správnou hodnotu obsahu daného kruhu. Tuto konstantu oznaˇc´ıme π. Existuje mnoho moˇznost´ı, jak tuto konstantu odhadnout. Snad nejjednoduˇsˇs´ım zp˚ usobem jak stanovit meze pro π, je vepsat do kruhu o polomˇeru r ˇctverec a stejnému kruhu opsat jiný ˇctverec, viz obr. 1.8. Délka strany vˇetˇs´ıho ˇctverce je menˇs´ıho ˇctverce je √ 2r a z Pythagorovy vˇety dále zjist´ıme, ˇze délka strany 2 2r, tud´ıˇz obsah vˇetˇs´ıho resp. menˇs´ıho ˇctverce je 4r resp. 2r 2. Jelikoˇz jsme si uˇz odvodili“, ˇze obsah kruhu o polomˇeru r je dán vzorcem πr 2 , pouhým ” porovnán´ım obsah˚ u ˇctverc˚ u a kruhu zjist´ıme, ˇze π ∈ (2, 4). Tuto metodu 8

12

Ludolph van Ceulen (1540–1610) – nˇemeck´ y matematik

r r

Obrázek 1.8: Exhaustn´ı metoda pro stanoven´ı mez´ı pro π nazýváme vyˇcerpávac´ı (exhaustn´ı) a pravdˇepodobnˇe prvn´ı ji pouˇzil Eudoxos.9 Neˇz se budeme vˇenovat dalˇs´ım odhad˚ um ˇc´ısla π, pod´ıvejme se krátce na výpoˇcet obvodu kruhu. Jistˇe v´ıme, ˇze obvod kruhu je pˇr´ımo u ´ mˇerný dvojnásobku jeho polomˇeru, ale abychom dostali správnou hodnotu, je nutno 2r vynásobit vhodnou konstantou k. Na obr. 1.9 provedeme pˇreuspoˇrádán´ı kruhu na u ´ tvar, který se pro zjemˇ nuj´ıc´ı se dˇelen´ı kruhu bl´ıˇz´ı obdéln´ıku. Porovnán´ım obsahu kruhu (πr 2 ) a obsahu vzniklého obdéln´ıku (kr 2 ) je názornˇe vidˇet, ˇze konstanta k je opˇet rovna π.

o = 2kr r

r kr

zjemnˇen´ı dˇelen´ı

r kr k=π

Obrázek 1.9: Stanoven´ı vzorce pro obvod kruhu Nyn´ı si ukáˇzeme nˇekolik zp˚ usob˚ u, jak nalézt pˇribliˇznou hodnotu ˇc´ısla π. 9

Eudoxos (410 nebo 408 pˇr. n. l. – 355 nebo 347 pˇr. n. l.) – ˇreck´ y astronom, matematik a fyzik, student Platóna

13

´ Ulohy z matematické analýzy Buffonova metoda ˇ sen´ım tzv. Buffonova10 problému s jehlou je aproximace ˇc´ısla π. Uloha ´ Reˇ spoˇc´ıvá v opakovaném házen´ı jehly o délce ℓ na rovinu (pap´ır), na které máme vyznaˇcenu s´ıt’ rovnobˇeˇzek se vzdálenost´ı 2ℓ. Jestliˇze jehlu hod´ıme n-krát a s-krát nám bˇehem tˇechto pokus˚ u po dopadu zkˇr´ıˇz´ı nˇekterou z rovnobˇeˇzek, pak ˇc´ıslo n s aproximuje ˇc´ıslo π. Nyn´ı si uved’me odhad týkaj´ıc´ı se pˇresnosti Buffonovy metody. Lze ukázat, ˇze s pravdˇ √ epodobnost´ı 95 procent nemá chyba aproximace hodnotu vˇetˇs´ı neˇz 9,011/ n. Tzn. napˇr´ıklad pro 104 pokus˚ u nám s pravdˇepodobnost´ı 95 procent chyba nepˇrekroˇc´ı hodnotu 0,09011. Pˇ r´ıklad 1 Implementujte Buffonovu metodu a pouˇzijte ji k aproximaci ˇc´ısla π. Porovnejte vaˇsi aproximaci se skuteˇcnou“ hodnotou ˇc´ısla π a urˇcete chybu ” aproximace. Metoda Monte Carlo Monte Carlo je tˇr´ıda výpoˇcetn´ıch algoritm˚ u zaloˇzená na provádˇen´ı náhodných experiment˚ u. Této metody se ˇcasto pouˇz´ıvá pro simulaci fyzikáln´ıch a matematických systém˚ u. Výsledkem proveden´ı velkého mnoˇzstv´ı experiment˚ u je obvykle pravdˇepodobnost urˇcitého jevu. Na základˇe z´ıskané pravdˇepodobnosti a známých vztah˚ u pak spoˇc´ıtáme potˇrebné výsledky. Protoˇze metoda vyˇzaduje generován´ı velkého souboru náhodných dat, je vhodné pro jej´ı implementaci pouˇzit´ı poˇc´ıtaˇce. Metod Monte Carlo se pouˇz´ıvá v pˇr´ıpadˇe, kdy je pˇr´ıliˇs pracné nebo nemoˇzné nalézt pˇresný výsledek jiným zp˚ usobem. Jej´ı výhodou je jednoduchá implementace, nevýhodou relativnˇe malá pˇresnost. Metoda byla vytvoˇrena skupinou fyzik˚ u pracuj´ıc´ıch na projektu jaderné pumy v Los Alamos, jméno metody bylo navrˇzeno v roce 1940 von Neumannem.11 V matematice se Monte Carlo pouˇz´ıvá zejména pro výpoˇcet urˇcitých integrál˚ u (zejména v´ıcenásobných urˇcitých integrál˚ u), které je obt´ıˇzné ˇci nemoˇzné vyˇc´ıslit analyticky nebo jinou numerickou metodou. Napˇr. obsah omezené plochy ohraniˇcené grafem funkce y = x2 , osou x a pˇr´ımkou x = 1 R1 E (tj. S(E) = 0 x2 dx, viz obr. 1.10) je moˇzné metodou Monte Carlo vypoˇc´ıst následuj´ıc´ım zp˚ usobem. Necht’ náˇs program generuje náhodnˇe dvojice ˇc´ısel 10 11

14

Georges Louis Leclerc de Buffon (1707–1788) – francouzsk´ y pˇr´ırodovˇedec John von Neumann (1903–1957) – v´ yznamn´ y mad’arsk´ y matematik

[x, y], pˇriˇcemˇz kaˇzdé z ˇc´ısel x a y je vybráno nezávisle z intervalu h0, 1i. Tuto dvojici budeme chápat jako souˇradnice bodu, který je náhodnˇe zvolen ve ˇctverci h0, 1i × h0, 1i. Pravdˇepodobnost toho, ˇze bod leˇz´ı uvnitˇr zadaného ˇctverce, je 1. Pravdˇepodobnost toho, ˇze bod leˇz´ı uvnitˇr podmnoˇziny E zadaného ˇctverce, je rovna obsahu plochy E, tj. S(E). Tedy obsah plochy, která je podmnoˇzinou zvoleného ˇctverce, m˚ uˇzeme odhadnout jako pravdˇepodobnost, ˇze náhodnˇe zvolený bod z daného ˇctverce leˇz´ı v této podmnoˇzinˇe.

1

1

S(E) =

y

y

x2 Z

1

x2 dx

0

0

0

1

x

0

0

1

x

Obrázek 1.10: Plocha E a Monte Carlo Pˇ 2 Implementujte metodu Monte Carlo pro pˇribliˇzný výpoˇcet R 1r´ıklad 2 x dx. Z´ıskané výsledky m˚ uˇzete porovnat s analytickým vyˇc´ıslen´ım in0 tegrálu. Pro odhad pˇresnosti metody Monte Carlo lze odvodit napˇr´ıklad následuj´ıc´ı tvrzen´ı. S pravdˇ nemá chyba aproximace hodnotu √ epodobnost´ı 75 procent 4 vˇetˇs´ı neˇz 1/ n. Tzn. napˇr´ıklad pro 10 pokus˚ u nám s pravdˇepodobnost´ı 75 procent chyba nepˇrekroˇc´ı hodnotu 0,01. Chceme-li vyuˇz´ıt metodu Monte Carlo k aproximaci ˇc´ısla π, vypoˇcteme pˇribliˇznˇe touto metodou (napˇr´ıklad) integrál Z 1√ 1 − x2 dx. 0

Snadno si uvˇedom´ıme (viz obr. 1.11), ˇze t´ımto zp˚ usobem z´ıskáme aproximaci hodnoty π/4. Pˇ r´ıkladR 3√Implementujte metodu Monte Carlo pro pˇribliˇzný výpoˇcet in1 tegrálu 0 1 − x2 dx a pouˇzijte ji k aproximaci ˇc´ısla π. Porovnejte vaˇsi

15

´ Ulohy z matematické analýzy aproximaci se skuteˇcnou“ hodnotou ˇc´ısla π a urˇcete chybu aproximace. ” 1

1

y

1 − x2

y

√

Z

0

0

1

√

1 − x2 dx

0

1

0

0

1

x

x

Obrázek 1.11: π a Monte Carlo Pˇripomeˇ nme jeˇstˇe, ˇze výˇse uvedené pouˇzit´ı metody Monte Carlo je pouze ilustrativn´ı; jednoduchý integrál z pˇr´ıkladu 3 lze velmi efektivnˇe spoˇc´ıst jinými numerickými integraˇcn´ımi metodami. Aproximace π pomoc´ı ˇ c´ıseln´ ych ˇ rad Posledn´ı metodou nalezen´ı aproximace ˇc´ısla π, kterou si v tomto pˇrehledu ukáˇzeme, je vyuˇzit´ı ˇc´ıselných ˇrad. K této aproximaci pouˇzijeme tzv. Leibnizovu ˇradu, která je rozvojem funkce arkustangens: ∞

X x3 x5 x7 x2n−1 arctg x = x − + − +··· = (−1)n+1 . 3 5 7 2n − 1 n=1 Tato ˇrada má koneˇcný souˇcet pro x ∈ h−1, 1i (ˇr´ıkáme, ˇze ˇrada konverguje), nav´ıc plat´ı, ˇze ˇc´ım v´ıce je |x| bliˇzˇs´ı 0, t´ım ménˇe ˇclen˚ u ˇrady potˇrebujeme pouˇz´ıt k nahrazen´ı arctg x s uspokojivou“ pˇresnost´ı. Pro x = 1 dostáváme ” z Leibnizovy ˇrady tzv. Gregoryho12 ˇradu ∞

X (−1)n+1 1 1 1 1− + − +··· = . 3 5 7 2n − 1 n=1 Jelikoˇz arctg 1 = 12

16

π , 4

James Gregory (1638–1675) – skotsk´ y matematik a astronom

m˚ uˇzeme Gregoryho ˇradu vyuˇz´ıt k aproximaci ˇc´ısla π. Aproximaci s rychlejˇs´ı konvergenc´ı pak z´ıskáme, pokud pouˇzijeme rovnost arctg 1 = arctg

1 1 + arctg , 2 3

(1.6)

viz obr. 1.12.

arctg 1

arctg

1 2

arctg

1 3

Obrázek 1.12: Ilustrace k odvozen´ı vztahu (1.6)

Pˇ r´ıklad 4 Implementujte metodu, která vyuˇzije Gregoryho, popˇr. Leibnizovu ˇradu k nalezen´ı aproximace ˇc´ısla π.

1.5

Numerick´ eˇ reˇsen´ı rovnic

Velmi ˇcasto se (nejen) matematik setkává s u ´ kolem, kdy je tˇreba vyˇreˇsit nˇejakou nelineárn´ı rovnici. Tento u ´ kol bývá vˇetˇsinou sloˇzitý a ˇreˇsen´ı nelineárn´ı rovnice ˇcasto nelze naj´ıt analyticky. Proto je uˇziteˇcné vˇedˇet, jak lze rovnice ˇreˇsit alespoˇ n pˇribliˇznˇe pomoc´ı numerických metod. V následuj´ıc´ım textu si ukáˇzeme dvˇe jednoduché iteraˇcn´ı metody, které umoˇzn ˇ uj´ı nalézt alespoˇ n pˇribliˇzné ˇreˇsen´ı nelineárn´ı rovnice. Metoda prost´ e iterace Bod x∗ ∈ R nazýváme pevn´ ym bodem funkce f : R 7→ R, pokud f (x∗ ) = x∗. 17

´ Ulohy z matematické analýzy Vˇ eta Bud’ ∅ = 6 I ⊂ R uzavˇrený omezený interval a f : R 7→ R. Necht’ f je spojitá v I a necht’ f má spojitou derivaci v I. At’ nav´ıc existuje λ ∈ R takové, ˇze13 |f ′(x)| ≤ λ < 1 pro vˇsechna x ∈ I. Zvolme bod x0 ∈ I. Jestliˇze kaˇzdý bod posloupnosti (xn ) definované pˇredpisem xn+1 := f (xn ),

n = 0, 1, 2, 3, . . . ,

leˇz´ı v I, pak existuje x∗ ∈ I takové, ˇze f (x∗ ) = x∗ = lim xn , tj. (xn ) konverguje k pevnému bodu f . Lze ukázat, ˇze ˇc´ım v´ıce je hodnota nejmenˇs´ıho horn´ıho odhadu mnoˇziny {|f ′ (x)| : x ∈ I} bliˇzˇs´ı 0, t´ım rychleji posloupnost (xn ) konverguje k pevnému bodu f . Pro dané x0 uvaˇzujme posloupnost (xn ) danou pˇredpisem 1 xn+1 := xn − (x2n − 2), n = 0, 1, 2, 3, . . . . (1.7) 4 √ Vˇsimnˇeme si nejprve, ˇze ± 2 jsou pevné body funkce f (x) := x − 41 (x2 − 2). √ Lze si rozmyslet (viz také obr. 1.13), ˇze posloupnost (1.7) konverguje k 2 pro libovolné x0 ∈ (0, 4). 2

5

f (x)

1

4

−1

2

|f ′ (x)|

y

3

y

0

−2

1

−3

0

−4

−1

−5

−2

0

2

4

6

−2

−2

x

0

2

4

6

x

Obrázek 1.13: Grafy funkc´ı f (x) = x − 41 (x2 − 2) a |f ′ (x)| = 1 − 12 x

13

V pravém, resp. levém krajn´ım bodˇe intervalu I nepoˇzadujeme existenci derivace, ale pouze existenci derivace zleva, resp. zprava. Derivaci zleva, resp. zprava funkce f (x) (x) v bodˇ e x ∈ R definujeme jako lim f (x+h)−f , resp. lim f (x+h)−f . h h h→0−

18

h→0+

Algoritmus (Metoda prost´ e iterace) 1. ε > 0 (ukonˇcuj´ıc´ı podm´ınka) x0 (poˇcáteˇcn´ı aproximace) n=0 x1 = f (x0 ) 2. while |xn+1 − xn | ≥ ε n=n+1 xn+1 = f (xn ) end 3. xn+1 aproximuje ˇreˇsen´ı rovnice f (x) = x Na√obr. 1.14 je ilustrována konvergence posloupnosti (1.7) k pevnému bodu 2 funkce f pro x0 = 1.

2

f (x)

y

1.5

1

0.5

x 0

0

0.5

x0 x1x21.5

2.5

3.5

x

Obr´ prosté iterace k pevnému bodu √ azek 1.14: Ilustrace konvergence metody 1 2 2= ˙ 1,41421356 funkce f (x) = x − 4 (x − 2) Naleznˇeme nyn´ı pomoc´ı metody prosté iterace pˇribliˇzné ˇreˇsen´ı rovnice x

10 e− 10 − x = 0.

(1.8)

ˇ sen´ı rovnice (1.8) pˇrevedeme na ˇreˇsen´ı rovnice Reˇ x

10 e− 10 = x, x

tj. budeme hledat pevný bod fukce f (x) := 10 e− 10 . Lze si rozmyslet (viz také obr. 1.15), ˇze posloupnost (xn ) daná pˇredpisem xn

xn+1 := 10 e− 10 ,

n = 0, 1, 2, 3, . . . , 19

´ Ulohy z matematické analýzy konverguje k pevnému bodu funkce f pro libovolné x0 > 0. Na obr. 1.16 je ilustrována konvergence této posloupnosti k pevnému bodu funkce f pro x0 = 1.

1.2 12 1

10 8

y

y

|f ′ (x)|

0.8

f (x) 6

0.6 4 0.4 2 0

5

10

15

0.2

x

0

5

10

15

x

x x Obrázek 1.15: Grafy funkc´ı f (x) = 10 e− 10 a |f ′ (x)| = −e− 10 10

8

x

f (x) y

6

4

2

0

0

x0

2

4 x2

x4 x

6

x3

8

x1

10

Obrázek 1.16: Ilustrace konvergence metody prosté iterace k pevnému bodu x x∗ = ˙ 5,67143290 funkce f (x) = 10 e− 10 Newtonova metoda Zabývejme se nyn´ı ˇreˇsen´ım rovnice g(x) = 0.

(1.9)

Pˇredpokládejme, ˇze g je spojitá a má spojitou a nenulovou derivaci v nˇejakém okol´ı14 U(x∗ ), kde x∗ ˇreˇs´ı rovnici (1.9). Budeme vycházet z metody prosté 14

20

Okol´ım U(x∗ ) budeme rozumˇet nˇejak´ y otevˇren´ y omezen´ y interval se stˇredem v x∗ .

iterace. Protoˇze x∗ ˇreˇs´ı rovnici (1.9), tak x∗ ˇreˇs´ı také rovnici x = f (x) := x −

g(x) , h(x)

(1.10)

kde h je libovolná nenulová a diferencovatelná funkce na okol´ı U(x∗ ). Pak dále plat´ı, ˇze f ′ (x) = 1 −

g ′ (x)h(x) − g(x)h′ (x) h2 (x)

pro vˇsechna x ∈ U(x∗ ).

Pro dobrou“ konvergenci budeme poˇzadovat f ′ (x∗ ) = 0, tj. ” g ′(x∗ ) 1− = 0. h(x∗ ) Odtud h(x∗ ) = g ′(x∗ ). Vol´ıme proto h(x) = g ′ (x). Dostáváme tak, ˇze posloupnost (xn ) definovaná pˇredpisem xn+1 := xn −

g(xn ) , g ′ (xn )

n = 0, 1, 2, 3, . . . ,

(1.11)

konverguje k ˇreˇsen´ı x∗ rovnice (1.9), pokud x0 je dostateˇcnˇe bl´ızko x∗ , funkce f je daná (1.10) a jsou splnˇeny pˇredpoklady vˇety ze str. 18. Algoritmus (Newtonova metoda) 1. ε > 0 (ukonˇcuj´ıc´ı podm´ınka) x0 (poˇcáteˇcn´ı aproximace) n=0 x1 = x0 − g(x0 )/g ′(x0 ) 2. while |xn+1 − xn | ≥ ε n=n+1 xn+1 = xn − g(xn )/g ′ (xn ) end 3. xn+1 aproximuje ˇreˇsen´ı rovnice g(x) = 0 Pˇripomeˇ nme si jeˇstˇe, ˇze rovnice teˇcny t grafu funkce g v dotykovém bodˇe (xn , g(xn )) má tvar t:

y = g ′ (xn )(x − xn ) + g(xn ).

Hledáme-li pr˚ useˇc´ık pˇr´ımky t a osy x, ˇreˇs´ıme rovnici g ′ (xn )(x − xn ) + g(xn ) = 0.

(1.12) 21

´ Ulohy z matematické analýzy Je-li g ′(xn ) 6= 0, pak lze ˇreˇsen´ı rovnice (1.12) psát ve tvaru x = xn −

g(xn ) . g ′(xn )

Porovnáme-li pak posledn´ı rovnost s iteraˇcn´ı formul´ı (1.11), zjist´ıme, ˇze kaˇzdé xn+1 vytvoˇr´ıme jako pr˚ useˇc´ık teˇcny sestrojené ke grafu funkce g v dotykovém bodˇe (xn , g(xn )) a osy x. Proto se Newtonova metoda ˇcasto nazývá metodou teˇ cen. Na obr. 1.17 ilustrujeme pˇribliˇzné ˇreˇsen´ı rovnice g(x) := e−x − 1/7 = 0 pomoc´ı Newtonovy metody pro x0 = 1.

0.8 0.6

y

g(x) 0.4 0.2 0

−0.2

0

0.5

1x0

x22

1.5 x1

2.5

3

3.5

x

Obrázek 1.17: Ilustrace konvergence Newtonovy x∗ = ˙ 1,94591015 rovnice g(x) := e−x − 1/7 = 0

metody

k ˇreˇsen´ı

Pˇ r´ıklad 1 Pomoc´ı Newtonovy metody aproximujte ˇreˇsen´ı rovnic 1) ln x + (x − 1)3 = 0, 2

2) x + e−x = 1. Pˇ r´ıklad 2 Pomoc´ı Newtonovy metody najdˇete pˇribliˇzné hodnoty ˇc´ısel a π.

1.6

√

2

Jak pracuje kalkulaˇ cka – k ˇ cemu je dobr´ y Taylor˚ uv polynom?

Pˇri ˇreˇsen´ı r˚ uzných matematických u ´ loh se setkáváme s potˇrebou vyˇc´ıslit hodnotu nˇejaké funkce v daném bodˇe. V nˇekterých pˇr´ıpadech to ale m˚ uˇze být 22

√ obt´ıˇzné. Uvaˇzujeme-li napˇ √r´ıklad funkci x, pak urˇcit hodnotu této funkce v bodˇe 1 je jednoduché: 1 = 1. Jak n √ ale pomoc´ı desetinného ˇc´ısla (alespoˇ pˇribliˇznˇe) zapsat hodnotu funkce x napˇr´ıklad v bodˇe 2? Nab´ız´ı se ˇreˇsen´ı: pouˇzijeme kalkulaˇcku. Jak ovˇsem ˇc´ıslo 1,414213 . . . spoˇcetla kalkulaˇcka? V této kapitole si ukáˇzeme postup, jak lze pˇribliˇznˇe vypoˇc´ıst hodnoty obt´ıˇznˇe ” vyˇc´ıslitelných“ funkc´ı v ruce a simulovat t´ım práci jednoduché kalkulaˇcky. Nejprve si zavedeme nˇekolik pojm˚ u, s nimiˇz budeme dále pracovat. Má-li funkce f v bodˇe x0 derivace aˇz do ˇrádu n ∈ N vˇcetnˇe, definujeme Taylor˚ uv15 polynom ˇ e x0 vztahem r´ adu n funkce f v bodˇ Tn (x) := f (x0 ) + f ′ (x0 )(x − x0 ) + =

n X f (k) (x0 ) k=0

k!

f ′′ (x0 ) f (n) (x0 ) (x − x0 )2 + · · · + (x − x0 )n 2! n!

(x − x0 )k .

Okol´ım bodu x0 ∈ R o polomˇeru δ > 0 budeme rozumˇet otevˇrený interval (x0 − δ, x0 + δ) a budeme jej znaˇcit Uδ (x0 ). Nyn´ı budeme cht´ıt aproximovat funkci f na okol´ı Uδ (x0 ) Taylorovým polynomem,16 ˇc´ımˇz se dopust´ıme urˇcité chyby. Výhodou lokáln´ı aproximace funkce pomoc´ı polynomu je napˇr´ıklad to, ˇze funkˇcn´ı hodnotu polynomu lze spoˇc´ıst pouze pomoc´ı operac´ı sˇc´ıtán´ı a násoben´ı. Obvykle ˇc´ım vyˇsˇs´ı stupeˇ n Taylorova polynomu pak budeme pˇri aproximaci uvaˇzovat, t´ım menˇs´ı se dopust´ıme chyby (tj. t´ım pˇresnˇejˇs´ı nahrazen´ı z´ıskáme). Taylorova vˇ eta Pˇredpokl´ adejme, ˇze funkce f m´ a v kaˇzdém bodˇe okol´ı Uδ (x0 ) vˇsechny derivace aˇz do ˇr´ adu n + 1 vˇcetnˇe, a uvaˇzujme x ∈ Uδ (x0 ). Pak existuje ˇc´ıslo ξ leˇz´ıc´ı mezi x0 a x takové, ˇze plat´ı f (x) = Tn (x) + Rn+1 (x), kde Rn+1 – tzv. zbytek po n-t´ em ˇ clenu – je daný vztahem Rn+1 (x) =

f (n+1) (ξ) (x − x0 )n+1 . (n + 1)!

Uvedená podoba zbytku se nazývá Lagrange˚ uv17 tvar zbytku. Existuje v´ıce (nˇekdy uˇziteˇcnˇejˇs´ıch) tvar˚ u zbytku. Pˇresnost aproximace hodnoty funkce 15

Brook Taylor (1685–1731) – anglick´ y matematik Toto nahrazen´ı budeme zapisovat takto: f (x) ≈ Tn (x). 17 Joseph Louis Lagrange (1736–1813) – v´ yznamn´ y italsko-francouzsk´ y matematik a mechanik 16

23

´ Ulohy z matematické analýzy f v x ∈ Uδ (x0 ) hodnotou Taylorova polynomu Tn v x lze zjistit odhadem velikosti zbytku Rn+1 (x). Vrat’me √ se na chv´ıli zpˇet k naˇsemu problému jak (pˇribliˇznˇe) vyjádˇr√it hodnotu ˇc´ısla 2. Pouˇzijeme výˇse uvedených u ´ vah a funkci f (x) = x nahrad´ıme na okol´ı bodu x0 = 1 Taylorovým polynomem napˇr´ıklad stupnˇe 4. Dostaneme tak (viz také obr. 1.18) √

x ≈ T4 (x) = f (1)+f ′(1)(x−1)+

f ′′′ (1) f (4) (1) f ′′ (1) (x−1)2 + (x−1)3 + (x−1)4 , 2 6 24

a proto pro x = 2 máme

2

f (x) 1.5

T4 (x)

y

1 0.5 0 −0.5

0

1

2 x

3

4

Obrázek 1.18: Nahrazen´ı funkce f Taylorovým polynomem T4 v bodˇe 1 √

2 ≈ T4 (2) = f (1) + f ′ (1) +

f ′′ (1) f ′′′ (1) f (4) (1) + + . 2 6 24

Poˇc´ıtejme: 1 1 • f ′ (x) = √ , f ′ (1) = , 2 2 x 1 1 • f ′′ (x) = − √ , f ′′ (1) = − , 4 4 x3 3 3 • f ′′′ (x) = √ , f ′′′ (1) = , 5 8 8 x • f (4) (x) = − 24

15 15 √ , f (4) (1) = − . 16 16 x7

Tedy √

2≈1+

3 15 1 1 − + − = 1,3984375. 2 8 48 384

Nen´ı tˇeˇzké v tomto pˇr´ıpadˇe odhadnout velikost zbytku R5 (2). Protoˇze f (5) (x) =

105 √ , 32 x9

plat´ı R5 (2) ≤

f (5) (1) 7 (2 − 1)5 = = 0,02734375. 5! 256

√ Lepˇs´ı pˇresnosti vyˇc´ıslen´ı hodnoty 2 bychom dosáhli, pokud bychom zvolili √ vyˇsˇs´ı ˇrád Taylorova polynomu funkce x v bodˇe 1.18,19 Pˇ r´ıklad 1 Pomoc´ı Taylorova polynomu urˇcete pˇribliˇznˇe hodnotu ˇc´ısla e.20 Pˇ r´ıklad 2 Pomoc´ı Taylorova polynomu urˇcete pˇribliˇznˇe hodnotu ˇc´ısla π.21 Zkusme se zamyslet, zda by neˇslo naˇsi jednoduchou kalkulaˇcku nˇejak zrychlit, tedy zda m˚ uˇzeme dosáhnout stejného pˇribl´ıˇzen´ı k poˇzadované hodnotˇe s vyuˇzit´ım menˇs´ıho poˇctu aritmetických operac´ı. Odpovˇed’ zn´ı ano – m˚ uˇzeme pouˇz´ıt metodu prosté iterace, kterou jsme si jiˇz popsali v kapitole 1.5. Pˇripomeˇ nme si v tuto chv´ıli alespoˇ n daný iteraˇcn´ı algoritmus pro nalezen´ı ˇreˇsen´ı rovnice f (x) = x, tj. pro nalezen´ı pevného bodu funkce f , za pˇredpokladu, ˇze f splˇ nuje pˇredpoklady vˇety ze str. 18.

1 1 1 1 Zkuste si rozmyslet, ˇze Tn (2) = 1 + 1! · 12 + 2! · − 41 + 3! · 38 + 4! · − 15 16 + · · · n P (−1)k −1·1·3·5· ... ·(2k−3) 1 1 1 3 5 2n−3 + n! =1+ · . 2 · −2 · −2 · −2 · . . . · − 2 k! 2k 18

19

k=1

Leckdy je nam´ısto zvyˇsován´ı ˇrádu n v´ yhodnˇejˇs´ı lépe zvolit x0 , napˇr. v tomto pˇr´ıpadˇe je rozumnou volbou x0 = 49/25. 20 Nápovˇeda: Pouˇzijte Taylor˚ uv polynom stupnˇe n funkce f (x) = ex ve vhodném bodˇe x0 . 21 Nápovˇeda: Pouˇzijte Taylor˚ uv polynom stupnˇe n funkce f (x) = arctg x ve vhodném bodˇe x0 .

25

´ Ulohy z matematické analýzy Algoritmus (Metoda prost´ e iterace) 1. ε > 0 (ukonˇcuj´ıc´ı podm´ınka) x0 (poˇcáteˇcn´ı aproximace) n=0 x1 = f (x0 ) 2. while |xn+1 − xn | ≥ ε n=n+1 xn+1 = f (xn ) end 3. xn+1 aproximuje ˇreˇsen´ı rovnice f (x) = x Pˇripomeˇ nme si dále (viz str. 18), ˇze posloupnost (xn ) daná pˇredpisem xn+1 := xn − k(x2n − 2),

n = 0, 1, 2, 3, . . . ,

1 k := , 4

(1.13)

√ √ konverguje k 2 pro libovolné x0 ∈ (0, 4). ( 2 je pevným bodem funkce f (x) := x − k(x2 − 2).) Pˇ r´ıklad 3 Zmˇ ˇ te koeficient k v (1.13) tak, aby daná posloupnost kon√ en vergovala k 2 rychleji neˇz stávaj´ıc´ı posloupnost (1.13), pokud poˇcáteˇcn´ı aproximace je x√0 = 1. Zmˇen ˇ te koeficient k v (1.13) tak, aby posloupnost konvergovala k 2 pro x0 = 5. Pˇ r´ıklad 4 Pomoc´ı metody prostých iterac´ı urˇcete pˇribliˇznˇe hodnotu ˇc´ısla π.

26

Kapitola 2 ´ Ulohy z line´ arn´ı algebry 2.1

Nˇ ekolik definic na u ´ vod: matice, souˇ cin matic a inverze matice

Uvaˇzujme prvky x1 , x2 , . . . , xn z dané mnoˇziny F ⊂ R. Aritmetick´ y vektor (n-rozmˇerný) je uspoˇrádaná n-tice ˇc´ısel, jej´ıˇz prvky se nazývaj´ı sloˇzky. Tyto uspoˇrádané n-tice budeme zapisovat do hranatých závorek do ˇrádk˚ u nebo sloupc˚ u:

x = [x1 , x2 , . . . , xn ]



 x1   nebo x =  ...  . xn

Struˇcnˇe p´ıˇseme téˇz x = [xi ]. Nˇekdy je uˇziteˇcné zapisovat i-tou sloˇzku vektoru x jako [x]i , tedy xi = [x]i . Pro libovolný skalár α a aritmetický vektor x dále definujeme jejich souˇcin jako αx := [αxi ] . Sˇc´ıtán´ı a odeˇc´ıtán´ı dvou aritmetických vektor˚ u x a y o n sloˇzkách je dáno pˇredpisem [x ± y]i := [x]i ± [y]i , i = 1, 2, . . . , n a jejich skal´ arn´ ı souˇ cin definujme vztahem xy := x1 y1 + x2 y2 + · · · + xn yn . Uvaˇzujme prvky a11 , a12 , . . . , amn z dané mnoˇziny F ⊂ R. Matice typu 27

´ Ulohy z line´ arn´ı algebry (m, n) je obdéln´ıková tabulka 

 a11 . . . a1n  ..  , .. A =  ... . .  am1 . . . amn

která má m · n prvk˚ u aij uspoˇrádaných do m ˇrádk˚ u riA a n sloupc˚ u sA ze j , takˇ 

 r1A   A A =  ...  = sA , . . . , s 1 n , A rm riA = [ai1 , . . . , ain ] ,



 a1j  ..  sA . . j =  amj

Struˇcnˇe p´ıˇseme téˇz A = [aij ] a ij-tou sloˇzku matice A m˚ uˇzeme alternativnˇe zapsat jako [A]ij , tedy aij = [A]ij . Opˇet pro libovolný skalár α a matici A definujeme jejich souˇcin jako αA := [αaij ]. Sˇc´ıtán´ı a odeˇc´ıtán´ı dvou matic téhoˇz typu je dáno pˇredpisem [A ± B]ij := [A]ij ± [B]ij . Nyn´ı si zavedeme operaci násoben´ı matic. Nejdˇr´ıve ovˇsem zadefinujme souˇ cin matice a vektoru. Souˇcinem matice A = [aij ] typu (m, n) a sloupcového vektoru x = [xi ] o n sloˇzkách nazýváme sloupcový vektor y o m sloˇzkách definovaný pˇredpisem A y = Ax := x1 sA 1 + · · · + xn sn .

Rozepsán´ım této definice po sloˇzkách dostaneme yi = [Ax]i = ai1 x1 + · · · + ain xn = riA x. Nyn´ı pˇrejdˇeme k násoben´ı matic. Jestliˇze A je matice typu (m, ℓ) a B je matice typu (ℓ, n), pak souˇ cin matic A a B je matice AB typu (m, n) definovaná pˇredpisem B AB := AsB 1 , . . . , Asn . 28

Rozep´ıˇseme-li si definici násoben´ı matic po sloˇzkách, dostaneme [AB]ij = ai1 b1j + · · · + aiℓ bℓj = riA sB j . Nakonec si zavedeme inverzn´ı matici. Necht’ A je ˇctvercová matice (tzn. jde o matici, která má stejný poˇcet ˇrádk˚ u a sloupc˚ u). Jestliˇze existuje matice B tak, ˇze AB = BA = I, 1 ˇ pak se matice B nazývá inverzn´ ı matic´ ı k matici A a znaˇc´ı se A−1 . Ctvercová matice, ke které existuje inverzn´ı matice, se nazývá regul´ arn´ ı. V opaˇcném pˇr´ıpadˇe takovou matici nazýváme singul´ arn´ ı. Dá se ukázat, ˇze ke kaˇzdé regulárn´ı matici existuje právˇe jedna inverzn´ı matice. Mnohem v´ıce detail˚ u o vektorech, matic´ıch a operac´ıch s nimi lze nalézt napˇr´ıklad v [5].

2.2

Jednoduch´ y model vyhled´ av´ an´ı v datab´ az´ıch

V této kapitole si ukáˇzeme jednoduchou ilustraci toho, jak lze pomoc´ı vektor˚ u a jejich skalárn´ıch souˇcin˚ u prohledávat databáze. Nejprve si pˇripomeˇ nme, ˇze pro skalárn´ı souˇcin dvou aritmetických vektor˚ u x = [xi ] a y = [yi ] o n sloˇzkách plat´ı vztah xy = ||x|| ||y|| cos ϕ, kde ||x|| :=

√

xx =

q

x21 + x22 + · · · + x2n

je velikost vektoru x (analogicky definujeme pro y) a ϕ je u ´ hel, který sv´ıraj´ı vektory x a y (viz obr. 2.1).2 Nyn´ı si uvedeme jednu z moˇzných definic´ı bl´ızkosti vektor˚ u. O dvou aritmetických vektorech ˇrekneme, ˇze jsou si bl´ ızk´ e“, pokud je jejich u ´ hel ϕ bl´ızký ” 0, tj. hodnota cos ϕ je bl´ızká 1. Napˇr. na obr. 2.2 je vektor x bl´ızký“ y a ” vektor z nen´ı bl´ızký“ ani x ani y. ” 1

Matici I ˇr´ıkáme jednotkov´ a matice  1 0  0 1  I :=  . .  .. .. 0

a definujeme ji takto  0 ··· 0 0 0 ··· 0 0   .. . . . . . . .. ..  .

0 0

···

0 1

2´

Uhel ϕ tedy leˇz´ı v intervalu h0, πi.

29

´ Ulohy z line´ arn´ı algebry x

ϕ

y

´ Obrázek 2.1: Uhel ϕ mezi vektory x a y y x ϕ1 ϕ2 z

Obrázek 2.2: Bl´ızkost“ vektor˚ u ” Právˇe bl´ızkosti vektor˚ u m˚ uˇzeme vyuˇz´ıt pˇri vyhledáván´ı v databáz´ıch. Ilustrujme si to na následuj´ıc´ım pˇr´ıkladˇe. Mˇejme kuchaˇrskou knihu, která obsahuje m recept˚ u, kde kaˇzdý recept obsahuje nˇekterou z k surovin. Takovou databázi recept˚ u m˚ uˇzeme popsat m vektory r1 , r2 , . . . , rm o k sloˇzkách. Kaˇzdý vektor ri bude obsahovat jen nuly a jedniˇcky (0, resp. 1 na dané pozici v pˇr´ıpadˇe, ˇze danou ingredienci neobsahuje, resp. obsahuje). Chceme-li vybrat recept se surovinami, které nejlépe splˇ nuj´ı naˇse poˇzadovaná kritéria, pak m˚ uˇzeme náˇs poˇzadavek zapsat jako vektor p o k sloˇzkách a naˇs´ım u ´ kolem bude naj´ıt vektor ri z vektor˚ u r1 , r2 , . . . , rm , který je nejbl´ ıˇ ze vektoru p, tj. plat´ı, ˇze ri p cos ϕi = ||ri || ||p||

je nejvˇetˇs´ı ze vˇsech

cos ϕℓ =

rℓ p , ||rℓ || ||p||

ℓ = 1, 2 . . . , m.

Na obr. 2.3 ilustrujeme databázi ˇctyˇr recept˚ u, pˇriˇcemˇz danému poˇzadavku p na suroviny je nejbl´ıˇze recept r2 . Vektor p naˇseho poˇzadavku na suroviny mus´ı m´ıt stejný formát“ jako naˇse recepty. Zejména si uvˇedomme, ˇze 0 na ” dané pozici vektoru p znamená, ˇze poˇzadujeme, aby recept danou surovinu 30

p

r1

r2

r3

r4

Obrázek 2.3: Vektorový model recept˚ u a daný poˇzadavek neobsahoval. Tato skuteˇcnost nám ˇcasto m˚ uˇze pˇrekáˇzet, a proto je vhodné provést projekci“ recept˚ u (a vektoru p) jen do tˇech sloˇzek, ve kterých je p ” nenulový.3 Ukáˇzeme si to na pˇr´ıkladu jednoduché databáze tˇr´ı recept˚ u r1 = [1, 1, 1, 0, 0],

r2 = [1, 0, 1, 1, 1],

r3 = [1, 0, 0, 1, 1],

kde jednotlivé sloˇzky postupnˇe odpov´ıdaj´ı ingredienc´ım kuˇre, vepˇrové, ˇcesnek, cibule, mrkev. Nyn´ı nás bude zaj´ımat, které recepty obsahuj´ı kuˇre a cibuli. Náˇs poˇzadavek m˚ uˇzeme zapsat jako p = [1, 0, 0, 1, 0], t´ım ale ˇr´ıkáme, ˇze nechceme, aby recept obsahoval vepˇrové, ˇcesnek a mrkev. Ve skuteˇcnosti nás vˇsak obsah ostatn´ıch ingredienc´ı nezaj´ımá a provedeme proto projekci“ na prvn´ı a ˇctvrtou sloˇzku4 ” r1 = [1, 0], r2 = [1, 1], r3 = [1, 1], p = [1, 1] a pomoc´ı kos´ın˚ uu ´ hl˚ u mezi projektovaným“ p a projektovanými“ recepty ” ” zjist´ıme, ˇze naˇsemu poˇzadavku (na obsah kuˇrete a cibule) jsou nebl´ıˇze recepty 2 a 3. Pˇ r´ıklad 1 Popiˇste následuj´ıc´ı recepty pomoc´ı aritmetických vektor˚ u: 1. Buchty - ingredience: 100 g Hery, 100 g cukru mouˇ cka, 2 vejce, 500 g hladk´ e mouky, 1/4 l ml´ eka, 30 g droˇ zd´ı ˇ ınsk´ 2. C´ e placiˇ cky - ingredience: 3 silnˇ ejˇs´ı kuˇrec´ı ˇr´ızky, 3 vejce, 1 drobnˇ e nakr´ ajen´ a cibule, 4 lˇ z´ıce Solamylu, s˚ ul, 4 lˇ z´ıce oleje, 1 lˇ z´ıce sojov´ e om´ aˇ cky, 1 lˇ z´ıce worcesterov´ e om´ aˇ cky ˇ ckov´ 3. Coˇ a pol´ evka - ingredience: 1 lˇ z´ıce olivov´ eho oleje, 1 mrkev, 1 cibule, 4 strouˇ zky ˇ cesneku, 50 g ˇ zitn´ e mouky, 450 g ˇ coˇ cky, s˚ ul, 1 kostka zeleninov´ eho buj´ onu, podle chuti chilli nebo cayensk´ y pepˇr, major´ anka, 10 g m´ asla, 2 lˇ z´ıce plnotuˇ cn´ e hoˇrˇ cice, 2 vejce, petrˇ zelka 3 4

To, zda projekci“ provedeme ˇci ne, tedy m˚ uˇze ovlivnit v´ ysledek vyhledáván´ı. ” Tedy na ty sloˇzky, ve kter´ ych je poˇzadavek p nenulov´ y.

31

´ Ulohy z line´ arn´ı algebry 4. Dom´ ac´ı buchty - ingredience: 600 g polohrub´ e mouky, 1/4 l ml´ eka, 2 vejce, 10 lˇ zic oleje, 40 g cukru, 1 pr´ aˇsek do peˇ civa, 1 kostiˇ cka droˇ zd´ı (42 g), ˇspetka soli 5. Ev´ıkova mˇ namka - ingredience: 4 kuˇrec´ı ˇr´ızky, 8 pl´ atk˚ u anglick´ e slaniny, s´ yr podle chuti (eidam, blat’´ ack´ e zlato, hermel´ın, niva..., ale vˇ zdy jen jeden druh), 1 cibule, s˚ ul, pepˇr, 50 g m´ asla, 1 lˇ z´ıce sojov´ e om´ aˇ cky 6. Chl´ eb - ingredience: 1,5 lˇ z´ıce octa, 3 lˇ z´ıce olivov´ eho oleje, 10 g cukru, 3 lˇ ziˇ cky soli, 360 g hladk´ e mouky pˇseniˇ cn´ e, 140 g ˇ zitn´ e mouky, 75 g celozrnn´ e mouky ˇ zitn´ e, 75 g celozrnn´ e mouky pˇseniˇ cn´ e, lˇ z´ıce km´ınu, 15 g suˇsen´ eho droˇ zd´ı 7. Chleb´ıˇ ckov´ a pomaz´ anka - ingredience: 100 g brambor, 1 cibule, 1 lˇ z´ıce tatarsk´ e om´ aˇ cky, s˚ ul, pepˇr 8. Kokosov´ a hrn´ıˇ ckov´ a b´ abovka - ingredience: 200 g hladk´ e mouky, 10 lˇ zic kokosu, 100 g cukru krupice, 1/4 l ml´ eka, 6 lˇ zic oleje, 1 pr´ aˇsek do peˇ civa, 1 vanilkov´ y cukr, 3 vejce 9. Kuˇrec´ı kousky v s´ yrov´ em tˇ est´ıˇ cku - ingredience: 4 kuˇrec´ıch prs´ıˇ cek, s˚ ul Tˇ est´ıˇ cko: 1 vejce, 0,05 l b´ıl´ eho v´ına, 80 g hladk´ e mouky, strouhan´ y s´ yr Dresink: 100 g b´ıl´ eho jogurtu, 5 lˇ zic tatarky, mlet´ y b´ıl´ y pepˇr, s˚ ul, 2 strouˇ zky ˇ cesneku 10. Paˇr´ıˇ zsk´ e kostky - ingredience: 440 g polohrub´ e mouky, 220 g cukru, 7 lˇ zic oleje, 1/4 l vlaˇ zn´ eho ml´ eka, 2-3 lˇ z´ıce kakaa, 1 pr´ aˇsek do peˇ civa, 3 cel´ a vejce 11. Pern´ık - ingredience: 1/2 kg polohrub´ e mouky, 350 g cukru krystal, 1/2 l ml´ eka, 10 lˇ zic oleje, 2 cel´ a vejce, 4 lˇ z´ıce rozˇredˇ en´ ych povidel, 2-3 lˇ z´ıce kakaa, 1 lˇ ziˇ cka jedl´ e sody, 1 pr´ aˇsek do peˇ civa, ˇspetka soli, 1 lˇ ziˇ cka mlet´ e skoˇrice, lze pˇridat najemno nastrouhan´ a 2-3 jablka 12. Pizza tˇ esto - ingredience: 0,5 kg hladk´ e mouky, 1 lˇ z´ıce olivov´ eho oleje, 1 lˇ z´ıce soli, 15 g droˇ zd´ı 13. Plnˇ en´ e kuˇre - ingredience: 1 kuˇre, 1 velk´ a cibule, 3 pl´ atky anglick´ e slaniny, 2 lˇ z´ıce oleje, s˚ ul 14. Plnˇ en´ y cop - ingredience: 500 g polohrub´ e mouky, 50 g mouˇ ckov´ eho cukru, 3 ˇ zloutky, 100 g rozpuˇstˇ en´ eho m´ asla, ˇspetka soli, 0,2 l ml´ eka, 40 g droˇ zd´ı N´ aplˇ n: 5 lˇ zic strouhan´ ych l´ıskov´ ych oˇr´ıˇsk˚ u, 100 g cukru, 3 vejce, ˇspetka skoˇrice, vejce na potˇren´ı, sekan´ e oˇr´ıˇsky na posyp´ an´ı 15. Tatransk´ e pracny - ingredience: 300 g hladk´ e mouky, 250 g Hery, 100 g mouˇ ckov´ eho cukru, 1 vejce, 1 lˇ z´ıce kakaa, 1 lˇ ziˇ cka skoˇrice, 6 lˇ zic mlet´ ych oˇrech˚ u, oˇr´ıˇsk˚ u nebo mandl´ı (m˚ uˇ ze se i nam´ıchat) 16. Tradiˇ cn´ı italsk´ e lasagne - ingredience: Boloˇ nsk´ a om´ aˇ cka: 1 lˇ z´ıce olivov´ eho oleje, 1 stˇredn´ı cibule, 1 strouˇ zek ˇ cesneku, 500 g mlet´ eho hovˇ ez´ıho (kuˇrec´ıho, s´ ojov´ eho masa, jak´ e kdo m´ a r´ ad), 250 g oloupan´ ych cel´ ych ˇ ci nakr´ ajen´ ych rajˇ cat v plechovce, 2 lˇ z´ıce rajˇ catov´ eho protlaku, 4 lˇ z´ıce ˇ cerven´ eho v´ına, s˚ ul, ˇ cerstvˇ e namlet´ y pepˇr S´ yrov´ a om´ aˇ cka: 50 g m´ asla, 50 g mouky, 0,6 l ml´ eka, s´ yr na strouh´ an´ı (ˇ cedar), 12 list˚ u vajeˇ cn´ ych tˇ estovin - lasagn´ı 17. Vepˇrov´ a k´ yta - ingredience: 6 ˇr´ızk˚ u z vepˇrov´ e k´ yty, 1 lˇ ziˇ cka soli, 1 lˇ ziˇ cka pepˇre, 1/2 lˇ z´ıce solamylu, 1 lˇ z´ıce worchestrov´ e om´ aˇ cky, 2 strouˇ zky ˇ cesneku, 1/2 lˇ ziˇ cky major´ anky, 2 vejce, 1 lˇ z´ıce raj. protlaku, 1 lˇ z´ıce oleje

Nyn´ı zkuste na základˇe výˇse popsané bl´ızkosti vektor˚ u vyhledat, které recepty obsahuj´ı ingredience nejlépe odpov´ıdaj´ıc´ı tˇemto poˇzadavk˚ um: 1) cukr, droˇzd´ı, Hera, kakao, mouka a skoˇrice, 2) s˚ ul, sýr a worchester, 3) ˇcesnek, hoˇrˇcice, kuˇrec´ı maso, majoránka, pepˇr, slanina a s˚ ul.

32

Pˇ r´ıklad 2 Vyhledejte recepty z naˇs´ı databáze, které neobsahuj´ı 1) mléko, jogurt, máslo, sýr ani tatarskou omáˇcku (recepty pro pacienty s alergi´ı na b´ılkovinu kravského mléka), 2) ˇzádné maso (hovˇez´ı, kuˇrec´ı, vepˇrové) ani slaninu (vegetariánské recepty).

2.3

Geometrick´ e transformace

Nejdˇr´ıve si pˇripomeneme nˇekteré pojmy z oblasti geometrických zobrazen´ı. Geometrick´ ym zobrazen´ ım rozum´ıme zobrazen´ı, které geometrickému u ´ t´ varu U pˇriˇrad´ı geometrický u ´ tvar U ′ . Utvar U je tzv. vzor a u ´ tvar U ′ se oznaˇcuje jako obraz. Pod´ıvejme se bl´ıˇze na tˇri konkrétn´ı typy zobrazen´ı, a to stejnolehlost, rotaci a posunut´ı. Stejnolehlost Uvaˇzujme bod S ∈ Rn , kde n ∈ {2, 3}. Geometrické zobrazen´ı, pˇri nˇemˇz obrazem bodu S je bod S a obrazem kaˇzdého bodu A ∈ Rn , A 6= S, je takový bod A′ ∈ Rn , ˇze pro vektor SA′ plat´ı SA′ = κSA, kde κ ∈ R \ {0, 1} je pevnˇe zvolené, se nazývá stejnolehlost´ ı (nebo také homoteti´ ı). Bod S nazýváme stˇ redem stejnolehosti a κ koeficientem stejnolehlosti. Stejnolehlost s κ = −1 je stˇredovou soumˇernost´ı. 5

4

y

3

2

1

0 0

1

2

3

4

5

x

Obrázek 2.4: Stejnolehlost (v rovinˇe) s koeficientem κ = 2 a stˇredem v poˇcátku (vzor zelenˇe, obraz modˇre) 33

´ Ulohy z line´ arn´ı algebry Rotace Rotace (otoˇ cen´ ı) v rovinˇe je geometrické zobrazen´ı, které je charakterizováno t´ım, ˇze spojnice vˇsech bod˚ u s pevnˇe zvoleným bodem S, tzv. stˇ redem otoˇ cen´ ı, se zmˇen´ı o stejný u ´ hel ϕ, tzv. ´ uhel otoˇ cen´ ı, a vzdálenost bod˚ u od stˇredu otoˇcen´ı z˚ ustává nezmˇenˇena.

3

y

2

1

0 0

1

2

3

x

Obrázek 2.5: Rotace s u ´ hlem ϕ = π/4 a stˇredem v poˇcátku (vzor zelenˇe, obraz modˇre) Posunut´ı Posunut´ ı (translace) je geometrické zobrazen´ı, které je charakterizováno t´ım, ˇze vˇsechny body transformované mnoˇziny bod˚ u zmˇen´ı své kartézské souˇradnice o stejnou hodnotu, tj. ke kaˇzdému bodu pˇriˇcteme stejný vektor posunut´ ı. Nyn´ı si ukáˇzeme, jak lze pomoc´ı elementárn´ıch maticových operac´ı zapsat výˇse uvedená zobrazen´ı. Zapiˇsme souˇradnice vzoru geometrického zobrazen´ı, tj. bodu z R2 ˇci R3 , do sloupcového vektoru. Pokud vzorem zobrazen´ı je v´ıce bod˚ u, zap´ıˇseme je jako sloupce matice P . Obraz bod˚ u ve stejnolehlosti s koeficientem κ a stˇredem v poˇcátku zap´ıˇseme jako souˇcin transformaˇcn´ı matice T typu (n, n) a matice P , pˇriˇcemˇz n je dimenze prostoru, ve kterém zobrazujeme (tj. 2 nebo 3). Matice T má vˇsechny prvky na hlavn´ı diagonále rovny koeficientu stejnolehlosti κ. Vˇsechny dalˇs´ı prvky jsou nulové.5 Pro n = 2 5

34

Tedy T = κI, kde I je jednotková matice.

4

y

3

2

1

0 0

1

2

3

4

x

Obrázek 2.6: Posunut´ı (v rovinˇe) o vektor [1, 1] (vzor zelenˇe, obraz modˇre) zap´ıˇseme obrazy bod˚ u v rotaci s u ´ hlem otoˇcen´ı ϕ a stˇredem otoˇcen´ı v poˇcátku jako souˇcin transformaˇcn´ı matice R a matice P , pˇriˇcemˇz cos ϕ − sin ϕ R= . sin ϕ cos ϕ Posunut´ı s vektorem posunut´ı p realizujeme tak, ˇze ke kaˇzdému sloupci matice P pˇriˇcteme sloupcový vektor p. Pˇ r´ıklad 1 Implementujte geometrická zobrazen´ı stejnolehlost, rotaci a posunut´ı v rovinˇe. Naleznˇete obraz troj´ uheln´ıku s vrcholy [1, −1], [2, 0] a [1, 1] ve stejnolehlosti se stˇredem v poˇcátku a s koeficientem stejnolehlosti κ = 2. Tento obraz transformujte rotac´ı se stˇredem otoˇcen´ı v poˇcátku a u ´ hlem otoˇcen´ı ϕ = π/2. Nakonec na tento nový obraz aplikujte posunut´ı s vektorem posunut´ı p = [−3, −2]. Zobrazte vzor (tj. p˚ uvodn´ı troj´ uheln´ık) a jednotlivé obrazy (pouˇzijte napˇr. matlabovskou funkci fill). Pˇ r´ıklad 2 Vytvoˇrte funkci s názvem transformace, která transformuje zadaný bod podle vstupn´ıch parametr˚ u: transformovany bod = transformace(bod, κ, ϕ, p). Transformaci budeme chápat ve smyslu pˇr´ıkladu 1, tj. nejdˇr´ıve aplikujeme stejnolehlost se stˇredem v poˇcátku a s koeficientem κ, poté rotaci se stˇredem v poˇcátku a u ´ hlem otoˇcen´ı ϕ a nakonec posunut´ı s vektorem posunut´ı p. Tuto funkci pouˇzijte pro vytvoˇren´ı rutiny, která transformuje n-´ uheln´ık zadaný matic´ı P (sloupce matice jsou tvoˇreny vektory souˇradnic vrchol˚ u).

35

´ Ulohy z line´ arn´ı algebry Pˇ r´ıklad 3 Vytvoˇrte rutinu, která v kaˇzdé iteraci s danou pravdˇepodobnost´ı provede transformaci bodu z pˇredchoz´ı iterace podle tabulky 2.1. Transformovaný bod v kaˇzdé iteraci vykreslete. Jako poˇcáteˇcn´ı bod volte [0,5, 0,5]. (Tento pˇr´ıklad byl pˇrevzat a upraven z [3].) pravdˇepodobnost 0,85 0,07 0,07 0,01

κ (stejnolehlost) 0,85 0,3 0,25

ϕ (rotace) −0,05 1 −1

p (posunut´ı) [0, 1,5] [0, 1,5] [0, 0,5] 0 0 transformace definovaná matic´ı 0 0,15

Tabulka 2.1: Krok iteraˇcnˇe definovaného geometrického zobrazen´ı

36

Kapitola 3 Poˇ c´ ateˇ cn´ı u ´lohy 3.1

Struˇ cn´ y u ´ vod do poˇ c´ ateˇ cn´ıch u ´ loh a jak je ˇ reˇsit

Na u ´ vod se nauˇc´ıme pracovat se speciáln´ım typem rovnice – s tzv. diferencia´ln´ı rovnic´ı. Pomoc´ı diferenciáln´ıch rovnic lze popsat celou ˇradu zákonitost´ı, které se objevuj´ı v pˇr´ırodn´ıch a spoleˇcenských vˇedách. Podrobnˇe se tomuto typu rovnic vˇenuje text [8]. V této ˇcásti si také ukáˇzeme, jak lze diferenciáln´ı rovnice numericky ˇreˇsit. A v následuj´ıc´ıch podkapitolách budeme pomoc´ı diferenciáln´ıch rovnic konstruovat modely popisuj´ıc´ı nˇekteré jevy, se kterými se bˇeˇznˇe setkáváme. Na závˇer poznamenejme, ˇze mnoho dalˇs´ıch zaj´ımavých aplikac´ı obyˇcejných diferenciáln´ıch rovnic lze nalézt v [8, 9]. Obyˇ cejnou diferenci´ aln´ ı rovnic´ ı 1. ˇ r´ adu rozum´ıme rovnici tvaru y ′ (t) = f (t, y(t)),

(3.1)

kde f: R2 7→ R je zadaná funkce a y je hledaná“ funkce. ” ˇ Reˇ sen´ ım t´ eto rovnice na otevˇ ren´ em intervalu (a, b) (a < b) rozum´ıme kaˇzdou funkci y¯: (a, b) 7→ R takovou, ˇze pro vˇsechna t ∈ (a, b) plat´ı y¯′ (t) = f (t, y¯(t)). Napˇr´ıklad funkce y¯(t) := t, t ∈ (0, +∞), je ˇreˇsen´ım diferenciáln´ı rovnice y ′(t) = y(t) na intervalu (0, +∞). Jiným ˇreˇsen´ım této rovnice na intervalu t (0, +∞) je napˇr´ıklad funkce y¯(t) := 2t, t ∈ (0, +∞), ˇci y¯(t) := 3t, t ∈ (0, +∞). Nen´ı tˇeˇzké ukázat, ˇze pro libovolné k ∈ R je funkce y¯k (t) := kt, t ∈ (0, +∞), ˇreˇsen´ım diferenciáln´ı rovnice y ′(t) = y(t) na intervalu (0, +∞). t Naˇse u ´ loha má nekoneˇcnˇe mnoho ˇreˇsen´ı. Zkusme nav´ıc pˇridat k naˇs´ı rovnici 37

Poˇc´ ateˇcn´ı u ´lohy napˇr´ıklad podm´ınku y(1) = 2, tj. chceme nalézt funkci, která ˇreˇs´ı naˇsi rovnici a nav´ıc jej´ı funkˇcn´ı hodnota v t = 1 je rovna 2. Lze ukázat, ˇze taková u ´ loha má na (0, +∞) pouze jediné ˇreˇsen´ı y¯(t) = 2t. ´ Ulohu, která se skládá z hledán´ı ˇreˇsen´ı obyˇcejné diferenciáln´ı rovnice 1. ˇrádu, které má nav´ıc splˇ novat tzv. poˇcáteˇcn´ı podm´ınku y(t0 ) = y0 , nazýváme Cauchyovou ´ ulohou1 a zapisujeme ji obecnˇe takto: ) y ′(t) = f (t, y(t)), (3.2) y(t0 ) = y0 . Pokud má funkce f (t, y(t)) speciáln´ı tvar,2 pak existuj´ı metody, jak anaˇ lyticky ˇreˇsit výˇse popsanou Cauchyovu u ´ lohu. Casto vˇsak analytické ˇreˇsen´ı nalézt nelze nebo by jeho nalezen´ı bylo pˇr´ıliˇs nároˇcné. V takovém pˇr´ıpadˇe se nab´ız´ı pouˇzit´ı nˇekteré z numerických metod pro pˇribliˇzné ˇreˇsen´ı Cauchyových u ´ loh. Pod´ıvejme se nyn´ı na jednu z tˇechto metod – Eulerovu metodu.3 O funkci f budeme dále pˇredpokládat, ˇze je v mnoˇzinˇe D := {(t, y) ∈ R2 : t0 ≤ t} spojitá a také jsou v D spojité derivace této funkce podle promˇenné y. Eulerova metoda Eulerova metoda je nejjednoduˇsˇs´ım zp˚ usobem numerického ˇreˇsen´ı Cauchyových u ´ loh. Výstup Eulerovy metody nám aproximuje ˇreˇsen´ı Cauchyovy u ´ lohy (3.2) na intervalu ht0 , tN i (tento interval specifikujeme n´ıˇze). Metoda vyuˇz´ıvá aproximace derivace4 funkce y v bodˇe t pomoc´ı tzv. dopˇ redn´ e diference y v tomto bodˇe y ′(t) ≈

y(t + h) − y(t) , h

kde h je malé“ a kladné. ” Po jednoduché u ´ pravˇe dostaneme y(t + h) ≈ y(t) + hy ′ (t).

Nahrad´ıme-li funkci y ′(t) funkc´ı f (t, y(t)) (viz diferenciáln´ı rovnici u ´ lohy (3.2)), z´ıskáme vztah y(t + h) ≈ y(t) + hf (t, y(t)). 1

(3.3)

Pojmenována podle francouzského matematika Augustina Louise Cauchyho (1789– 1857). 2 Napˇr´ıklad f závis´ı lineárnˇe na y, tj. f (t, y) = a(t)y + b(t), kde a a b jsou reálné funkce. 3 Publikoval ji v´ yznamn´ y ˇsv´ ycarsk´ y matematik a fyzik Leonhard Euler (1707–1783) v roce 1768. 4 Pro jistotu pˇripomeˇ nme, ˇze derivace funkce y v bodˇe t je definována jako lim y(t+h)−y(t) . h

h→0

38

Dále zvolme dostateˇcnˇe malé“ pevné h > 0 a sestrojme posloupnost ” t0 , t1 := t0 + h, t2 := t0 + 2h, . . . , tN := t0 + Nh. Oznaˇcme pomoc´ı yn aproximaci hodnoty pˇresného ˇreˇsen´ı y(tn ). Po dodán´ı poˇcáteˇcn´ı podm´ınky dostáváme z (3.3) rekurzivn´ı vztah y0 = y(t0 ), (3.4) yn+1 = yn + hf (tn , yn ), n = 0, 1, . . . , N − 1, který pouˇzijeme pro numerické ˇreˇsen´ı Cauchyovy u ´ lohy (3.2). Dá se ukázat, ˇze odhad chyby aproximace ˇreˇsen´ı Cauchyovy u ´ lohy pomoc´ı Eulerovy metody v bodˇe t1 je pˇr´ımo u ´ mˇerný druhé mocninˇe kroku h. Odhad chyby aproximace ˇreˇsen´ı Cauchyovy u ´ lohy v bodˇe tN je pˇr´ımo u ´ mˇerný kroku h. Plat´ı totiˇz, ˇze absolutn´ı hodnotu chyby aproximace ˇreˇsen´ı Cauchyovy u ´ lohy v bodˇe tN m˚ uˇzeme odhadnout shora výrazem M K(tN −t0 ) e −1 h 2K

(3.5)

(tento odhad je moˇzné pouˇz´ıt jen pˇri splnˇen´ı podm´ınky Kh < 1), kde ∂f (t, y) K := max t∈ht0 ,tN i, y∈hy0 ,yN i ∂y a

∂f 5 ∂f . M := max + f (t, y) t∈ht0 ,tN i, y∈hy0 ,yN i ∂t ∂y

3.2

Modelov´ an´ı zmˇ eny koncentrace l´ ek˚ u v krvi

V této ˇcásti se pokus´ıme vytvoˇrit matematický model, který popisuje, jakým zp˚ usobem se mˇen´ı koncentrace léku v krvi v závislosti na ˇcase. Necht’ funkce c(t) ≥ 0 udává okamˇzitou koncentraci látky (vhodného léku) v krvi v ˇcase t (v µg/mℓ). Lékaˇrskými pokusy bylo zjiˇstˇeno, ˇze pokles koncentrace této látky v krvi je pˇr´ımo u ´ mˇerný jej´ı samotné koncentraci, tj. plat´ı: c′ (t) = −kc(t), 6 ∂f Symboly ∂f cuj´ı parciáln´ı derivaci funkce f podle promˇenné y, resp. t. ∂y a ∂t oznaˇ O parciáln´ıch derivac´ıch se v´ıce dozv´ıte v kapitole 5.1. 6 Uvˇedomme si, ˇze koncentrace léku v krvi klesá. Proto zmˇena koncentrace léku popsan´ a funkc´ı c′ je záporná. 5

39

Poˇc´ ateˇcn´ı u ´lohy kde k > 0 je konstanta.7 Tato konstanta popisuj´ıc´ı u ´ bytek dané látky je urˇcena dvˇema farmakokinetickými parametry, kterými jsou clearance (m´ıra schopnosti organismu eliminovat látku) a distribuˇ cn´ ı objem (m´ıra kapacity zdánlivého prostoru, který je v organismu pro tuto látku k dispozici). Clearance budeme dále znaˇcit cℓ a jej´ı jednotky budou mℓ/min, distribuˇcn´ı objem oznaˇc´ıme vd a jeho jednotky budou ℓ. Pokud zvol´ıme jako jednotku ˇcasu hodiny, m˚ uˇzeme konstantu k popsat následuj´ıc´ı závislost´ı k=

60 cℓ 1000 . vd

Dále pro jednoduchost pˇredpokládejme, ˇze látka je distribuována do krve intravenózn´ı injekc´ı a rozˇsiˇruje se do krve okamˇzitˇe. Pˇredpokládejme, ˇze t´ımto zp˚ usobem bylo v ˇcase t = 0 dodáno do krve takové mnoˇzstv´ı látky, ˇze jej´ı koncentrace v krvi mˇela hodnotu c0 . T´ım jsme z´ıskali jednoduchý model, popisuj´ıc´ı hodnoty koncentrace látky v krvi po jej´ı intravenózn´ı aplikaci: ) c′ (t) = −kc(t), (3.6) c(0) = c0 . U léˇciv je dalˇs´ım významným farmakokinetickým parametrem ´ uˇ cinn´ a koncentrace. Ta udává hodnotu koncentrace ˇci interval hodnot koncentrac´ı, pˇri kterých látka p˚ usob´ı prospˇeˇsnˇe na organismus. Pokud známe hodnotu výˇse zm´ınˇených farmakokinetických parametr˚ u pro konkrétn´ı léˇcivo, m˚ uˇzeme vyuˇz´ıt ˇreˇsen´ı u ´ lohy (3.6) pro odpovˇed’ na otázku, jak ˇcasto je tˇreba lék obsahuj´ıc´ı tuto látku pacientovi aplikovat pro zajiˇstˇen´ı u ´ spˇeˇsné léˇcby. V následuj´ıc´ıch pˇr´ıkladech pˇredpokládáme, ˇze pacientem je pr˚ umˇerná osoba s tˇelesnou hmotnost´ı 70 kg. Hodnoty farmakokinetických parametr˚ u pro léˇciva z tˇechto pˇr´ıklad˚ u byly pˇrevzaty z [6], kde se lze seznámit s oblast´ı farmakokinetiky daleko detailnˇeji. Pˇ r´ıklad 1 Pro léˇcbu nemocného s pr˚ uduˇskovým astmatem se pouˇz´ıvá theofylin. Jeho farmakokinetické parametry jsou: cℓ = 48 mℓ/min, vd = 35 ℓ, u ´ˇcinná koncentrace se pohybuje v intervalu h10, 20i µg/mℓ. Toto léˇcivo mus´ıme pacientovi podávat v pravidelných ˇcasových intervalech tak, aby jeho koncentrace v krvi léˇcené osoby nepˇresáhla horn´ı mez u ´ˇcinné koncentrace a neklesla pod jej´ı doln´ı mez. Na zaˇcátku léˇcby byla aplikována 7

40

Hodnota k závis´ı na léku a pacientovi.

Obrázek 3.1: Molekula theofylinu zavádˇec´ı dávka, která zp˚ usobila, ˇze v ˇcase t = 0 byla koncentrace theofylinu v krvi c0 = 20µg/mℓ. Dalˇs´ı dávky chceme aplikovat vzhledem k pohodl´ı pacienta tak, aby interval podáván´ı léku byl co nejvˇetˇs´ı. Zjistˇete, po jakém ˇcase T (v hodinách) je nutné lék obsahuj´ıc´ı theofylin znovu podat pacientovi. Pˇripomeˇ nme, ˇze aby byla léˇcba u ´ˇcinná, je tˇreba lék pacientovi podat dˇr´ıve, neˇz koncentrace theofylinu klesne pod doln´ı mez u ´ˇcinné koncentrace. Odhadnˇete chybu aproximace koncentrace léku v ˇcase T . Pˇ r´ıklad 2 Pro léˇcbu horeˇcnatých stav˚ u, bolesti hlavy, sval˚ u ˇci kloub˚ u se pouˇz´ıvá kyselina acetylsalicylová.8 Jej´ı farmakokinetické parametry jsou:

O

OH O O

Obrázek 3.2: Molekula kyseliny acetylsalicylové cℓ = 650 mℓ/min, vd = 11 ℓ, u ´ˇcinná koncentrace se pohybuje v intervalu h150, 300i µg/mℓ. Toto léˇcivo mus´ıme pacientovi podávat v pravidelných ˇcasových intervalech tak, aby jeho koncentrace v krvi léˇcené osoby nepˇresáhla horn´ı mez u ´ˇcinné koncentrace a neklesla pod jej´ı doln´ı mez. Na zaˇcátku léˇcby byla aplikována 8

Kyselina acetylsalicylová je hlavn´ı sloˇzkou lék˚ u jako je Aspirin, Acylpirin ˇci Anopyrin.

41

Poˇc´ ateˇcn´ı u ´lohy zavádˇec´ı dávka, která zp˚ usobila, ˇze v ˇcase t = 0 byla koncentrace kyseliny acetylsalicylové v krvi c0 = 300µg/mℓ. Dalˇs´ı dávky chceme aplikovat vzhledem k pohodl´ı pacienta tak, aby interval podáván´ı léku byl co nejvˇetˇs´ı. Zjistˇete, po jakém ˇcase T (v hodinách) je nutné lék obsahuj´ıc´ı kyselinu acetylsalicylovou znovu podat pacientovi. Aby byla léˇcba u ´ˇcinná, je tˇreba lék pacientovi podat dˇr´ıve, neˇz koncentrace kyseliny acetylsalicylové klesne pod doln´ı mez u ´ˇcinné koncentrace. Odhadnˇete chybu aproximace koncentrace léku v ˇcase T .

3.3

Tak je to padˇ elek nebo to nen´ı padˇ elek aneb jak poznat st´ aˇ r´ı nˇ ekter´ ych Vermee” rov´ ych“ obraz˚ u?

V této ˇca´sti se pokus´ıme vytvoˇrit matematický model, který popisuje rozpad atom˚ u radioaktivn´ıch prvk˚ u. Neˇz se zaˇcneme zabývat tvorbou dalˇs´ıho matematického modelu, vrát´ıme se v ˇcase do doby krátce po konci druhé svˇetové války. Tˇesnˇe po válce zjistila nizozemská policie, ˇze bˇehem války bylo prodáno nˇekolik Vermeerových9 obraz˚ u nˇemeckému ministrovi letectv´ı Hermannu Göringovi. Tuto transakci zprostˇredkoval nizozemský mal´ıˇr Han van Meegeren. Na základˇe tˇechto zjiˇstˇených fakt˚ u byl 29. 5. 1945 van Meegeren zadrˇzen a obvinˇen z kolaborace s nepˇr´ıtelem. 12. 7. 1945 van Meegeren vydal prohláˇsen´ı, ˇze Göringovi nikdy ˇzádný Vermeer˚ uv obraz neprodal. Naopak Göringa napálil, protoˇze obrazy, které mu prodal, jsou podvrhy Vermeerových obraz˚ u a sám je vytvoˇril. A aby dokázal své tvrzen´ı, zaˇcal jeden z Vermeerových“ obraz˚ u10 napo” dobovat. Van Meegeren pˇrizvaným znalc˚ um pˇredvedl zp˚ usob, jakým vytváˇr´ı barvy, jak pˇripravuje plátno, ˇci jak zaˇr´ıd´ı, aby povrch malby vypadal jako u nˇekolik set let starého obrazu. Tˇesnˇe pˇred dokonˇcen´ım podvrhu Vermeerova obrazu se van Meegeren dozvˇedˇel, ˇze obvinˇen´ı z kolaborace bude nahrazeno obvinˇen´ım z padˇelatelstv´ı, a tak odm´ıtl tuto kopii dokonˇcit. I tak ale vˇetˇsina pˇrizvaných odborn´ık˚ u uznala, ˇze obrazy prodané Göringovi jsou pravdˇepodobnˇe falzum a van Meegeren byl 12. 11. 1947 odsouzen za padˇelatelstv´ı na rok do vˇezen´ı, ve kterém 30. 12. 1947 na infarkt zemˇrel. I pˇresto, ˇze komise, která posuzovala pravost Vermeerových“ obraz˚ u ” uznala, ˇze to jsou pravdˇepodobnˇe podvrhy vytvoˇrené van Meegerenem, z˚ ustávali odborn´ıci u nˇekterých obraz˚ u, k jejichˇz autorstv´ı se také van Meegeren pˇrihlásil, na pochybách. Zejména zpochybˇ nován´ı pravosti obrazu Emauzˇst´ı uˇcedn´ıci, který zakoupilo muzeum v Rotterdamu za 170 000 dolar˚ u, vy9 10

42

Jan Vermeer (1632–1675) byl nizozemsk´ y mal´ıˇr. Konkrétnˇe ˇslo o obraz Jeˇz´ıˇs mezi znalci P´ısma.

volávalo velké spory. Proto se pˇristoupilo u tohoto obrazu v roce 1967 k metodˇe radioaktivn´ıho datován´ı, která mˇela tyto pochyby rozhodnout. Metoda radioaktivn´ıho datován´ı vyuˇz´ıvá toho, ˇze nˇekteré tzv. radioaktivn´ı prvky jsou nestabiln´ı a ˇcást jejich atom˚ u se samovolnˇe rozpadá na atomy jiných prvk˚ u. Experimenty bylo zjiˇstˇeno, ˇze rychlost rozpadu atom˚ u radioaktivn´ıch prvk˚ u je pˇr´ımo u ´ mˇerná poˇctu tˇechto atom˚ u. Pokud funkci udávaj´ıc´ı poˇcet atom˚ u radioaktivn´ıho prvku v ˇcase t v gramu látky oznaˇc´ıme jako N(t), pak výˇse zm´ınˇenou závislost m˚ uˇzeme popsat diferenciáln´ı rovnic´ı N ′ (t) = −λN(t), 11

(3.7)

kde λ je konstanta, která popisuje rychlost rozpadu atom˚ u daného radioaktivn´ıho prvku. Tato konstanta je dána pro kaˇzdý radioaktivn´ı prvek t´ımto vztahem ln 2 λ= . poloˇcas rozpadu prvku v roc´ıch ˇ t v naˇsem modelu budeme mˇeˇrit v roc´ıch a jednotkou konstanty λ je Cas rok −1 . Metoda radioaktivn´ıho datován´ı je zaloˇzena na jednoduchém pozorován´ı. Pokud bychom vˇedˇeli, kolik atom˚ u radioaktivn´ıho prvku mˇela látka v jednom svém gramu pˇri svém vzniku (tzn. známe hodnotu N0 , pro kterou plat´ı N(0) = N0 ), a znali bychom také aktuáln´ı poˇcet tˇechto atom˚ u v gramu látky, mohli bychom ˇreˇsen´ım u ´ lohy ) N ′ (t) = −λN(t), (3.8) N(0) = N0 zjistit, jak je tato látka stará. Neˇz se zaˇcneme zabývat datován´ım Vermeerových“ obraz˚ u, uvˇedomme ” si, ˇze vˇsechny horniny na Zemi obsahuj´ı malé mnoˇzstv´ı radioaktivn´ıho uranu, který se rozpadá na atomy dalˇs´ıho prvku. Tyto atomy se opˇet samovolnˇe mˇen´ı na dalˇs´ı atomy atd. Viz obr. 3.3 (ˇcasy u ˇsipek udávaj´ı poloˇcasy rozpadu12 jednotlivých radioaktivn´ıch prvk˚ u). Dále je známo, ˇze olovnatá bˇeloba pouˇz´ıvaná na malbách obsahuje oxid olovnatý, který obsahuje malé mnoˇzstv´ı olova-210 a jeˇstˇe menˇs´ı mnoˇzstv´ı radia-226. V okamˇziku, kdy je barva obsahuj´ıc´ı oxid olovnatý vyrobena, zaˇcnou se atomy olova-210 velmi rychle rozpadat s poloˇcasem rozpadu 22 11

Uvˇedomme si, ˇze poˇcet atom˚ u radioktivného prvk˚ u klesá v d˚ usledku jejich samovolného rozpadu. Proto je zmˇena poˇctu atom˚ u popsaná funkc´ı N ′ záporná. 12 Poloˇcas rozpadu je doba, za kterou poˇcet atomov´ ych jader ve vzorku klesne na polovinu.

43

Poˇc´ ateˇcn´ı u ´lohy Uran - 234

Uran - 238 1.2 minuty 4.5 miliardy let

Thorium - 234

250 tisíc let Protaktinium - 234

24 dní

Thorium - 230

80 tisíc let

Radium - 226

1 600 let

Radon - 222

3.8 dne

Polonium - 218

3 minuty

Olovo - 214

20 minut Bismut - 214

27 minut

Polonium - 214

5 dní

necelá sekunda Olovo - 210

Bismut - 210

22 let

Polonium - 210

138 dní

Olovo - 206

Obrázek 3.3: Uranová rozpadová ˇrada let a mnoˇzstv´ı olova-210 v této barvˇe klesá. Na druhé stranˇe vzniká malé mnoˇzstv´ı olova-210 rozpadem radia-226 (a prvk˚ u, které následuj´ı v rozpadové ˇradˇe za n´ım). Tento proces m˚ uˇzeme popsat následuj´ıc´ı poˇcáteˇcn´ı u ´ lohou ) N ′ (t) = −λN(t) + r(t), (3.9) N(0) = N0 , kde N(t) je funkce udávaj´ıc´ı poˇcet atom˚ u olova-210 v ˇcase t v gramu látky, r(t) je funkce udávaj´ıc´ı poˇcet atom˚ u olova-210, které vzniknou v ˇcase t v gramu oxidu olovnatého za rok. Protoˇze poloˇcas rozpadu radia-226 je 1600 let a metodu radioaktivn´ıho datován´ı chceme pouˇz´ıt pro rozpoznán´ı stáˇr´ı obraz˚ u, které mˇely v roce 1967 pˇribliˇznˇe bud’ 300 let nebo 20 let, m˚ uˇzeme funkci r(t) povaˇzovat za konstantn´ı. Pak r(t) = r = konst. a u ´ lohu (3.9) m˚ uˇzeme nahradit poˇcáteˇcn´ı u ´ lohou ) N ′ (t) = −λN(t) + r, (3.10) N(0) = N0 . Mnohem v´ıce podrobnost´ı o metodˇe radioaktivn´ıho datován´ı m˚ uˇze ˇctenáˇr nalézt v [2]. Také v pˇr´ıpadˇe u ´ lohy (3.10) jsme schopni, pokud známe poˇcet atom˚ u olova-210 v gramu oxidu olovnatého v dobˇe výroby olovnaté bˇeloby, urˇcit stáˇr´ı obrazu, na kterém je tato barva pouˇzita. K ˇreˇsen´ı této u ´ lohy m˚ uˇzeme 44

opˇet pouˇz´ıt Eulerovu metodu, se kterou jste se seznámili v kapitole 3.1. V naˇs´ı u ´ loze poˇcet atom˚ u olova-210 v gramu oxidu olovnatého v dobˇe výroby barvy bohuˇzel neznáme. I pˇresto jsme schopni rozliˇsit obraz, jehoˇz stáˇr´ı je 300 let, od obrazu, který má 20 let. Je totiˇz známo, jaké bývaj´ı koncentrace radioaktivn´ıho olova-210 v rudách, ze kterých se vyráb´ı oxid olovnatý. Je naprosto nemoˇzné, aby poˇcet atom˚ u olova-210 v gramu rudy, ze které se oxid olovnatý vyrobil pˇresáhl poˇcet 5 · 1011 . Proto m˚ uˇzeme zjistit, pokud známe potˇrebné parametry, zda je moˇzné, aby bylo stáˇr´ı obrazu 300 let. Pˇ r´ıklad 1 Urˇcete, zda je moˇzné,13 aby byl obraz Emauzˇst´ı uˇcedn´ıci opravdu starý 300 let a byl tedy pravý, pokud bylo mˇeˇren´ım zjiˇstˇeno, ˇze v ˇcase mˇeˇren´ı t (tzn. v roce 1967) plat´ı N(t) = 1,42 · 108 , r = 420480. Odhadnˇete chybu aproximace N(t) v ˇcase vytvoˇren´ı obrazu. Pˇ r´ıklad 2 Urˇcete, zda je moˇzné, aby byl obraz Krajkáˇrka opravdu starý 300 let a byl tedy pravý, pokud bylo mˇeˇren´ım zjiˇstˇeno, ˇze v ˇcase mˇeˇren´ı t (tzn. v roce 1967) plat´ı N(t) = 0,25 · 108 , r = 735840. Odhadnˇete chybu aproximace N(t) v ˇcase vytvoˇren´ı obrazu.

Obrázek 3.4: Obrazy Emauzˇst´ı uˇcedn´ıci a Krajkáˇrka

13

V tomto i v následuj´ıc´ım pˇr´ıkladˇe nám m˚ uˇze pomoci, pokud pro ˇreˇsen´ı pouˇzijeme y(t)−y(t−h) ′ tzv. zpˇ etnou diferenci y (t) ≈ , kde h je malé“ a kladné. h ”

45

Poˇc´ ateˇcn´ı u ´lohy

3.4

Populaˇ cn´ı modely

Nyn´ı se pokus´ıme o odhad vývoje poˇctu obyvatel. K tomu se potˇrebujeme nauˇcit pracovat s jednoduchými populaˇcn´ımi modely. Necht’ y(t) oznaˇcuje velikost populace v ˇcase t a y0 popisuje velikost populace v ˇcase t0 . Necht’ p je kladná konstanta, která udává pˇr´ır˚ ustek populace. Pak jednoduchý model vývoje populace nám poskytne následuj´ıc´ı populaˇcn´ı rovnice s poˇcáteˇcn´ı podm´ınkou ) y ′ (t) = py(t), (3.11) y(t0 ) = y0 . ˇ sen´ı u Reˇ ´ lohy (3.11) pomˇernˇe dobˇre aproximuje vývoj populace, která má dostateˇcnˇe velké zásoby potravy a dalˇs´ıch zdroj˚ u a m˚ uˇze neomezenˇe r˚ ust. Lepˇs´ı model dává následuj´ıc´ı populaˇcn´ı rovnice s poˇcáteˇcn´ı podm´ınkou, ve které se nav´ıc objevuje ˇclen qy 2(t) (q je kladná konstanta) y ′ (t) = py(t) − qy 2(t),

y(t0 ) = y0 .

)

(3.12)

ˇ sen´ı u Reˇ ´ lohy (3.12) dobˇre aproximuje vývoj populace, která uˇz je dostateˇcnˇe velká, má omezené zásoby potravy i dalˇs´ıch zdroj˚ u a mezi ˇcleny populace docház´ı k soupeˇren´ı o tyto zdroje (to popisuje ˇclen −qy 2 (t)). Pˇ r´ıklad 1 Pouˇzijte populaˇcn´ı rovnici s poˇcáteˇcn´ı podm´ınkou (3.12) k modelován´ı vývoje poˇctu obyvatel USA v letech 1790–1950. Konstanty p, q byly ˇ t je v roc´ıch. Nav´ıc odhadnuty takto: p = 0,03134, q = 1,5887 · 10−10 . Cas v´ıte, ˇze poˇcet obyvatel v USA v roce 1790 byl 3 929 000. Spoˇc´ıtané hodnoty m˚ uˇzete porovnat se skuteˇcnými hodnotami v následuj´ıc´ı tabulce. Odhadnˇete chybu aproximace poˇctu obyvatel USA v roce 1950.

Rok 1790 1800 1850 1900 1950

Poˇcet obyvatel 3 929 000 5 308 000 23 192 000 75 995 000 150 697 000

Tabulka 3.1: Populace USA v letech 1790–1950

46

3.5

Zmˇ ena teploty tˇ elesa

ˇ Nastal ˇcas na ˇsálek kávy. Cerstvˇ e uvaˇrená káva má teplotu Tk a je ochlazovaná vzduchem v m´ıstnosti o konstantn´ı teplotˇe Tv . Je rozumné modelovat teplotu

kávy jako funkci Tk (t), která se mˇen´ı14 pˇr´ımo u ´ mˇernˇe rozd´ılu teploty vzduchu a kávy s koeficientem k > 0. Necht’ teplota kávy v ˇcase t0 je dána hodnotou Tk0 . Dostáváme tedy následuj´ıc´ı model: ) Tk′ (t) = k(Tv − Tk (t)), (3.13) Tk (t0 ) = Tk0 . Pˇ r´ıklad 1 Pouˇzijte ˇreˇsen´ı u ´ lohy (3.13) k modelován´ı ochlazován´ı kávy, která mˇela v ˇcase t0 = 0 teplotu 100◦ C. Teplota vzduchu má hodnotu Tv = 20◦ C ˇ t je v minutách. Zjistˇete, za jak dlouho se káva a konstanta k = 0,04. Cas ochlad´ı na teplotu niˇzˇs´ı neˇz 50◦ C. Odhadnˇete chybu aproximace ˇreˇsen´ı u ´ lohy (3.13) v ˇcase, kdy se káva ochlad´ı na poˇzadovanou teplotu.

3.6

RC obvod

V této kapitole vytvoˇr´ıme matematický model, který popisuje, jakým zp˚ usobem se mˇen´ı náboj v jednoduchém obvodu s rezistorem a kondenzátorem v závislosti na ˇcase. Uvaˇzujme nyn´ı elektrický obvod, ve kterém je sériovˇe zapojen rezistor a kondenzátor. V obvodu je také v sérii zapojen elektrický zdroj, který do s´ıtˇe dodává elektromotorické napˇet´ı. Schéma obvodu je na obr. 3.5. Nyn´ı vytvoˇrme jednoduchý model, který pop´ıˇse hodnotu elektrického náboje q v obvodu v závislosti na ˇcase t. Odpor rezistoru oznaˇcme R a kapacitu kondenzátoru popiˇsme parametrem C. Elektrický zdroj dodává do obvodu elektromotorické napˇet´ı E(t) (ˇcasto E(t) = A sin (ωt), kde ω, A ∈ R). 14

Tuto zmˇenu matematicky popisuje derivace funkce Tk (t).

47

Poˇc´ ateˇcn´ı u ´lohy

i

R

~ C

Obrázek 3.5: Schéma elektrického obvodu Matematicky lze uvedený systém popsat následuj´ıc´ı Cauchyovou u ´ lohou Rq ′ (t) + C1 q(t) = E(t), q(t0 ) = q0 ,

)

(3.14)

kde q0 oznaˇcuje poˇcáteˇcn´ı hodnotu elektrického náboje. Protoˇze elektrický proud i v ˇcase t je definován jako okamˇzitá zmˇena náboje q a tento vztah m˚ uˇzeme matematicky zapsat jako i(t) = q ′ (t), lze derivován´ım diferenciáln´ı rovnice v (3.14) snadno z´ıskat rovnici pro hodnotu elektrického proudu i procházej´ıc´ıho obvodem v závislosti na ˇcase t: Ri′ (t) + C1 i(t) = E ′ (t). Pˇ r´ıklad 1 Pomoc´ı u ´ lohy (3.14) pop´ıˇseme náboj v elektrickém obvodu, ve kterém je sériovˇe zapojen rezistor s odporem R = 2 Ω, kondenzátor s kapacitou C = 1 F a elektrický zdroj, který do s´ıtˇe dodává elektromotorické napˇet´ı E(t) = 10 V. Hodnota elektrického náboje je na poˇcátku nulová, tj. q(0) = 0 C. Pomoc´ı Eulerovy metody aproximujte chován´ı elektrického náboje v ˇcasovém intervalu h0, 10i s. Odhadnˇete chybu aproximace elektrického náboje v ˇcase t = 10 s. Pˇ r´ıklad 2 Uvaˇzujme obvod bez elektrického zdroje se sériovˇe zapojeným rezistorem a kondenzátorem (se stejnými parametry jako v pˇredchoz´ım pˇr´ıkladˇe). Kondenzátor je na poˇcátku nabit nábojem o hodnotˇe 10 C, tj. q(0) = 10 C. S vyuˇzit´ım modelu (3.14) aproximujte pomoc´ı Eulerovy metody chován´ı elektrického náboje v ˇcasovém intervalu h0, 10i s. Odhadnˇete chybu aproximace elektrického náboje v ˇcase t = 10 s.

48

3.7

Matematick´ e kyvadlo

V této ˇcásti vytvoˇr´ıme matematický model, který popisuje pohyb matematického kyvadla v závislosti na ˇcase. Abychom mohli takovýto model sestrojit, potˇrebujeme pracovat se speciáln´ım typem rovnice – s tzv. diferenciáln´ı rovnic´ı druhého ˇrádu. Nav´ıc si pov´ıme, jak lze diferenciáln´ı rovnici druhého ˇrádu numericky ˇreˇsit. Pro práci s dalˇs´ımi typy diferenciáln´ıch rovnic si zaved’me pojem vektorové funkce. Vektorovou funkc´ ı n re´ aln´ ych promˇ enn´ ych o m sloˇ zk´ ach n m nazýváme kaˇzdé zobrazen´ı z R do R , m, n ∈ N. Vektorová funkce f = (f1 , . . . , fm ) n reálných promˇenných o m sloˇzkách je tedy pˇredpis, který kaˇzdému x = [x1 , x2 , . . . , xn ] ∈ Df ⊂ Rn

pˇriˇrad´ı jednoznaˇcnˇe hodnotu

f (x) = [f1 (x1 , x2 , . . . , xn ), f2 (x1 , x2 , . . . , xn ), . . . , fm (x1 , x2 , . . . , xn )] ∈ Hf ⊂ Rm , kde fi : Rn 7→ R pro i = 1, . . . , m. Skuteˇcnost, ˇze f je vektorovou funkc´ı n reálných promˇenných, zapisujeme pomoc´ı f : Rn 7→ Rm . Soustava diferenci´ aln´ıch rovnic 1. ˇ r´ adu V pˇr´ıpadˇe, ˇze v rovnici (3.1) jsou funkce y = (y1 , . . . , yn ) a f = (f1 , . . . , fn ) vektorové, tzn., ˇze y : (a, b) 7→ Rn a f : Rn+1 7→ Rn , pak m˚ uˇzeme tuto rovnici pˇrepsat do tvaru  y1′ (t) = f1 (t, y1 (t), y2 (t), . . . , yn (t)),     y2′ (t) = f2 (t, y1 (t), y2 (t), . . . , yn (t)),  (3.15) ..  .     yn′ (t) = fn (t, y1 (t), y2 (t), . . . , yn (t)),

který nazýváme soustavou obyˇ cejn´ ych diferenci´ aln´ ıch rovnic 1. ˇ r´ adu. Pokud má nav´ıc systém (3.15) splˇ novat poˇcáteˇcn´ı podm´ınky y(t0 ) = [y10 , . . . , yn0 ] ∈ Rn , tzn. y1 (t0 ) = y10 , y2 (t0 ) = y20 , . . . , yn (t0 ) = yn0, nazveme takovou u ´ lohu opˇet Cauchyovou u ´ lohou. K z´ıskán´ı pˇribliˇzného ˇreˇsen´ı takové u ´ lohy lze pouˇz´ıt rekurzivn´ı vztah (3.4), kde y0 , y1 , . . . , yk , . . . jsou n-rozmˇerné vektory a f : Rn+1 7→ Rn je vektorová funkce. 49

Poˇc´ ateˇcn´ı u ´lohy Obyˇ cejn´ a diferenci´ aln´ı rovnice n-t´ eho ˇ r´ adu Obyˇ cejnou diferenci´ aln´ ı rovnic´ ı n-t´ eho ˇ r´ adu (v explicitn´ ım tvaru) rozum´ıme rovnici tvaru y (n) (t) = f t, y(t), y ′(t), . . . , y (n−1) (t) , (3.16)

kde f : Rn+1 7→ R je zadaná funkce. Pokud k rovnici (3.16) pˇridáme n poˇcáteˇcn´ıch podm´ınek y(t0 ) = y0 , y ′ (t0 ) = y1 , . . . , y (n−1) (t0 ) = yn−1 ,

(yk ∈ R pro k = 0, . . . , n − 1, t0 ∈ R), z´ıskáme Cauchyovu u ´ lohu pro diferenciáln´ı rovnici n-tého ˇrádu. Tuto u ´ lohu m˚ uˇzeme s vyuˇzit´ım substituce z1 (t) = y ′ (t), z2 (t) = y ′′ (t), . . . , zn−1 (t) = y (n−1) (t) pˇrepsat na systém  y ′(t) = z1 (t),     ′  z1 (t) = z2 (t),   .. (3.17) .    ′  zn−2 (t) = zn−1 (t),    ′ zn−1 (t) = f (t, y(t), z1(t), . . . , zn−1 (t)) s poˇcáteˇcn´ımi podm´ınkami

y(t0 ) = y0 , z1 (t0 ) = y1 , . . . , zn−1 (t0 ) = yn−1.

Matematick´ e kyvadlo Nyn´ı sestavme obyˇcejnou diferenciáln´ı rovnici pro popis matematického kyvadla. Zabývejme se zkoumán´ım chován´ı hmotného bodu zavˇeˇseného na tenkém, pevném a nepr˚ utaˇzném vláknu zanedbatelné hmotnosti, zanedbává se odpor vzduchu pˇri pohybu kyvadla i tˇren´ı v závˇesu a gravitaˇcn´ı pole se povaˇzuje za homogenn´ı. T´ıhové zrychlen´ı p˚ usobené t´ımto polem oznaˇc´ıme g. Délka vlákna je ℓ, hmotný bod zavˇeˇsený na vláknu má hmotnost m. Náˇcrtek kyvadla je na obr. 3.6. Vytvoˇrme matematický model kyvadla popisuj´ıc´ı jeho výchylku y od rovnováˇzné polohy v závislosti na ˇcase t. Výchylku y budeme mˇeˇrit v radiánech. S vyuˇzit´ım Newtonova zákona15 lze uvedený systém po urˇcitém zjed15

Pˇresnˇeji druhého Newtonova pohybového zákona, kter´ y je matematicky vyjádˇren vztahem F = ma.

50

y

ℓ

mg sin y

y mg

Obrázek 3.6: Matematické kyvadlo

noduˇsen´ı16 popsat následuj´ıc´ı Cauchyovou u ´ lohou y ′′ (t) = (−g/ℓ) y(t), y(t0 ) = y0 , y ′ (t0 ) = dy0,

)

(3.18)

kde y0 ∈ R oznaˇcuje poˇcáteˇcn´ı výchylku závaˇz´ı (v radiánech) a dy0 ∈ R znaˇc´ı poˇcáteˇcn´ı (´ uhlovou) rychlost závaˇz´ı. Obdobnˇe jako rovnici (3.16) m˚ uˇzeme pomoc´ı substituce pˇrepsat na soustavu (3.17), lze s vyuˇzit´ım substituce z(t) = y ′(t) transformovat u ´ lohu (3.18) na následuj´ıc´ı systém:  y ′ (t) = z(t),   ′ z (t) = (−g/ℓ) y(t), (3.19)   y(t0 ) = y0 , z(t0 ) = dy0 . K pˇribliˇznému ˇreˇsen´ı této soustavy m˚ uˇzeme vyuˇz´ıt dˇr´ıve uvedenou Eulerovu metodu, tj. vektorový tvar formule (3.4).

Pˇ r´ıklad 1 Pomoc´ı u ´ lohy (3.19) pop´ıˇseme kyvadlo, jehoˇz délka je ℓ = 1 m a t´ıhové zrychlen´ı g povaˇzujme rovné 10 m/s−2 . Výchylka kyvadla je na poˇcátku rovna 0,1 (tj. y(0) = 0,1 rad) a poˇcáteˇcn´ı u ´ hlová rychlost kyvadla 16

Toto zjednoduˇsen´ı spoˇc´ıvá v aproximaci y ≈ sin y pro malé“ y. ”

51

Poˇc´ ateˇcn´ı u ´lohy je nulová (tj. z(0) = 0 rad/s). Pomoc´ı Eulerovy metody aproximujte pohyb kyvadla v ˇcasovém intervalu h0, 10i s. Odhadnˇete chybu aproximace výchylky kyvadla v ˇcase t = 10 s. Pˇ r´ıklad 2 Uvaˇzujme opˇet matematické kyvadlo (se stejnými parametry a stejným modelem jako v pˇredchoz´ım cviˇcen´ı). Výchylka kyvadla je na poˇcátku nulová (tj. y(0) = 0 rad) a poˇcáteˇcn´ı u ´ hlová rychlost kyvadla je rovna 1 (tj. z(0) = 1 rad/s). Pomoc´ı Eulerovy metody aproximujte pohyb kyvadla v ˇcasovém intervalu h0, 10i s. Odhadnˇete chybu aproximace výchylky kyvadla v ˇcase t = 10 s.

52

Kapitola 4 Okrajov´ eu ´lohy 4.1

Jednorozmˇ ern´ e (1D) okrajov´ eu ´ lohy a 1D metoda s´ıt´ı

V této kapitole si ukáˇzeme jednu z numerických metod pro pˇribliˇzné ˇreˇsen´ı ´ lohy v jedné dimenzi 1D okrajových u ´ loh, a to tzv. metodu s´ ıt´ ı. Okrajové u sestávaj´ı z pˇr´ısluˇsné diferenciáln´ı rovnice pˇredepsané na daném otevˇreném intervalu a dále z okrajových podm´ınek pˇredepsaných v krajn´ıch bodech tohoto intervalu. V dalˇs´ım výkladu se omez´ıme pouze na tzv. Dirichletovu okrajovou u ´ lohu 2. ˇrádu, kdy nejvyˇsˇs´ı ˇrád derivace vyskytuj´ıc´ı se v diferenciáln´ı rovnici je roven dvˇema a okrajové podm´ınky jsou pˇredepsány ve tvaru hodnoty hledaného ˇreˇsen´ı v krajn´ıch bodech intervalu. 1D Dirichletova u ´ loha At’ k a ℓ jsou libovolné kladné reálné konstanty a f je reálná spojitá funkce. Uvaˇzujme následuj´ıc´ı Dirichletovu1 u ´ lohu: najdˇete funkci u takovou, ˇze ) −k u′′ (x) = f (x) pro x ∈ (0, ℓ), (4.1) u(0) = u(ℓ) = 0. ˇ sen´ı u tohoto problému si lze pˇredstavit jako pr˚ Reˇ uhyb struny délky ℓ, která je uchycena na obou konc´ıch a na n´ıˇz pˇr´ıˇcnˇe p˚ usob´ı s´ıla s hustotou f . Konstanta k souvis´ı s tuhost´ı struny. Takovouto intepretaci u ´ lohy (4.1) si m˚ uˇzeme prohlédnout na obr. 4.1. 1

Johann Peter Gustav Lejeune Dirichlet (1805–1859) – v´ yznamn´ y nˇemeck´ y matematik

53

Okrajové u ´lohy f

u ℓ

Obrázek 4.1: Pr˚ uhyb pˇr´ıˇcnˇe zat´ıˇzené a na obou konc´ıch uchycené struny

Metoda s´ıt´ı pro Dirichletovu u ´ lohu (4.1) Na intervalu (0, ℓ) zvol´ıme pravidelnou s´ıt’ uzl˚ u x0 , x1 , . . . , xn−1 , xn , n ∈ N, viz obr. 4.2. Vzdálenost dvou sousedn´ıch uzl˚ u vyjádˇr´ıme ˇc´ıslem h = ℓ/n. Oznaˇc´ıme si dále fi = f (xi ) pro i = 1, . . . , n − 1. h x0

x1

x2

...

x3

xn−2 xn−1

xn

ℓ

Obrázek 4.2: Pravidelná s´ıt’ uzl˚ u na strunˇe délky ℓ Nyn´ı se budeme snaˇzit aproximovat hodnoty funkce u v daných uzlech, tj. budou nás zaj´ımat ˇc´ısla ui ≈ u(xi ): a) poˇcáteˇcn´ı uzel: u0 = 0, b) vnitˇrn´ı uzly: derivaci druhého ˇrádu vyskytuj´ıc´ı se v diferenciáln´ı rovnici nahrad´ıme diferenˇ cn´ ımi pod´ ıly2 u′ (xi + h2 ) − u′ (xi − h2 ) u (xi ) ≈ ≈ h ′′

=

u(xi +h)−u(xi ) h

− h

u(xi )−u(xi −h) h

=

u(xi+1 ) − 2u(xi ) + u(xi−1 ) ; h2

2

V Eulerovˇe metodˇe jsme na str. 38 pouˇzili k nahrazen´ı derivace tzv. dopˇredné diference. Zde pouˇzijeme tzv. centr´ aln´ ı diferenci: u x + h2 − u x − h2 u′ (x) ≈ , h pˇriˇcemˇz centráln´ı diference pouˇzijeme pro nahrazen´ı druhé i prvn´ı derivace funkce u v dan´ ych uzlech. V´ ysledná aproximace u′′ (xi ) pomoc´ı diferenˇcn´ıch pod´ıl˚ u je pˇr´ımo u ´2 mˇerná h .

54

rovnosti −ku′′ (xi ) = f (xi ) tak nahrad´ıme systémem rovnic k [−ui+1 + 2ui − ui−1 ] = fi , h2

i = 1, 2, . . . , n − 1,

c) koncový uzel: un = 0. Dostáváme tak soustavu n − 1 lineárn´ıch algebraických rovnic o n − 1 neznámých, kterou si m˚ uˇzeme maticovˇe zapsat takto:      2 −1 0 ··· 0 u1 f1  −1     2 −1 0    u2   f2  k   .. . . . .     . . . .. ..   ..  =  ..   . , . . . 2     h   0 −1 2 −1   un−2   fn−2  0 ··· 0 −1 2 un−1 fn−1 | {z } | {z } | {z } =: u

=: K

=: f

zkrácenˇe pak

Ku = f. Matice K se ˇcasto nazývá matic´ ı tuhosti. Vˇsimnˇeme si také, ˇze matice K je tˇr´ıdiagonáln´ı.3 Poznamenejme jeˇstˇe, ˇze jelikoˇz metoda s´ıt´ı vyuˇz´ıvá vˇzdy k nahrazen´ı derivac´ı vhodné diferenˇcn´ı pod´ıly, ˇr´ıká se j´ı také ˇcasto metoda koneˇ cn´ ych diferenc´ ı. Analytick´ eˇ reˇ sen´ı u ´ lohy (4.1) pro konstantn´ı f Pˇredpokládejme, ˇze f (x) := c ∈ R pro vˇsechna x ∈ (0, ℓ) a ˇze u je ˇreˇsen´ım u ´ lohy (4.1). Pak pro vˇsechna x ∈ (0, ℓ) plat´ı: −ku′′ (x) = c, c u′′ (x) = − , k Z c c u′ (x) = − dx = − x + a, a ∈ R, k k Z c c u(x) = − x + a dx = − x2 + ax + b, a, b ∈ R. k 2k

Nyn´ı vezmˇeme v u ´ vahu okrajové podm´ınky a najdˇeme hodnoty konstant a 3

Tˇr´ıdiagonáln´ı matice má nenulové prvky pouze na hlavn´ı diagonále, prvn´ı diagon´ ale pod hlavn´ı diagonálou a prvn´ı diagonále nad hlavn´ı diagonálou.

55

Okrajové u ´lohy a b. u(0) = 0 : u(ℓ) = 0 :

c · 02 + a · 0 + b ⇒ b = 0 2k c c 0 = − · ℓ2 + a · ℓ ⇒ a = ℓ 2k 2k

0=−

Z´ıskali jsme analytické ˇreˇsen´ı Dirichletovy u ´ lohy (4.1) ve tvaru u(x) =

c x(ℓ − x) pro x ∈ h0, ℓi, 2k

viz obr. 4.3. 0

y

u(x) = 21 x(x − 1) −0.1 0

0.2

0.4

0.6

0.8

1

x

Obrázek 4.3: Analytické ˇreˇsen´ı u ´ lohy (4.1) pro volbu c = −1, k = 1 a ℓ = 1 Pˇ r´ıklad 1 Naimplementujte metodu s´ıt´ı pro modelován´ı pr˚ uhybu struny o tuhosti k a délce ℓ, která je uchycená na obou konc´ıch a zat´ıˇzená pˇr´ıˇcnˇe silou s konstantn´ı hustotou c. Své výsledky z´ıskané pomoc´ı metody s´ıt´ı m˚ uˇzete srovnat se známým analytickým ˇreˇsen´ım. Pˇ r´ıklad 2 Svou implementaci metody s´ıt´ı z pˇredchoz´ıho pˇr´ıkladu jednoduˇse modifikujte pro pˇr´ıpad nekonstantn´ıho zat´ıˇzen´ı s hustotou f (x) := (x − 1/3)2 − 1. Vyzkouˇsejte si, jak implementace funguje napˇr´ıklad pro volby k = 2, ℓ = 2 a n = 300. Pˇ r´ıklad 3 Uvaˇzujte strunu délky ℓ = 2, která je uchycená na obou konc´ıch a pod n´ıˇz se ve vzdálenosti d = 0,2 nacház´ı rovinná pˇrekáˇzka. S pˇresnost´ı 10−1 urˇcete nejmenˇs´ı moˇznou tuhost struny4 tak, aby pˇri pˇr´ıˇcném zat´ıˇzen´ı silou o hustotˇe f (x) := (x − 1/3)2 − 1 nedoˇslo ke kontaktu s pˇrekáˇzkou. K ˇreˇsen´ı pouˇzijte metodu s´ıt´ı, pˇriˇcemˇz budete volit n = 300.

4

56

Urˇcete tedy nejmenˇs´ı moˇznou hodnotu koeficientu k.

Kapitola 5 Optimalizaˇ cn´ı u ´lohy 5.1

Nˇ ekolik definic na u ´ vod: gradient, Hessi´ an

Neˇz se zaˇcneme vˇenovat optimalizaˇcn´ım u ´ lohám, seznám´ıme se s nˇekterými pojmy matematické analýzy funkc´ı v´ıce promˇenných. Nejprve si zavedeme gradient a Hessián pro funkci v´ıce promˇenných. Bud’ f : Rn 7→ R,1 x ∈ Rn je vnitˇrn´ı bod Df a k ∈ {1, 2, . . . , n}. Má-li funkce gk : R 7→ R definovaná gk (t) := g(x1 , x2 , . . . , xk−1 , t, xk+1 , . . . , xn ) vlastn´ı derivaci v bodˇe xk , nazýváme ji parci´ aln´ ı derivac´ ı funkce f podle k-t´ e promˇ enn´ e v bodˇ e x a znaˇc´ıme ji ∂f (x) ∂f (x) ′ tj. := gk (xk ) . ∂xk ∂xk Poznamenejme, ˇze podobnˇe m˚ uˇzeme definovat parciáln´ı derivace vyˇsˇs´ıch ˇrád˚ u. Bud’ f : Rn 7→ R funkce, která má v bodˇe x spojité vˇsechny parciáln´ı derivace prvn´ıho ˇrádu. Vektor ∂f (x) ∂f (x) ∂f (x) ∇f (x) := , ,..., ∂x1 ∂x2 ∂xn nazýváme gradient funkce f v bodˇ e x. 1

Takové funkci ˇr´ıkáme re´ aln´ a funkce v´ ıce re´ aln´ ych promˇ enn´ ych.

57

Optimalizaˇcn´ı u ´lohy Bud’ f : Rn 7→ R funkce, která má v bodˇe x spojité vˇsechny parciáln´ı derivace druhého ˇrádu. Matici  ∂ 2 f (x) ∂ 2 f (x) ∂ 2 f (x)  . . . 2 ∂x1 ∂x2 ∂x1 ∂xn ∂x1  ∂ 2 f (x) ∂ 2 f (x) ∂ 2 f (x)    . . . 2 ∂x ∂x ∂x ∂x 2 1 2 ∂xn  2 ∇2 f (x) :=    .. .. .. .. .   . . . ∂ 2 f (x) ∂ 2 f (x) ∂ 2 f (x) ... ∂xn ∂x1 ∂xn ∂x2 ∂x2 n

nazýváme Hessi´ anem nebo Hessovou matic´ ı funkce f v bodˇ e x.

5.2

Rovnov´ aˇ zn´ a poloha tˇ elesa

Nyn´ı se seznám´ıme s tzv. optimalizaˇcn´ımi u ´ lohami, ukáˇzeme si jejich praktickou aplikaci a také jak je lze numericky ˇreˇsit. Mnohem v´ıce mohou ˇctenáˇri nalézt v textu [4]. Optimalizaˇcn´ı u ´ lohy vznikaj´ı ˇcasto pˇri ˇreˇsen´ı praktických u ´ loh. Matematicky jsou tyto u ´ lohy obvykle formulovány jako problémy hledán´ı extrém˚ u n cenov´ e funkce, coˇz je funkce f : R 7→ R. Extrémy cenové funkce ˇcasto hledáme na pˇ r´ ıpustn´ e mnoˇ zinˇ e Ω ⊂ Df . Jelikoˇz bod, v nˇemˇz funkce f nabývá svého minima, je stejný jako bod, v nˇemˇz funkce −f nabývá svého maxima, m˚ uˇzeme za obecnou u ´ lohu optimalizace povaˇzovat problém naj´ıt x ∈ Ω tak, aby platilo f (x) ≤ f (x),

x ∈ Ω.2

(5.1)

Pokud je Ω = Df , mluv´ıme o optimalizaci bez omezen´ ı. Pokud pro Ω plat´ı Ω ⊂ Df a Ω 6= Df , mluv´ıme o optimalizaci s omezen´ ım. Pˇr´ıkladem optimalizaˇcn´ı u ´ lohy bez omezen´ı je napˇr´ıklad problém nalezen´ı rovnováˇzné polohy kuliˇcky o hmotnosti m volnˇe zavˇeˇsené na pruˇzinˇe. Pˇr´ıkladem optimalizaˇcn´ı u ´ lohy s omezen´ım je napˇr´ıklad problém nalezen´ı rovnováˇzné polohy kuliˇcky o hmotnosti m zavˇeˇsené na pruˇzinˇe, která vis´ı nad pˇrekáˇzkou. Uvaˇzujme, ˇze pruˇzina je uchycena v obou pˇr´ıpadech v bodˇe [0, 0]. Tyto problémy m˚ uˇzeme vyjádˇrit jako minimalizaci funkce potenciáln´ı energie tohoto systému. V prvn´ım pˇr´ıpadˇe je Ω = R2 , v druhém pˇr´ıpadˇe pˇr´ıpustná mnoˇzina Ω popisuje oblast, ve které se kuliˇcka m˚ uˇze pohybovat a která leˇz´ı mimo pˇrekáˇzku. Tyto u ´ lohy pak zap´ıˇseme jako u ´ lohu naj´ıt x ∈ Ω tak, aby platila nerovnost (5.1), kde f (x) = x21 + x22 + mx2 a Ω = R2 nebo Ω = {x ∈ R2 : x1 + x2 + 1 ≥ 0, −x1 + x2 + 3 ≥ 0} (pro druhou u ´ lohu viz obr. 5.1). 2

58

ˇ ste u Tuto u ´lohu m˚ uˇzeme zadat také takto: Reˇ ´lohu min f (x).“ ” x∈Ω

x2 −1 −x1 − 1

0

1

m −1

2

3

x1

x1 − 3

−2

Obrázek 5.1: Pˇr´ıklad rovnováˇzného stavu kuliˇcky na pruˇzinˇe (´ uloha s omezen´ım) Odkazy na obecnˇe formulované optimalizaˇcn´ı u ´ lohy najdeme jiˇz ve sta3 rovˇeku. Patˇr´ı mezi nˇe i D´ıdonina u ´ loha. D´ıdó u ´ dajnˇe dostala po pˇristán´ı na africkém pobˇreˇz´ı od m´ıstn´ıho vládce k˚ uˇzi z jednoho býka. Se svým doprovodem mˇela právo usadit se na pozemku, který lze ohraniˇcit touto k˚ uˇz´ı. K˚ uˇzi rozˇrezala na u ´ zké prouˇzky a pomoc´ı nich ohraniˇcili u ´ zem´ı, na kterém zaloˇzila mˇesto Kartágo. Snaha z´ıskat co nejvˇetˇs´ı pozemek vede na problém naj´ıt spojitou kˇrivku dané délky, která ohraniˇcuje co nejvˇetˇs´ı plochu. Dále se budeme zabývat numerickým ˇreˇsen´ım optimalizaˇcn´ıch u ´ loh. Metody ˇ reˇ sen´ı u ´ loh bez omezen´ı Pro u ´ lohu (5.1) existuje ˇrada pˇr´ıstup˚ u, jak ji numericky ˇreˇsit. Pro pˇekné“ ” funkce je velmi efektivn´ı Newtonova metoda. Neˇz se zaˇcneme vˇenovat odvozen´ı Newtonovy metody pro ˇreˇsen´ı optimalizaˇcn´ı u ´ lohy (5.1), odvod´ıme si iteraˇcn´ı metodu pro hledán´ı ˇreˇsen´ı rovnice. Tato metoda se také nazývá Newtonova. Pˇripomeˇ nme, ˇze jsme jiˇz tuto metodu odvodili v kapitole 1.5. Tam jsme ale pro odvozen´ı pouˇzili metodu prostých iterac´ı. Nyn´ı pˇredpokládejme, ˇze v k-tém kroce iteraˇcn´ı metody máme bod aproximuj´ıc´ı ˇreˇsen´ı x∗ rovnice f (x) = 0. Tento bod oznaˇc´ıme xk . Dále pˇredpokládejme, ˇze známe funkˇcn´ı hodnotu funkce f v bodˇe xk a derivaci funkce f v bodˇe xk . Nyn´ı chceme uˇcinit (k + 1)-n´ı krok metody a nalézt bod xk+1 , který lépe“ aproximuje ˇreˇsen´ı ” rovnice f (x) = 0. Vyuˇzit´ım vlastnosti derivace dostaneme z pravo´ uhlého troj´ uheln´ıku v obr. 5.2 následuj´ıc´ı vztah 3

Dle legendy byla D´ıdó fénická princezna, sestra krále Pygmalióna a zakladatelka starovˇekého Kartága.

59

Optimalizaˇcn´ı u ´lohy x2 f (x)

tg α = f ′ (xk ) α ∗

x

xk+1

xk x1

Obrázek 5.2: Newtonova metoda f ′ (xk ) =

f (xk ) . xk − xk+1

(5.2)

Rovnici (5.2) snadno uprav´ıme do tvaru xk+1 = xk −

f (xk ) . f ′ (xk )

(5.3)

T´ım dostáváme pˇredpis, jehoˇz opakovaným pouˇzit´ım dostáváme zpˇresˇ nuj´ıc´ı ∗ se aproximaci ˇreˇsen´ı x rovnice f (x) = 0. Pomoc´ı tohoto pˇredpisu z´ıskáme následuj´ıc´ı algoritmus pro iteraˇcn´ı ˇreˇsen´ı rovnice f (x) = 0. Algoritmus (Newtonova metoda (pro ˇ reˇ sen´ ı rovnice s jednou nezn´ amou)) 1. ε > 0 (ukonˇcuj´ıc´ı podm´ınka) x0 (poˇcáteˇcn´ı bod iteraˇcn´ıho procesu) k=0 x1 = x0 − f (x0 )/f ′ (x0 ) 2. while |xk+1 − xk | ≥ ε (|f (xk+1 )| ≥ ε) k =k+1 xk+1 = xk − f (xk )/f ′ (xk ) end 3. xk+1 aproximuje ˇreˇsen´ı rovnice f (x) = 0 Nyn´ı vyuˇzijeme pˇredpisu (5.3) pro hledán´ı extrému funkce f . Protoˇze pro funkce f , které maj´ı v celém svém definiˇcn´ım oboru derivaci, plat´ı, ˇze extrémy této funkce jsou uvnitˇr definiˇcn´ıho oboru v bodech splˇ nuj´ıc´ıch f ′ (x) = 0,4 4

60

Této podm´ınce ˇr´ıkáme nutná podm´ınka existence lokáln´ıho extrému.

budeme numericky hledat ˇreˇsen´ı rovnice f ′ (x) = 0. Pouˇzit´ım pˇredpisu (5.3) obdrˇz´ıme

xk+1 = xk −

f ′ (xk ) . f ′′ (xk )

(5.4)

Pokud bude funkce f nav´ıc konvexn´ı v celém svém definiˇcn´ım oboru, pak uvedený postup nalezne minimum funkce. Vˇse shrneme do následuj´ıc´ıho algoritmu.

Algoritmus (Newtonova metoda (pro minimalizaci funkce jedn´ e promˇ enn´ e)) 1. ε > 0 (ukonˇcuj´ıc´ı podm´ınka) x0 (poˇcáteˇcn´ı bod iteraˇcn´ıho procesu) k=0 x1 = x0 − f ′ (x0 )/f ′′ (x0 ) 2. while |xk+1 − xk | ≥ ε (|f ′ (xk+1 )| ≥ ε) k =k+1 xk+1 = xk − f ′ (xk )/f ′′ (xk ) end 3. xk+1 aproximuje minimum funkce f (x)

Zobecnˇen´ım pˇredpisu (5.4) pro funkce v´ıce promˇenných z´ıskáme metodu pro minimalizaci funkce v´ıce promˇenných. Zobecnˇen´ım prvn´ı derivace je vektor vˇsech prvn´ıch parciáln´ıch derivac´ı, který nazýváme gradient (oznaˇcujeme jej, jak jsme jiˇz uvedli v u ´ vodu této kapitoly, jako ∇f (x)). Zobecnˇen´ım druhé derivace je matice vˇsech druhých parciáln´ıch derivac´ı, kterou nazýváme Hessova matice nebo Hessián (znaˇc´ımepjej ∇2 f (x)). Absolutn´ı hodnotu |x| zobecn´ıme na velikost vektoru kxk := x21 + x22 + · · · + x2n . Vˇse opˇet pˇrep´ıˇseme do algoritmu. 61

Optimalizaˇcn´ı u ´lohy Algoritmus (Newtonova metoda (pro minimalizaci funkce v´ ıce promˇ enn´ ych)) 1. ε > 0 (ukonˇcuj´ıc´ı podm´ınka) x0 (poˇcáteˇcn´ı bod iteraˇcn´ıho procesu) k=0 −1 x1 = x0 − (∇2 f (x0 )) · ∇f (x0 ) 2. while kxk+1 − xk k ≥ ε (k∇f (xk+1 )k ≥ ε) k =k+1 −1 xk+1 = xk − (∇2 f (xk )) · ∇f (xk ) end 3. xk+1 aproximuje minimum funkce f (x) Poznamenejme, ˇze (∇2 f (xk ))

−1

znaˇc´ı inverzn´ı matici k matici ∇2 f (xk ).

Metody ˇ reˇ sen´ı u ´ loh s omezen´ımi Velmi jednoduchá metoda pˇreveden´ı u ´ lohy s omezen´ım na u ´ lohu bez omezen´ı spoˇc´ıvá v zahrnut´ı nerovnostn´ıch podm´ınek obsaˇzených v popisu pˇr´ıpustné mnoˇziny Ω do cenové funkce tak, ˇze k n´ı pˇriˇcteme vhodný ˇclen, který penalizuje poruˇsen´ı nerovnostn´ıch podm´ınek. Tuto metodu nazýváme metodou kvadratick´ e penalty a pouˇzijeme ji pro u ´ lohu s omezen´ım, která je ve tvaru min f (x), Ω = {x ∈ Rn : g(x) ≤ o} , (5.5) x∈Ω

kde g je vektorová funkce5 a o oznaˇcuje nulový vektor [0, . . . , 0]. Metoda kvadratické penalty aproximuje ˇreˇsen´ı x u ´ lohy (5.5) ˇreˇsen´ım x̺ u ´ lohy minn f̺ (x),

x∈R

1 f̺ (x) = f (x) + ̺kα(x)k2 , 2

α(x) = max{g(x), o}, 6

(5.6)

kde ̺ > 0 je penalizaˇ cn´ ı parametr a kα(x)k2 je penalizaˇ cn´ ı funkce. Intuitivnˇe je zˇrejmé, ˇze kdyˇz je penalizaˇcn´ı parametr ̺ velký, tak ˇreˇsen´ı x̺ , ve kterém je dosaˇzeno minimum penalizované funkce f̺ , nem˚ uˇze být daleko od pˇr´ıpustné mnoˇziny. Dokonce plat´ı, ˇze kdyby ̺ = ∞, pak by bylo minimum f̺ ˇreˇsen´ım p˚ uvodn´ı u ´ lohy. M˚ uˇzeme tedy oˇcekávat, ˇze pro dostateˇcnˇe velké hodnoty penalizaˇcn´ıho parametru ̺ bude ˇreˇsen´ı x u ´ lohy (5.5) bl´ızko x̺ . Je zˇrejmé, ˇze ˇreˇsen´ı penalizované u ´ lohy je typicky bl´ızko pˇr´ıpustné 5

Tento pojem jsme si zavedli jiˇz v u ´vodu kapitoly 3.7. Funkce max znamená maximum po sloˇzkách“, tzn. funkˇcn´ı hodnotou této funkce je ” vektor, kter´ y v i-té sloˇzce obsahuje max{gi (x), 0}. 6

62

mnoˇziny, avˇsak nepatˇr´ı do n´ı. Proto se naˇse penalizaˇcn´ı metoda nazývá také metoda extern´ ı penalty. Pˇ r´ıklad 1 Naleznˇete rovnováˇznou polohu kuliˇcky hmotnosti m volnˇe zavˇeˇsené na pruˇzinˇe (pruˇzina je uchycena v bodˇe [0, 0]), tj. najdˇete ˇreˇsen´ı následuj´ıc´ı minimalizaˇcn´ı u ´ lohy min2 x21 + x22 + mx2 , (5.7) x∈R

kde hmotnost m volte postupnˇe jako m = 1, m = 2, m = 3, m = 4, m = 5 a m = 10. Pˇ r´ıklad 2 Naleznˇete rovnováˇznou polohu kuliˇcky hmotnosti m zavˇeˇsené na pruˇzinˇe (pruˇzina je uchycena v bodˇe [0, 0]), která vis´ı nad pˇrekáˇzkou, tj. najdˇete ˇreˇsen´ı následuj´ıc´ı minimalizaˇcn´ı u ´ lohy min x21 + x22 + mx2 ,

(5.8)

Ω = {x ∈ R2 : x1 + x2 + 1 ≥ 0, −x1 + x2 + 3 ≥ 0},

(5.9)

x∈Ω

kde hmotnost m volte postupnˇe jako m = 1, m = 2, m = 3, m = 4, m = 5 a m = 10 (viz obr. 5.1).

63

64

Kapitola 6 ˇ sen´ı u Reˇ ´loh 6.1

´ Ulohy kapitoly 1

ˇ sen´ı ke kapitole 1.2 Reˇ Pˇ r´ıklad 1 Z obr. 1.5 vid´ıme, ˇze pro obsah S jednotkového kruhu (S = π) plat´ı Z 1√ 1 − x2 dx. S=4 0

Výˇse uvedený integrál nyn´ı pˇribliˇznˇe vypoˇcteme pomoc´ı obdéln´ıkového a lichobˇeˇzn´ıkového pravidla. Nejprve si rozdˇelme interval h0, 1i na n stejných d´ılk˚ u. h = 1 / n; % delka dilku x = 0 : h : 1; % krajni body jednotlivych dilku Nyn´ı si uved’me kód, který pomoc´ı obdéln´ıkového pravidla vypoˇcte aproximaci Sobd obsahu S. c = 0.5 * (x(1:n) + x(2:n+1)); % stredy jednotlivych dilku f = sqrt(1 - c.^2); S obd = 4 * h * sum(f); Kód n´ıˇze vypoˇcte pomoc´ı lichobˇeˇzn´ıkového pravidla aproximaci Slich obsahu S. f = sqrt(1 - x.^2); S lich = 4 * h/2 * ( f(1) + f(n+1) + 2*sum(f(2:n)) ); 65

ˇ sen´ı u Reˇ ´loh

Následuj´ıc´ı tabulka uvád´ı z´ıskané chyby aproximace pro r˚ uzná n: n 102 103 104

|Sobd − π| 3,44 · 10−4 1,09 · 10−5 3,44 · 10−7

|Slich − π| 1,18 · 10−3 3,72 · 10−5 1,18 · 10−6

Pˇ r´ıklad 2 Nejprve se pod´ıvejme na obr. 6.1. Pro obsah S srdce plat´ı

1

S1

S1 0

y

S1 −1

S2

S2

π −2

S2 −3

−3

−2

−1

0

1

2

3

x

Obrázek 6.1: Obsah matematického srdce S = 2S1 + 2S2 = 3π. |{z} |{z} π

2π

Nyn´ı si S vyjádˇreme jako1 Z 2 Z p S = 2 (S1 + S2 ) = 2 1 − (x − 1)2 dx − 0

0

2

π arcsin(x − 1) − dx , 2

pˇriˇcemˇz pro pˇribliˇzné vyˇc´ıslen´ı tˇechto integrál˚ u pouˇzijeme obdéln´ıkové a lichobˇeˇzn´ıkové pravidlo. Nejprve si rozdˇelme interval h0, 2i na n stejných d´ılk˚ u. h = 2 / n; % delka dilku x = 0 : h : 2; % krajni body jednotlivych dilku 1

66

Uvˇedomme si, ˇze plocha S2 leˇz´ı pod osou x.

Nyn´ı vypoˇcteme pomoc´ı obdéln´ıkového pravidla pˇribliˇzný obsah srdce Sobd : c = 0.5 * (x(1:n) + x(2:n+1)); % stredy jednotlivych dilku f1 = sqrt(1 - (c-1).^2); S1 = h * sum(f1); f2 = asin(c-1) - pi/2; S2 = - h * sum(f2); S obd = 2 * (S1 + S2); Pomoc´ı n´ıˇze uvedeného kódu z´ıskáme pˇribliˇzný obsah srdce Slich vypoˇctený lichobˇeˇzn´ıkovým pravidlem: f1 = sqrt(1 - (x-1).^2); S1 = h/2 * ( f1(1) + f1(n+1) + 2*sum(f1(2:n)) ); f2 = asin(x-1) - pi/2; S2 = - h/2 * ( f2(1) + f2(n+1) + 2*sum(f2(2:n)) ); S lich = 2 * (S1 + S2); V následuj´ıc´ı tabulce uvád´ıme z´ıskané chyby aproximace pro r˚ uzná n: n 102 103 104

|Sobd − 3π| |Slich − 3π| 9,73 · 10−4 3,32 · 10−3 3,08 · 10−5 1,05 · 10−5 9,74 · 10−7 3,33 · 10−6

Pˇ r´ıklad 3 K nalezen´ı poˇctu potˇrebných d´ılk˚ u n lze vyuˇz´ıt odhad (1.4). ′ 2 ′′ Jelikoˇz plat´ı (x ) = (2x) = 2, m˚ uˇzeme spoˇc´ıtat n jako ˇreˇsen´ı nerovnice (1 − 0)3 · 2 ≤ 10−4, 24n2 tedy n≥ √

1

= ˙ 28,9. 12 · 10−4 Dˇelen´ı intervalu h0, 1i na 29 (a v´ıce) stejných d´ılk˚ u nám tak zaruˇc´ı, ˇze chyba −4 aproximace bude nejvýˇse 10 . K nalezen´ı poˇctu potˇrebných d´ılk˚ u n m˚ uˇzeme vyuˇz´ıt i Matlab. Nejprve si uvˇedomme, ˇze d´ıky Newtonovu-Leibnizovu vzorci máme 3 1 Z 1 x 1 2 x dx = = . 3 0 3 0 N´ıˇze uvedený kód najde minimáln´ı poˇcet d´ılk˚ u n, pro který chyba aproxi−4 mace err nebude vˇetˇs´ı neˇz 10 . 67


n = 0; err = 1; % inicializace while err > 1e-4 n = n + 1; h x c f S

= 1 / n; % delka dilku = 0 : h : 1; % krajni body jednotlivych dilku = 0.5 * (x(1:n) + x(2:n+1)); % stredy jednotlivych dilku = c.^2; obd = h * sum(f);

err = abs(S obd - 1/3); end Po pr˚ uchodu smyˇckou programu dostáváme n = 29 a err = 9,91 · 10−5 . Je dobré si uvˇedomit, ˇze výˇse uvedený program sám o sobˇe nezaruˇc´ı, ˇze pro n > 29 nebude chyba aproximace vˇetˇs´ı neˇz 10−4 . Tuto informaci nám poskytne, jak uˇz jsme uvedli, odhad (1.4).

ˇ sen´ı ke kapitole 1.3 Reˇ Pˇ r´ıklad 1 Následuj´ıc´ı kód pro dané n vypoˇcte ˇcásteˇcný souˇcet sn ˇrady ∞ X k=1

1 . k(k + 1)

(6.1)

s n = 0; % inicializace for k = 1 : n s n = s n + 1/k/(k + 1); end Nˇekolik spuˇstˇen´ı programu nám napov´ıdá, ˇze by se souˇcet ˇrady (6.1) mohl rovnat jedné - dostáváme napˇr´ıklad: n 102 103 104 105 68

sn 0,990099 0,999001 0,999900 0,999990

Uvedeme si nyn´ı d˚ ukaz, ˇze tomu tak opravdu je. Nejprve si uvˇedomme, ˇze 1 1+k−k 1 1 = = − . k(k + 1) k(k + 1) k k+1 Takˇze

n X

n X 1 1 1 sn = = − = k(k + 1) k k + 1 k=1 k=1 1 1 1 1 1 1 1 1 1 = − + − + − +···+ − =1− . 1 2 2 3 3 4 n n+1 n+1

Jelikoˇz lim [1/(n + 1)] = 0, dostáváme ∞ X k=1

1 = lim sn = 1. k(k + 1)

Pˇ r´ıklad 2 Nejprve si pˇripomeˇ nme, ˇze tˇeleso pohybuj´ıc´ı se rovnomˇernˇe pˇr´ımoˇcaˇre rychlost´ı v uraz´ı za ˇcas t vzdálenost s = vt. Nahlédnˇeme nyn´ı na obr. 6.2. Moucha vyletˇela rychlost´ı vm = 100 km/h z mˇesta A souˇcasnˇe s vlakem a vv t1

vm t1 A′

A vv t1

B′

B

d2

Obrázek 6.2: Vlaky a moucha let´ı vstˇr´ıc druhému vlaku, který vyjel z mˇesta B. Oba vlaky jedou rychlost´ı vv = 10 km/h. Poˇcáteˇcn´ı vzdálenost mezi vlaky (tj. vzdálenost mˇest A a B) je d1 = 90 km. Za ˇcas t1 se moucha setká s vlakem z mˇesta B, pˇriˇcemˇz m´ısto setkán´ı mouchy a vlaku oznaˇc´ıme B ′ . Poˇcáteˇcn´ı vzdálenost d1 mezi vlaky si m˚ uˇzeme rozepsat jako d1 = vm t1 + vv t1 , pˇriˇcemˇz prvn´ı ˇclen souˇctu znaˇc´ı dráhu prvn´ıho letu mouchy (tj. vzdálenost m´ıst A a B’) a druhý ˇclen souˇctu znaˇc´ı dráhu vlaku, kterou ujel z mˇesta B za ˇcas t1 (tj. vzdálenost m´ıst B ′ a B). Odtud t1 =

90 9 d1 = = . vm + vv 100 + 10 11

Vlak z mˇesta A dojel za ˇcas t1 do m´ısta A′ . Vzdálenost m´ıst A a A′ je rovna vv t1 . Nyn´ı se moucha odráˇz´ı od vlaku v m´ıstˇe B ′ a let´ı vstˇr´ıc druhému 69

ˇ sen´ı u Reˇ ´loh vlaku. Vzdálenost vlak˚ u v tomto okamˇziku (tedy vzdálenost m´ıst A′ a B ′ ) lze vyjádˇrit jako d2 = d1 − 2vv t1 = 90 − 2 · 10 ·

9 810 = . 11 11

Pouˇzit´ım stejných u ´ vah jako výˇse si odvod´ıme ˇcas druhého letu mouchy, a to 810 d2 11 t2 = = = vm + vv 100 + 10

9 11

2

.

Obdobnˇe bychom postupovali dále a zjistili bychom, ˇze pro dobu n-tého letu mouchy plat´ı n 9 . tn = 11 Délku n-tého letu mouchy pak vyjádˇr´ıme jako n 9 sn = vm tn = 100 , 11 a proto s15

9 = 100 11

15

= ˙ 4,93,

tj. moucha mezi 15. a 16. odraˇzen´ım uletˇela pˇribliˇznˇe 4,93 km. Celkovou délku letu mouchy z´ıskáme jako souˇcet délek vˇsech2 let˚ u mouchy, tj. jako souˇcet geometrické ˇrady3 ∞ X

n ∞ 9 X 9 sn = 100 = 100 11 9 = 450. 11 1 − 11 n=1 n=1

Celkem tedy moucha nalétala 450 km. Pro zaj´ımavost je na obr. 6.3 znázor∞ P nˇena posloupnost ˇcásteˇcných souˇct˚ u ˇrady sn . n=1

Poznamenejme jeˇstˇe, ˇze pokud by nás nezaj´ımala délka d´ılˇc´ıch let˚ u mouchy, ale jen celková vzdálenost, kterou moucha nalétala, u ´ loha by se dala ˇreˇsit jednoduˇse následuj´ıc´ım zp˚ usobem. Jelikoˇz se vlaky sraz´ı pˇresnˇe v polovinˇe trati (tj. oba ujedou 45 km), dostáváme, ˇze se vlaky sraz´ı za ˇcas ts = 45/vv = 4,5 hodiny. Celkem tedy moucha nalétala vm ts = 100 · 4,5 = 450 km. 2

Let˚ u mouchy je nekoneˇcnˇe mnoho.

3

V následuj´ıc´ım v´ ypoˇctu vyuˇzijeme toho, ˇze

∞ P

n=1

70

qn =

q 1−q

pro |q| < 1.

500

ˇ asteˇcné souˇcty ˇrady C´

450 400 350 300 250 200 150 100 50

0

5

10

15

20

25

30

Lety mouchy

Obrázek 6.3: Celková uraˇzená dráha mouchy v závislosti na letech Pˇ r´ıklad 3 Uvád´ıme program, který vypoˇcte, kolik ˇclen˚ u harmonické ˇrady mus´ıme nejménˇe seˇc´ıst, aby ˇcásteˇcný souˇcet sn harmonické ˇrady mˇel hodnotu alespoˇ n h, kde h je malé“. ” s n = 0; n = 0; % inicializace while s n < h n = n + 1; s n = s n + 1/n; end Po ukonˇcen´ı cyklu je v promˇenné n uloˇzen poˇcet ˇclen˚ u harmonické ˇrady, pro nˇejˇz dostáváme nejmenˇs´ı ˇcásteˇcný souˇcet s hodnotou alespoˇ n h. Posloupnost ˇcásteˇcných souˇct˚ u této ˇrady roste jen velmi pomalu: h 5 10 15 20

n 83 12367 1835421 272400600

Je dokázáno, ˇze napˇr´ıklad pro h = 100 bychom potˇrebovali seˇc´ıst alespoˇ n 15092688622113788323693563264538101449859497 ˇclen˚ u harmonické ˇrady. Tento výsledek ovˇsem nelze z´ıskat pomoc´ı výˇse uvedeného kódu. Pokud bychom se totiˇz tolik ˇclen˚ u harmonické ˇrady pokusili 71

ˇ sen´ı u Reˇ ´loh seˇc´ıst na poˇc´ıtaˇci, museli bychom na výsledek ˇcekat ˇrádovˇe 1017 miliard let (za pˇredpokladu, ˇze naˇsemu poˇc´ıtaˇci trvá pˇriˇcten´ı kaˇzdého nového ˇclenu 10−9 sekund). Uved’me si jeˇstˇe pro zaj´ımavost náznak d˚ ukazu divergence harmonické ˇrady: ∞ X 1 1 1 1 1 1 1 1 1 1 1 = 1+ + + + + + + + + +···+ +··· n 2 3 4 5 6 7 8 9 10 16 n=1

> 1+

1 1 1 1 1 1 1 1 1 1 +··· + + + + + + + + +···+ 2 |4 {z 4} |8 8 {z 8 8} |16 16 {z 16} 1/2

= ∞.

1/2

1/2

Pˇ r´ıklad 4 Následuj´ıc´ı kód pro dané n vypoˇcte ˇcásteˇcný souˇcet sn ˇrady ∞ X (−1)k+1 4 . 2k − 1 k=1

(6.2)

s n = 0; % inicializace for k = 1 : n s n = s n + (-1)^(k+1)/(2*k - 1); end s n = 4 * s n; Nˇekolik spuˇstˇen´ı programu nám napov´ıdá, ˇze by se souˇcet ˇrady (6.2) mohl rovnat π – dostáváme napˇr´ıklad: n 103 104 105 106

sn 3,140593 3,141493 3,141583 3,141592

V´ıce o této aproximaci ˇc´ısla π pomoc´ı ˇc´ıselné ˇrady najdete na konci kapitoly 1.4.

ˇ sen´ı ke kapitole 1.4 Reˇ Pˇ r´ıklad 1 Uvaˇzujme rovinu a na n´ı vyznaˇcenu s´ıt’ rovnobˇeˇzek o vzdálenosti 2ℓ = 1. Na rovinu ház´ıme n-krát jehlu o délce ℓ = 1/2, viz obr. 6.4. Nejdˇr´ıve si uvˇedomme, ˇze se v dalˇs´ı analýze m˚ uˇzeme bez u ´ jmy na obecnosti omezit pouze na dvˇe sousedn´ı rovnobˇeˇzky (volme napˇr´ıklad pˇr´ımky x = 0 a 72

y

x

Obrázek 6.4: Házen´ı jehel na s´ıt’ rovnobˇeˇzek x = 1) a také ˇze y-ová souˇradnice jehly nebude ovlivˇ novat to, zda jehla protne ˇci neprotne nˇekterou z pˇr´ımek. Inspirujme se nyn´ı obr. 6.5 a konstruujme náhodný hod. Nejprve si vygenerujeme náhodné ˇc´ıslo xs z intervalu h0, 1).

α

0

xℓ

xs

xp x

1

Obrázek 6.5: Náhodnˇe hozená jehla Toto ˇc´ıslo bude pˇredstavovat x-ovou souˇradnici stˇredu jehly. Dále si vygenerujeme náhodné ˇc´ıslo α z intervalu h−π/2, π/2), coˇz bude u ´ hel, který sv´ırá jehla s kladným smˇerem osy x. Vypoˇctˇeme nyn´ı hodnoty xℓ a xp pˇredstavuj´ıc´ı x-ové souˇradnice levého a pravého konce jehly: xℓ = xs −

1 1 cos α a xp = xs + cos α. 4 4 73

ˇ sen´ı u Reˇ ´loh Jehla protne pˇr´ımku x = 0 v pˇr´ıpadˇe, ˇze xℓ ≤ 0, nebo pˇr´ımku x = 1 v pˇr´ıpadˇe, ˇze xp ≥ 1. Nyn´ı si uved’me program, který implementuje výˇse zm´ınˇené u ´ vahy. Do promˇenné s se naˇc´ıtá poˇcet hod˚ u, pˇri kterých jehla protla nˇekterou z pˇr´ımek, a z´ıskaná aproximace ˇc´ısla π se uloˇz´ı do promˇenné πaprox . s = 0; % inicializace for i = 1 : n x s = rand; alfa = -pi/2 + rand * pi; x l = x s - 1/4*cos(alfa); x p = x s + 1/4*cos(alfa); if x l <= 0 | x p >= 1 s = s + 1; % jehla protne nekterou z primek end end pi aprox = n / s; Na závˇer poznamenejme, ˇze náˇs program vytvoˇrený k aproximaci π je pouze ilustrativn´ı, jelikoˇz sám π obsahuje. Uved’me si jeˇstˇe tabulku výsledk˚ u pro jedno konkrétn´ı spuˇstˇen´ı programu.4 √ n |πaprox − π| 9,011/ n 102 2 · 10−1 9,01 · 10−1 103 3,3 · 10−2 2,85 · 10−1 4 10 8,89 · 10−4 9,01 · 10−2 105 1,61 · 10−2 2,85 · 10−2 106 7,56 · 10−3 9,01 · 10−3 107 7,19 · 10−4 2,85 · 10−3 108 8,67 · 10−4 9,01 · 10−4 Pˇ r´ıklad 2 Nejprve poˇc´ıtejme pomoc´ı Newtonova-Leibnizova vzorce: Z

0

4

1

x3 x dx = 3 2

1 0

1 = . 3

Uvˇedomme si, ˇze pˇri dalˇs´ım spuˇstˇen´ı obdrˇz´ıme, vzhledem k pˇr´ıtomnosti náhodnˇe generovan´ ych ˇc´ısel, v´ ysledky jiné; jak jsme si vˇsak uvedli √ v kapitole 1.4, chyba aproximace nen´ı s pravdˇepodobnost´ı 95 procent vˇetˇs´ı neˇz 9,011/ n.

74

Nyn´ R 1 2ı si uvedeme kód, který pomoc´ı metody Monte Carlo pˇribliˇznˇe vypoˇcte x dx a spoˇctenou hodnotu uloˇz´ı do promˇenné obsah. Promˇenná n udává 0 poˇcet náhodnˇe generovaných bod˚ u v jednotkovém ˇctverci a do promˇenné s se naˇc´ıtá, kolik z tˇechto bod˚ u leˇz´ı pod grafem ˇci na grafu funkce y = x2 . s = 0; % inicializace for i = 1 : n x = rand; y = rand; % nahodny bod [x,y] z (0,1)x(0,1) if y <= x^2 s = s + 1; end end obsah = s / n; Uved’me si jeˇstˇe tabulku výsledk˚ u pro jedno konkrétn´ı spuˇstˇen´ı programu.5 n 102 103 104 105 106 107 108

√ |obsah − 1/3| 1/ n 3,33 · 10−3 1 · 10−1 6,33 · 10−3 3,16 · 10−2 1,03 · 10−3 1 · 10−2 1,44 · 10−3 3,16 · 10−3 2,04 · 10−4 1 · 10−3 −4 1,23 · 10 3,16 · 10−4 7,32 · 10−6 1 · 10−4

Pˇ r´ıklad 3 Následuj´ıc´ı program R 1 √ implementuje metodu Monte Carlo pro pˇribliˇzný výpoˇcet integrálu 0 1 − x2 dx a vypoˇcte aproximaci πaprox ˇc´ısla π. Jedná se o jednoduchou modifikaci programu z pˇr´ıkladu 2. s = 0; % inicializace for i = 1 : n x = rand; y = rand; % nahodny bod [x,y] z (0,1)x(0,1) if x^2 + y^2 <= 1 s = s + 1; end end 5

Pˇri dalˇs´ım spuˇstˇen´ı obdrˇz´ıme (vzhledem k pˇr´ıtomnosti náhodnˇe generovan´ ych √ ˇc´ısel) v´ ysledky jiné. Chyba aproximace nen´ı s pravdˇepodobnost´ı 75 procent vˇetˇs´ı neˇz 1/ n.

75


pi aprox = 4 * s / n; Uved’me si jeˇstˇe tabulku výsledk˚ u pro jedno konkrétn´ı spuˇstˇen´ı programu.6 n 102 103 104 105 106 107 108

|πaprox − π| 5,84 · 10−2 2,24 · 10−2 5,21 · 10−3 2,83 · 10−3 3,61 · 10−4 9,51 · 10−5 1,27 · 10−4

Pˇ r´ıklad 4 Uvedeme kód, který do promˇenné πaprox uloˇz´ı aproximaci ˇc´ısla π vypoˇctenou pomoc´ı m prvn´ıch ˇclen˚ u Gregoryho ˇrady. s = 0; % inicializace for n = 1 : m s = s + (-1)^(n+1)/(2*n - 1); end pi aprox = 4 * s; V následuj´ıc´ı tabulce jsou výsledky po spuˇstˇen´ı programu: m 102 103 104 105 106

|πaprox − π| 10−2 10−3 10−4 10−5 10−6

M˚ uˇzeme také pouˇz´ıt vztah (1.6) pro výrazné zrychlen´ı konvergence a nahradit tak cyklus ve výˇse uvedeném programu cyklem for n = 1 : m k = 2*n - 1; s = s + (-1)^(n+1)/k * ( 1/2^k + 1/3^k ); end 6

76

Opˇet si uvˇedomme, ˇze pˇri dalˇs´ım spuˇstˇen´ı obdrˇz´ıme v´ ysledky jiné.

Z´ıskané výsledky uvád´ıme v následuj´ıc´ı tabulce: m |πaprox − π| 5 1,49 · 10−4 10 7,4 · 10−8 15 4,87 · 10−11 20 3,69 · 10−14

ˇ sen´ı ke kapitole 1.5 Reˇ Pˇ r´ıklad 1 Máme za u ´ kol nejprve vyˇreˇsit rovnici g(x) := ln(x) + (x − 1)3 = 0 pomoc´ı Newtonovy metody. Po chvilce pˇremýˇslen´ı nás moˇzná napadne, ˇze ˇreˇsen´ım naˇs´ı rovnice je ˇc´ıslo 1. Spoˇctˇeme nyn´ı derivaci funkce g v bodˇe x: g ′(x) =

1 + 3(x − 1)2 . x

Následuj´ıc´ı program pak realizuje ˇreˇsen´ı zadané rovnice Newtonovou metoˇ sen´ı rovnice program nakonec dou pro dané ε a poˇcáteˇcn´ı aproximaci x0 . Reˇ uloˇz´ı do promˇenné x. x = x 0; g = log(x) + (x-1)^3; dg = 1/x + 3*(x-1)^2; x n = x - g/dg; while abs(x n - x) >= epsilon x = x n; g = log(x) + (x-1)^3; dg = 1/x + 3*(x-1)^2; x n = x - g/dg; end x = x n;

77

ˇ sen´ı u Reˇ ´loh Napˇr´ıklad pro ε = 10−4 a x0 = 5 dostáváme ˇreˇsen´ı x = ˙ 0,999999999997277, −12 pro které plat´ı g(x) = ˙ − 2,72 · 10 . Dalˇs´ım u ´ kolem je vyˇreˇsit pomoc´ı Newtonovy metody rovnici 2

g(x) := 1 − x − e−x = 0.

Zamysl´ıme-li se opˇet nad ˇreˇsen´ım naˇs´ı rovnice, mohlo by nás napadnout, ˇze ˇreˇsen´ım je ˇc´ıslo 0. Spoˇctˇeme dále derivaci funkce g v bodˇe x: 2

g ′(x) = −1 + 2xe−x .

Program pro numerickou realizaci pouˇzijeme stejný jako výˇse, aˇz na rozd´ılné definice g(x) a g ′(x): g = 1 - x - exp(-x^2); dg = -1 + 2*x*exp(-x^2); Napˇr´ıklad opˇet pro ε = 10−4 a x0 = 5 dostáváme ˇreˇsen´ı x = ˙ − 4,22 · 10−10 , pro které plat´ı g(x) = ˙ 4,22 · 10−10 . Pˇ r´ıklad 2 Budeme postupovat velmi podobnˇe jako v ˇreˇsen´ı pˇ r´ıkladu 1. Nej√ prve pomoc´ı Newtonovy metody zkusme aproximovat ˇc´ıslo 2. Uvaˇzujme (napˇr´ıklad) rovnici g(x) := x2 − 2 = 0

a poˇc´ıtejme g ′ (x) = 2x. Program pouˇzijeme stejný jako v pˇr´ıkladu 1, aˇz na rozd´ılné definice g(x) a g ′ (x): g = x^2 - 2; dg = 2*x; −4 Pro volby ε = 10√ a x0 = 1 dostáváme ˇreˇsen´ı x = ˙ 1,414213562374690, −12 pro nˇeˇz plat´ı |x − 2| = ˙ 1,59 · 10 . Nyn´ı pomoc´ı Newtonovy metody zkusme aproximovat ˇc´ıslo π. Uvaˇzujme (napˇr´ıklad) rovnici g(x) := sin x = 0

a poˇc´ıtejme g ′(x) = cos x. Program opˇet pouˇzijeme stejný jako v pˇr´ıkladu 1, aˇz na rozd´ılné definice g(x) a g ′ (x): g = sin(x); dg = cos(x); Pro volby ε = 10−3 a x0 = 3 dostáváme ˇreˇsen´ı x = ˙ 3,141592653300477, −10 pro nˇeˇz plat´ı |x − π| = ˙ 2,89 · 10 . 78

ˇ sen´ı ke kapitole 1.6 Reˇ Pˇ r´ıklad 1 Funkci ex nahrad´ıme na okol´ı bodu x0 = 0 Taylorovým polynomem stupnˇe n ∈ N: 7 e0 e0 (x − 0)2 + · · · + (x − 0)n 2! n! 2 n x x = 1+x+ +···+ . 2! n!

ex ≈ Tn (x) = e0 + e0 (x − 0) +

Proto pro x = 1 dostáváme e ≈ Tn (1) = 1 + 1 +

1 1 +···+ . 2! n!

Nyn´ı si uvedeme kód pro aproximaci ˇc´ısla e dle výˇse uvedeného vztahu. % inicializace T n = 1; fact = 1; % promenna pro vypocet faktorialu for k = 1 : n fact = fact * k; T n = T n + 1/fact; end Po dobˇehnut´ı programu máme (pro dané n) v promˇenné Tn uloˇzenu hodnotu Tn (1) aproximuj´ıc´ı ˇc´ıslo e. Následuj´ıc´ı tabulka ilustruje z´ıskané výsledky

n Rn+1 (1) = e − Tn (1) 5 1,62 · 10−3 10 2,73 · 10−8 15 5,06 · 10−14 a na obr. 6.6 m˚ uˇzeme vidˇet, jak se na okol´ı bodu 0 Taylorovy polynomy postupnˇe pˇrimykaj´ı k funkci ex .

7

Pˇripomeˇ nme, ˇze (ex )′ = ex .

79


7

T1 (x) = 1 + x T2 (x) = T1 (x) +

x2 2!

T3 (x) = T2 (x) +

x3 3!

6

5

y

4

3

2

1

ex

0

−1 −2

−1

0

x

1

2

Obrázek 6.6: Nahrazen´ı funkce ex Taylorovými polynomy v bodˇe x0 = 0

Pˇ r´ıklad 2 Naˇs´ım u ´ kolem bude nahradit funkci arctg x na okol´ı bodu x0 = 0 Taylorovým polynomem stupnˇe n ∈ N: arctg x ≈ Tn (x), pˇriˇcemˇz

Tn (x) = arctg (0)+arctg′ (0)(x−0)+

arctg′′ (0) arctg(n) (0) (x−0)2 +· · ·+ (x−0)n . 2! n!

Plat´ı, ˇze arctg (0) = 0, a derivován´ım zjist´ıme, ˇze

80

arctg(n) (x) 1 1 + x2 2x − (1 + x2 )2 1 − 3x2 −2 (1 + x2 )3 x(1 − x2 ) 24 (1 + x2 )4 1 − 10x2 + 5x4 24 (1 + x2 )5 x(3 − 10x2 + 3x4 ) −240 (1 + x2 )6 1 − 21x2 + 35x4 − 7x6 −720 (1 + x2 )7 .. .

n 1 2 3 4 5 6 7 .. .

arctg(n) (0) 1 0 -2 = -(2!) 0 24 = 4! 0 -720 = -(6!) .. .

Odtud lze odvodit, ˇze arctg

(n)

(0) =

0 pro n sudé, n+1 +1 2 (−1) [(n − 1)!] pro n liché.

Pro sudá n tedy plat´ı, ˇze Tn (x) = Tn−1 (x), a pro lichá n dostáváme8 n+1 (n − 1)! xn 2! x3 4! x5 6! x7 + − + · · · + (−1) 2 +1 3! 5! 7! n! 5 7 n 3 n+1 x x x x = x− + − + · · · + (−1) 2 +1 3 5 7 n

Tn (x) = x −

n+1

=

2 X

k=1

(−1)k+1

x2k−1 . 2k − 1

(6.3)

Proto pro x = 1 a lichá n máme n+1 2 X π (−1)k+1 = arctg (1) ≈ Tn (1) = 4 2k − 1 k=1

a tedy n+1

π≈4 8

2 X (−1)k+1

k=1

2k − 1

.

Srovnejte s Leibnizovou ˇradou na str. 16.

81

ˇ sen´ı u Reˇ ´loh Nyn´ı si uvedeme program, který pro dané liché n vypoˇcte aproximaci πaprox ˇc´ısla π dle výˇse uvedeného vztahu. T n = 0; % inicializace for k = 1 : (n + 1)/2 T n = T n + (-1)^(k + 1)/(2*k - 1); end pi aprox = 4 * T n; Následuj´ıc´ı tabulka ilustruje z´ıskané výsledky n |πaprox − π| 10 − 1 2 · 10−3 4 10 − 1 2 · 10−4 105 − 1 2 · 10−5 106 − 1 2 · 10−6 3

a na obr. 6.7 m˚ uˇzeme vidˇet, jak se na okol´ı bodu 0 Taylorovy polynomy Tn (x) postupnˇe pˇrimykaj´ı k funkci arctg x.

1

y

0.5

0

T1 (x) = x

arctg x

T3 (x) = T1 (x) −

−0.5

T5 (x) = T3 (x) + T7 (x) = T5 (x) −

−1 −1.5

−1

−0.5

0

0.5

1

x3 3 x5 5 x7 7

1.5

x

Obrázek 6.7: Nahrazen´ı funkce arctg x Taylorovými polynomy v bodˇe x0 = 0 Pˇ r´ıklad 3 Uvaˇzujme funkci f (x) := x − k(x2 − 2), 82

k > 0.

Plat´ı f ′ (x) = 1 − 2kx a také, ˇze existuje ˇc´ıslo λ takové, ˇze9 ∀x ∈ I :

|f ′ (x)| ≤ λ < 1,

kde I ⊂ (0, 1/k) je libovolný uzavˇrený interval. Aby pevný bod náleˇzel intervalu (0, 1/k), pˇredep´ıˇseme

√

2 funkce f

k ∈ (0, 1/2i. Pouˇzijeme-li pak vˇetu na str. 18, dostáváme,10 ˇze posloupnost (xn ) definovaná pˇredpisem xn+1 := f (xn ) = xn − k(x2n − 2), n = 0, 1, 2, 3, . . . , (6.4) √ konverguje k 2 pro libovolné x0 ∈ (0, 1/k). Nyn´ı si uvedeme kód metody prostých iterac´ı, který pro dané k,√danou ukonˇcuj´ıc´ı podm´ √ınku ε a poˇcáteˇcn´ı aproximaci x0 vydá aproximaci 2aprox pevného bodu 2 funkce f . Promˇenná n slouˇz´ı k naˇc´ıtán´ı poˇctu potˇrebných iterac´ı, tj. k mˇeˇren´ı rychlosti konvergence. % n x x

incializace = 1; = x 0; n = x - k*(x^2 - 2);

while n x x end

abs(x n - x) >= epsilon = n + 1; = x n; n = x - k*(x^2 - 2);

sqrt2 aprox = x n; V následuj´ıc´ı tabulce pak uvád´ıme rychlost konvergence (tj. poˇcet potˇrebných iterac´ı n) pro ε = 10−8 a x0 = 1. 9

Vyˇreˇste si nerovnici |f ′ (x)| < 1. 1 Plat´ı, ˇze f ′ (x) = 1 − 2kx = 0 pro kvadratické funkce f je parabola. 1x =1 2k . Grafem Vrchol této paraboly má souˇradnice 2k , 4k + 2k a jej´ı ramena smˇeˇruj´ı dol˚ u. Dále plat´ı f (0) = f (1/k) = 2k, a proto 1 f ((0, 1/k)) = 2k, + 2k ⊂ (0, 1/k). 4k 1 Vol´ıme-li tedy uzavˇren´ y interval I ⊂ (0, 1/k) tak, aby 2k, 4k + 2k ⊂ I a x0 ∈ I, pak posloupnost (6.4) leˇz´ı v I. 10

83

ˇ sen´ı u Reˇ ´loh √ √ k n | 2aprox − 2| 1/2 21 2,02 · 10−9 1/3 8 2,48 · 10−10 1/4 16 2 · 10−9 Tedy napˇr´ıklad pro k = 1/3 dostáváme rychlejˇs´ı konvergenci neˇz v pˇr´ıpadˇe k = 1/4. Pro konvergenci i pro x0 = 5 vol´ıme11 napˇr´ıklad k = 1/6, aby interval (0, 1/k) = (0, 6) obsahoval x0 . Pro ε = 10−8 dostáváme: √ √ k n | 2aprox − 2| 1/6 28 9,83 · 10−9 1/7 34 9,65 · 10−9 1/8 41 1,32 · 10−8 1/9 47 1,89 · 10−8 1/10 53 2,35 · 10−8 Pˇ r´ıklad 4 Budeme postupovat podobnˇe jako v ˇreˇsen´ı pˇr´ıkladu 3. Pevným bodem funkce f (x) := x + sin x je ˇc´ıslo π. Plat´ı f ′ (x) = 1 + cos x a také, ˇze existuje ˇc´ıslo λ takové, ˇze ∀x ∈ I :

|f ′ (x)| ≤ λ < 1,

kde I ⊂ π2 , 32 π je libovolný uzavˇrený interval. Pouˇzijeme-li pak vˇetu na str. 18, dostáváme,12 ˇze posloupnost (xn ) definovaná pˇredpisem xn+1 := f (xn ) = xn + sin xn , n = 0, 1, 2, 3, . . . , (6.5) konverguje k π pro libovolné x0 ∈ π2 , 32 π . Nyn´ı si uvedeme kód metody prostých iterac´ı, který pro dané ε a danou poˇcáteˇcn´ı aproximaci x0 vydá aproximaci πaprox pevného bodu π funkce f . Promˇenná n slouˇz´ı k naˇc´ıtán´ı poˇctu potˇrebných iterac´ı. % incializace 11

Experimentálnˇe z´ıskáme dobré v´ ysledky i pro volbu k = 1/5 ˇci k = 1/4, konvergence metody vˇsak jiˇz nen´ı zaruˇcena vˇetou na str. 18. 12 Plat´ı π 3 π 3 π 3 π 3 f , π = f ,f π = + 1, π − 1 ⊂ , π . 2 2 2 2 2 2 2 2 Vol´ıme-li tedy uzavˇren´ y interval I ⊂ π2 , 32 π tak, aby π2 + 1, 32 π − 1 ⊂ I a x0 ∈ I, pak posloupnost (6.5) leˇz´ı v I.

84

n = 1; x = x 0; x n = x + sin(x) while n x x end

abs(x n - x) >= epsilon = n + 1; = x n; n = x + sin(x);

pi aprox = x n; V následuj´ıc´ı tabulce pak uvád´ıme výsledky pro ε = 10−1 : x0 2 3 4

6.2

n 3 2 3

|πaprox − π| 1,51 · 10−9 1,76 · 10−11 8,89 · 10−13

´ Ulohy kapitoly 2

ˇ sen´ı ke kapitole 2.2 Reˇ Pˇ r´ıklad 1 Nejprve si oˇc´ıslujme jednotlivé suroviny (ingredience) pouˇzité v naˇsich receptech: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

brambory bujón cibule cukr ˇcesnek ˇcoˇcka droˇzd´ı Hera hoˇrˇcice hovˇez´ı maso jablko jogurt kakao km´ın kokos

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

kuˇrec´ı maso lasagne majoránka máslo mléko mouka mrkev ocet olej oˇr´ıˇsky pepˇr petrˇzelka povidla práˇsek do peˇciva rajˇcata

31 32 33 34 35 36 37 38 39 40 41 42 43 44

rajský protlak skoˇrice slanina soda sójová omáˇcka solamyl s˚ ul sýr tatarka vanilkový cukr vejce vepˇrové maso v´ıno worchesterová omáˇcka

85

ˇ sen´ı u Reˇ ´loh Nyn´ı si uvedeme kód pro vytvoˇren´ı naˇs´ı databáze recept˚ u. Ta je reprezentována incidenˇcn´ı matic´ı inc, jej´ıˇz ˇrádky odpov´ıdaj´ı jednotlivým recept˚ um a sloupce odpov´ıdaj´ı jednotlivým ingredianc´ım. Matice inc obsahuje na ij-té pozici jedniˇcku, pokud i-tý recept obsahuje ingredienci s ˇc´ıslem j; v opaˇcném pˇr´ıpadˇe je na ij-té pozici nula. function inc = recepty() % function inc = recepty() % % inc ... matice incidence receptu a surovin pocet surovin = 44; pocet receptu = 17; inc = zeros(pocet receptu,pocet surovin); inc(1,[4,7,8,20,21,41]) = 1; inc(2,[3,16,24,35,36,37,41,44]) = 1; inc(3,[2,3,5,6,9,18,19,21,22,24,26,27,37,41]) = 1; inc(4,[4,7,20,21,24,29,37,41]) = 1; inc(5,[3,16,19,26,33,35,37,38]) = 1; inc(6,[4,7,14,21,23,24,37]) = 1; inc(7,[1,3,26,37,39]) = 1; inc(8,[4,15,20,21,24,29,40,41]) = 1; inc(9,[5,12,16,21,26,37,38,39,41,43]) = 1; inc(10,[4,13,20,21,24,29,41]) = 1; inc(11,[4,11,13,20,21,24,28,29,32,34,37,41]) = 1; inc(12,[7,21,24,37]) = 1; inc(13,[3,16,24,33,37]) = 1; inc(14,[4,7,19,20,21,25,32,37,41]) = 1; inc(15,[4,8,13,21,25,32,41]) = 1; inc(16,[3,5,10,17,19,20,21,24,26,30,31,37,38,43]) = 1; inc(17,[5,18,24,26,31,36,37,41,42,44]) = 1; Dále si uvedeme program, který posoud´ı bl´ızkost dotazu s recepty v databázi a vyp´ıˇse recept (popˇr. recepty), který je nejbl´ıˇze danému dotazu. V programu je pouˇzita projekce do tˇech sloˇzek (surovin), ve kterých je vektor poˇzadavku (dotazu) nenulový.13 13

Uvaˇzujme napˇr. projekci do n = 4 sloˇzek. Projektovan´ y recept zapiˇsme jako rp = [s1 , s2 , s3 , s4 ] a projektovan´ y poˇzadavek jako pp = [1, 1, 1, 1]. Kos´ınus jejich u ´hlu pak

86

function blizkost(inc) % function blizkost(inc) % % Zjistuje blizkost pozadavku (vektoru danych surovin) % a receptu z databaze % % inc ... matice incidence projektovana na zadane % slozky (suroviny)14 % radky matice normujeme15 for i = 1 : size(inc,1) norm radek = norm(inc(i,:)); if norm radek > 0 inc(i,:) = inc(i,:) / norm radek; end end % velikost projektovaneho pozadavku norm pozadavek = sqrt(size(inc,2)); % urceni blizkosti projektovaneho pozadavku % a projektovanych receptu kosinus uhlu = sum(inc,2) / norm pozadavek;16 shoda = max(kosinus uhlu); indexy = find(kosinus uhlu == shoda); disp([’Dotazu nejlepe vyhovuje recept cislo: ’, int2str(indexy’)]) disp([’Mira shody je: ’, num2str(shoda), ’ (max: 1, min: 0)’]) A koneˇcnˇe dotazy nad naˇs´ı databáz´ı recept˚ u provedeme pomoc´ı následuj´ıc´ıho vypoˇcteme jako cos ϕ =

rp pp s1 + s2 + s3 + s4 √ = . ||rp || ||pp || ||rp || n

14

Hledáme-li napˇr´ıklad recepty obsahuj´ıc´ı brambory a bujón, projektujeme inc na prvn´ı dva jej´ı sloupce. 15 Tj. kaˇzd´ y nenulov´ y ˇrádek matice vynásob´ıme pˇrevrácenou hodnotou jeho velikosti. Velikost ˇrádku matice z´ıskáme pomoc´ı funkce norm. 16 V promˇenné kosinus uhlu je uloˇzen vektor vˇsech kos´ın˚ uu ´hl˚ u projektovan´ ych recept˚ u a projektovaného poˇzadavku.

87

ˇ sen´ı u Reˇ ´loh skriptu: inc = recepty(); % dotaz 1 suroviny = [4,7,8,13,21,32]; blizkost(inc(:,suroviny)) % dotaz 2 suroviny = [37,38,44]; blizkost(inc(:,suroviny)) % dotaz 3 suroviny = [5,9,16,18,26,33,37]; blizkost(inc(:,suroviny)) Uved’me si také výsledky naˇseho prohledáván´ı: • dotazu 1 (cukr, droˇzd´ı, Hera, kakao, mouka a skoˇrice) odpov´ıdá nejlépe recept ˇc. 15 (m´ıra schody je 0,91), • dotazu 2 (s˚ ul, sýr a worchester) odpov´ıdaj´ı nejlépe recepty ˇc. 2, 5, 9, 16 a 17 (m´ıra schody je 0,82), • dotazu 3 (ˇcesnek, hoˇrˇcice, kuˇrec´ı maso, majoránka, pepˇr, slanina a s˚ ul) odpov´ıdá nejlépe recept ˇc. 3 (m´ıra schody je 0,85). Pˇ r´ıklad 2 Postupujeme podobnˇe jako v pˇr´ıkladu 1. Recepty neobsahuj´ıc´ı zadané ingredience vyhledáme pomoc´ı negovaných hodnot matice incidence: inc = recepty(); % dotaz 1 (ABKM) suroviny = [8,12,19,20,38,39]; blizkost(double(∼(inc(:,suroviny)))) % dotaz 2 (vegetarianske) suroviny = [10,16,33,42]; blizkost(double(∼(inc(:,suroviny)))) Výsledky prohledáván´ı jsou: 88

• mléko, Heru, jogurt, máslo, sýr a tatarskou omáˇcku neobsahuj´ı recepty ˇc. 2, 6, 12, 13 a 17, • ˇzádné maso ani slaninu neobsahuj´ı recepty ˇc. 1, 3, 4, 6, 7, 8, 10, 11, 12, 14 a 15.

ˇ sen´ı ke kapitole 2.3 Reˇ Pˇ r´ıklad 1 Uvád´ıme kód ˇreˇsen´ı: a = [1,-1]; b = [2,0]; c = [1,1]; % vrcholy trojuhelniku abc k = 2; % stejnolehlost se stredem v [0,0] a koeficientem k fi = pi / 2; % rotace se stredem v [0,0] a uhlem fi p = [-3,-2]’; % posunuti s vektorem posunuti p P = [a’ b’ c’]; figure; hold on; axis equal; grid on; box on % vykresleni puvodniho trojuhelniku abc fill(P(1,:),P(2,:),’red’); % stejnolehlost Ps = k * P;17 fill(Ps(1,:),Ps(2,:),’yellow’); % rotace R = [cos(fi),-sin(fi); sin(fi),cos(fi)]; % matice rotace Pr = R * Ps; fill(Pr(1,:),Pr(2,:),’blue’); % posunuti Pp = Pr + [p p p]; fill(Pp(1,:),Pp(2,:),’green’); Na obr. 6.8 je znázornˇen grafický výstup programu. k 0 = kI, kde I je jednotkov´ a 0 k matice, pak Ps = T P = kIP = kP . (Plat´ı totiˇz IP = P I = P .) 17

Jelikoˇz pro transformaˇcn´ı matici plat´ı T =

89


4

3

y

2

1

0

−1

−2 −5

−4

−3

−2

−1

0

1

2

3

4

x

Obrázek 6.8: Vzorový troj´ uheln´ık (ˇcervenˇe) a jednotlivé obrazy

Pˇ r´ıklad 2 Uved’me nejprve kód funkce transformace: function at = transformace(a,k,fi,p) % function at = transformace(a,k,fi,p) % % Funkce transformuje zadany bod a. % % a ... a = [x;y] zadany bod % k ... koeficient stejnolehlosti % fi ... uhel otoceni % p ... vektor posunuti % at ... at = [xt;yt] transformovany bod % matice rotace R = [cos(fi),-sin(fi);sin(fi),cos(fi)]; % aplikace transformaci at = R*(k*a) + p; A dále uvedeme kód transformuj´ıc´ı podle zadaných parametr˚ u k, fi a p n-´ uheln´ık zadaný matic´ı P. figure; hold on; axis equal; grid on % zjisteni poctu vrcholu n-uhelniku 90

n = size(P,2); % prealokace vysledne matice Pt = zeros(2,n); % vykresleni vzoru fill(P(1,:),P(2,:),’red’); % aplikace transformaci for i = 1 : n Pt(:,i) = transformace(P(:,i),k,fi,p); end % vykresleni vysledku fill(Pt(1,:),Pt(2,:),’green’); Pˇ r´ıklad 3 Uvád´ıme program, který pro daný poˇcet iterac´ı n provád´ı geometrickou transformaci dle tabulky 2.1. Tento program vyuˇz´ıvá funkci transformace z pˇr´ıkladu 2. x = [0.5; 0.5]; % pocatecni bod figure; hold on; h = plot(x(1),x(2)); set(h,’markersize’,3,’color’,[0 0.3 0],’erasemode’,’none’); axis([-3 3 0 10]) for i = 1 : n r = rand; if r < 0.85 x = transformace(x,0.85,-0.05,[0;1.5]); elseif r < 0.92 x = transformace(x,0.3,1,[0;1.5]); elseif r < 0.99 x = transformace(x,0.25,-1,[0;0.5]); else x = [0 0;0 0.15]*x; end 91


% vykresleni bodu set(h,’xdata’,x(1),’ydata’,x(2)); drawnow end Na obr. 6.9 je znázornˇen výsledný obrazec pro 30000 iterac´ı.18

10 9 8 7

y

6 5 4 3 2 1 0 −3

−2

−1

0

1

2

3

x

Obrázek 6.9: Matematická kapradina

6.3

´ Ulohy kapitoly 3

ˇ sen´ı ke kapitole 3.2 Reˇ Pˇ r´ıklad 1 Následuj´ıc´ı kód nalezne pˇribliˇznˇe pomoc´ı Eulerovy metody ˇcas v hodinách, po kterém je nutno podat dalˇs´ı dávku léku pacientovi trp´ıc´ım astmatem. Za konstantu cl dosad´ıme pro theofylin 48 mℓ/min, za v d dosad´ıme 35 ℓ, za c0 vol´ıme 20 µg/mℓ a za c dolni dosad’me 10 µg/mℓ. Výsledný ˇcas je v promˇenné cas. 18

Jelikoˇz se jedná o transformaci zadanou pomoc´ı pravdˇepodobnosti, dostaneme i pˇri nemˇenném poˇctu iterac´ı pˇri kaˇzdém spuˇstˇen´ı lehce jin´ y v´ ysledek.

92

tN = 10;19 n = 1e5;20 k = cl*60/1000/v d; h = tN / n; t = 0 : h : tN; i = 1; c(1) = c0; while c(i) >= c dolni c(i+1) = c(i) - h*k*c(i); i = i+1; end cas = t(i-1); Pro dostateˇcnˇe velké n (napˇr. pro námi zvolené n = 105 ), zjist´ıme, ˇze dalˇs´ı dávku léku mus´ıme pacientovi podat po 8,423 hodinách (tj. cca po 8 hodinách 25 minutách 23 sekundách). Jinou moˇznost´ı, jak vyˇreˇsit náˇs pˇr´ıklad, je pouˇzit´ı vztahu (3.4) pro nalezen´ı hodnot funkce c(t), vykreslen´ı grafu této funkce a následný odhad ˇcasu, ve kterém je tˇreba podat dalˇs´ı lék (viz obr. 6.10). Pro u ´ plnost si uved’me jeˇstˇe analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı zadaných hodnot dostaneme Cauchyovu u ´ lohu ) c(t), c′ (t) = − 2,88 35 (6.6) c(0) = 20. Derivován´ım se lze pˇresvˇedˇcit, ˇze ˇreˇsen´ım u ´ lohy (6.6) je funkce c¯(t) = 20e−

2,88 t 35

.

Z této funkce sestav´ıme rovnici s neznámou t, která popisuje ˇcas do dalˇs´ı aplikace léku 2,88 10 = 20e− 35 t . Po dalˇs´ı u ´ pravˇe z´ıskáme rovnici ln

1 2,88 =− t. 2 35

19 ˇ

Cas tN jsme uhádli“ tak, aby pˇr´ıliˇs nepˇrevyˇsoval potˇrebn´ y ˇcas pro podán´ı léku. Toto ” hádán´ı“ se dá vyˇreˇsit napˇr. prvotn´ım spuˇstˇen´ım tohoto programu s dostateˇcnˇe velk´ ym ” tN. 20 Poznamenejme, ˇze n je rozumné volit dostateˇcnˇe velké, protoˇze h je nepˇr´ımo u ´mˇerné n. V kapitole 3.1 jsme si uvedli, ˇze chyba aproximace ˇreˇsen´ı pomoc´ı Eulerovy metody z´ avis´ı pˇr´ımo u ´mˇernˇe na kroku h.

93


koncentrace léku v krvi (v µg/mℓ)

20 18 16 14 12 10 8 6 4 2 0

0

1

2

3

4

5

6

7

8

9

10

ˇcas (v hodinách)

Obrázek 6.10: Pokles koncentrace theofylinu v krvi v závislosti na ˇcase Z n´ı uˇz snadno vypoˇc´ıtáme, ˇze dalˇs´ı dávku theofylinu je nutno podat za pˇribliˇznˇe 8,424 hodin, tedy za 8 hodin 25 minut 26 sekund. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace koncentrace theofylinu v ˇcase T = 8,423 hodin z´ıskaná Eulerovou metodou na volbˇe n, a tuto chybu porovnáme s odhadem (3.5). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri. n 104 105 106

h 10−2 10−3 10−4

|caprox (T ) − c¯(T )| 2,9 · 10−3 2,9 · 10−4 2,9 · 10−5

M 2K

eK(T −t0 ) − 1 h 8,2 · 10−3 8,2 · 10−4 8,2 · 10−5

Pˇ r´ıklad 2 Pomoc´ı stejného kódu jako v pˇr´ıkladu 1 najdeme pˇribliˇznˇe ˇcas v hodinách, po kterém je nutno podat dalˇs´ı dávku léku pacientovi s horeˇckami, bolestmi hlavy, sval˚ u ˇci kloub˚ u. Za konstantu cl dosad´ıme pro kyselinu acetylsalicylovou 650 mℓ/min, za v d dosad´ıme 11 ℓ, za c0 vol´ıme 300 µg/mℓ a za c dolni dosad’me 150 µg/mℓ. Výsledný ˇcas je opˇet v promˇenné cas. Jedinou zmˇenu provedeme v nastaven´ı promˇenné tN.

94

tN = 2;21 . . . cas = t(i-1); Pro dostateˇcnˇe velké n (napˇr. pro námi zvolené n = 105 ), zjist´ıme, ˇze dalˇs´ı dávku léku mus´ıme pacientovi podat po 0,195 hodinách (tj. cca po 11 minutách 42 sekundách). Jinou moˇznost´ı, jak vyˇreˇsit náˇs pˇr´ıklad, je pouˇzit´ı vztahu (3.4) pro nalezen´ı hodnot funkce c(t), vykreslen´ı grafu této funkce a následný odhad ˇcasu, ve kterém je tˇreba podat dalˇs´ı lék (viz obr. 6.11).

koncentrace léku v krvi (v µg/mℓ)

300

250

200

150

100

50

0

0

1

2

ˇcas (v hodinách)

Obrázek 6.11: Pokles koncentrace kyseliny acetylsalicylové v krvi v závislosti na ˇcase Nakonec si uved’me opˇet analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı zadaných hodnot dostaneme Cauchyovu u ´ lohu ) c′ (t) = − 39 c(t), 11 (6.7) c(0) = 300. Derivován´ım se lze pˇresvˇedˇcit, ˇze ˇreˇsen´ım u ´ lohy (6.7) je funkce 39

c¯(t) = 300e− 11 t . 21 ˇ

Cas tN jsme opˇet uhádli“ tak, aby pˇr´ıliˇs nepˇrevyˇsoval potˇrebn´ y ˇcas pro podán´ı léku. ”

95

ˇ sen´ı u Reˇ ´loh Z této funkce sestav´ıme rovnici s neznámou t, která popisuje ˇcas do dalˇs´ı aplikace léku 39 150 = 300e− 11 t . Po dalˇs´ı u ´ pravˇe z´ıskáme rovnici ln

1 39 = − t. 2 11

Z n´ı uˇz snadno vypoˇc´ıtáme, ˇze dalˇs´ı dávku kyseliny acetylsalicylové je nutno podat za pˇribliˇznˇe 0,196 hodin, tedy za 11 minut 46 sekund. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace koncentrace kyseliny acetylsalicylové v ˇcase T = 0,195 hodin z´ıskaná Eulerovou metodou na volbˇe n, a tuto chybu porovnáme s odhadem (3.5). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri. M n h |caprox (T ) − c¯(T )| eK(T −t0 ) − 1 h 2K 104 10−2 1,9 5,1 5 10 10−3 1,9 · 10−1 5,3 · 10−1 106 10−4 1,8 · 10−2 5,3 · 10−2

ˇ sen´ı ke kapitole 3.3 Reˇ Pˇ r´ıklad 1 Abychom numericky vyˇreˇsili náˇs pˇr´ıklad, budeme se nyn´ı zabývat u ´ pravou Eulerovy metody, kterou jsme se zabývali v kapitole 3.1. Ale zat´ımco rekurzivn´ı vztah (3.4) nám umoˇzn ˇ uje nalézt hodnoty funkce y pro ˇcasy t > t0 , potˇrebujeme nyn´ı odvodit obdobný vztah pro ˇcasy t < t0 . Hodnota t0 je pro nás ˇcas, ve kterém známe hodnotu hledané funkce popsané diferenciáln´ı rovnic´ı. Napˇr. v naˇsem pˇr´ıpadˇe je t0 rok 1967, ve kterém jsme zmˇeˇrili poˇcet atom˚ u olova-210 v gramu rudy v barvˇe na zkoumaném obrazu. K ˇreˇsen´ı u ´ lohy na intervalu ht−N , t0 i pouˇzijeme opˇet aproximaci derivace funkce y v bodˇe t pomoc´ı tzv. diference y v tomto bodˇe. Ale m´ısto tzv. dopˇredné diference, kterou jsme pouˇzili v kapitole 3.1, nyn´ı pouˇzijeme tzv. zpˇetnou diferenci y(t) − y(t − h) y ′ (t) ≈ , h kde h je malé“ a kladné. ” Po u ´ pravˇe dostaneme y(t − h) ≈ y(t) − hy ′ (t). 96

Pouˇzijeme-li diferenciáln´ı rovnici z (3.2), z´ıskáme vztah y(t − h) ≈ y(t) − hf (t, y(t)).

(6.8)

Zvolme dostateˇcnˇe malou“ pevnou velikost h > 0 a sestrojme posloupnost ” t0 , t−1 := t0 − h, t−2 := t0 − 2h, . . . , t−N := t0 − Nh. Oznaˇcme pomoc´ı yn aproximaci hodnoty pˇresného ˇreˇsen´ı y(tn ). Po pˇridán´ı poˇcáteˇcn´ı podm´ınky dostáváme z (6.8) rekurzivn´ı vztah y0 = y(t0 ), (6.9) yn−1 = yn − hf (tn , yn ), n = 0, −1, . . . , −N. Pomoc´ı následuj´ıc´ıho kódu lze zjistit, zda je moˇzné, aby byl obraz Emauzˇst´ı uˇcedn´ıci opravdu cca 300 let starý. Za konstantu N 0 dosad´ıme 1,42 · 108 a za r dosad´ıme 420480. Pokud nen´ı splnˇena nerovnost N(n+1) <= 5e11, pak obraz nem˚ uˇze být starý 300 let. V tom pˇr´ıpadˇe m˚ uˇzeme ovˇeˇrit, zda je moˇzné, aby vznikl pˇred 20 lety. K tomu staˇc´ı pouˇz´ıt stejný kód s jedinou u ´ pravou, pˇri které za tN dosad´ıme m´ısto hodnoty 300 hodnotu 20. p r = 22;22 lambda = log(2)/(p r); n = 1e5;23 tN = 300; h = tN / n; N(1) = N 0; for i = 1:n N(i+1) = N(i) - h*(-lambda*N(i) + r); end if N(n+1) <= 5e11 disp([’Obraz je pravdepodobne stary 300 let.’]) return end 22

Konstanta p r popisuje poloˇcas rozpadu olova-210. Poznamenejme, ˇze n je rozumné volit dostateˇcnˇe velké, protoˇze h je nepˇr´ımo u ´mˇerné n. V kapitole 3.1 jsme si uvedli, ˇze chyba aproximace ˇreˇsen´ı pomoc´ı Eulerovy metody z´ avis´ı pˇr´ımo u ´mˇernˇe na kroku h. 23

97


koncentrace olova-210 (v atomech na gram rudy)

Pro uvaˇzovaný obraz bychom pomoc´ı uvedeného kódu zjistili, ˇze obraz nem˚ u. 12 ˇze být starý 300 let, protoˇze v tom pˇr´ıpadˇe by N(−300) = 1,6375 · 10 a to nen´ı moˇzné. Obraz je tedy dle uvedených fakt˚ u v kapitole 3.3 starý pˇribliˇznˇe . 20 let. Pro takové stáˇr´ı z´ıskáme N(−20) = 2,5494 · 108 , coˇz je v souladu s bˇeˇznými koncentracemi atom˚ u olova-210. Dalˇs´ı moˇznost´ı, jak vyˇreˇsit náˇs pˇr´ıklad je pouˇzit´ı vztahu (6.9) pro nalezen´ı hodnot funkce N(t), kde t < t0 , vykreslen´ı grafu této funkce a následný odhad, zda je moˇzné, aby byl obraz starý 300 ˇci jen 20 let. Z obr. 6.12 plyne, ˇze obraz nem˚ uˇze být starý 300 let. Zat´ımco z obr. 6.13 je vidˇet, ˇze obraz 20 let m´ıt m˚ uˇze. Stupnice na svislé ose je pro lepˇs´ı pˇrehlednost logaritmická. ˇ Cervená ˇcára na obou obrázc´ıch oznaˇcuje hodnotu 5 · 1011 . Pˇripomeˇ nme, ˇze je naprosto nemoˇzné, aby poˇcet atom˚ u olova-210 v gramu rudy, ze které se oxid olovnatý vyrobil, pˇresáhl tuto hodnotu.

12

10

11

10

10

10

9

10

8

10

7

10

0

50

100

150

200

250

300

ˇcas (v roc´ıch)

Obrázek 6.12: Pokles koncentrace radioaktivn´ıho olova-210 v gramu rudy v závislosti na ˇcase Pro u ´ plnost si uved’me jeˇstˇe analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı zadaných hodnot a vyuˇzit´ım faktu, ˇze poloˇcas rozpadu olova-210 je 22 let, dostaneme Cauchyovu u ´ lohu N ′ (t) = −0,0315N(t) + 420480, N(0) = 1,42 · 108. 98

)

(6.10)


12

10

11

10

10

10

9

10

8

10

7

10

0

5

10

15

20

ˇcas (v roc´ıch)

Obrázek 6.13: Pokles koncentrace radioaktivn´ıho olova-210 v gramu rudy v závislosti na ˇcase Derivován´ım si m˚ uˇzeme ovˇeˇrit, ˇze ˇreˇsen´ım rovnice (6.10) je funkce 420480 −0,0315t 420480 8 ¯ N(t) = 1,42 · 10 − e + . 0,0315 0,0315 Zkusme s vyuˇzit´ım ˇreˇsen´ı u ´ lohy (6.10) urˇcit, zda je moˇzné, aby byl obraz pravý. Pokud je t = −300 rok˚ u, pak . ¯ (−300) = N 1,6383 · 1012 .

Z toho je zˇrejmé, ˇze obraz je podvrh. Pokud zvol´ıme t = −20 rok˚ u, pak . ¯ (−20) = N 2,5494 · 108 .

To odpov´ıdá bˇeˇzným koncentrac´ım atom˚ u olova-210. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace poˇctu atom˚ u olova-210 N(t) v ˇcase vytvoˇren´ı obrazu z´ıskaná Eulerovou metodou na volbˇe n, a tuto chybu porovnáme s odhadem (3.5) (tento odhad plat´ı i pro zpˇetné diference). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri.

99

ˇ sen´ı u Reˇ ´loh n 104 105 106

h 2 · 10−3 2 · 10−4 2 · 10−5

Naprox (−20) − N ¯ (−20) 4796,2 479,6 47,96

M 2K

eK(t0 −(−20)) − 1 h 6682,0 668,2 66,82

Pˇ r´ıklad 2 Pro numerické ˇreˇsen´ı naˇseho pˇr´ıkladu pouˇzijeme obdobné u ´ pravy Eulerovy metody jako v pˇredcházej´ıc´ım pˇr´ıkladu, tj. pouˇzijeme rekurzivn´ı vztah (6.9). Pomoc´ı stejného kódu jako v pˇr´ıkladu 1 m˚ uˇzeme ovˇeˇrit, zda je moˇzné, aby byl obraz Krajkáˇrka opravdu cca 300 let starý. Za konstantu N 0 dosad´ıme 0,25·108 a za r dosad´ıme 735840. Pokud by nebyla splnˇena nerovnost N(n+1) <= 5e11, pak obraz nem˚ uˇze být starý 300 let. V tom pˇr´ıpadˇe m˚ uˇzeme ovˇeˇrit, zda je moˇzné, aby vznikl pˇred 20 lety. K tomu staˇc´ı pouˇz´ıt stejný kód s jedinou u ´ pravou, pˇri které za tN dosad´ıme m´ısto hodnoty 300 hodnotu 20. Pro uvaˇzovaný obraz bychom pomoc´ı uvedeného kódu zjistili, ˇze obraz . m˚ uˇze být starý 300 let, protoˇze N(−300) = 2,096 · 1010 a to je v souladu s bˇeˇznými koncentracemi atom˚ u olova-210. Dalˇs´ı moˇznost´ı, jak vyˇreˇsit náˇs pˇr´ıklad, je pouˇzit´ı vztahu (6.9) pro nalezen´ı hodnot funkce N(t), kde t < t0 , vykreslen´ı grafu této funkce a následný odhad, zda je moˇzné, aby byl obraz starý 300 ˇci jen 20 let. Z obr. 6.14 plyne, ˇze obraz m˚ uˇze být starý 300 let. Stupnice na svislé ose je pro lepˇs´ı pˇrehlednost ˇ logaritmická. Cerven´ a ˇcára na obrázku oznaˇcuje hodnotu 5 · 1011 . Pro u ´ plnost si uved’me jeˇstˇe analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı zadaných hodnot a vyuˇzit´ım faktu, ˇze poloˇcas rozpadu olova-210 je 22 let, dostaneme Cauchyovu u ´ lohu ) N ′ (t) = −0,0315N(t) + 735840, (6.11) N(T ) = 0,25 · 108 . Derivován´ım si m˚ uˇzeme ovˇeˇrit, ˇze ˇreˇsen´ım rovnice (6.11) je funkce 735840 735840 8 ¯ N(t) = 0,25 · 10 − e−0,0315t + . 0,0315 0,0315 Zkusme s vyuˇzit´ım ˇreˇsen´ı u ´ lohy (6.11) urˇcit, zda je moˇzné, aby byl obraz pravý. Pokud je t = −300 rok˚ u, pak . ¯ (−300) = N 2,097 · 1010 .

Z toho je zˇrejmé, ˇze obraz m˚ uˇze být skuteˇcnˇe pravý. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace poˇctu atom˚ u olova-210 N(t) v ˇcase vytvoˇren´ı obrazu z´ıskaná Eulerovou metodou 100


12

10

11

10

10

10

9

10

8

10

7

10

0

50

100

150

200

250

300

ˇcas (v roc´ıch)

Obrázek 6.14: Pokles koncentrace radioaktivn´ıho olova-210 v gramu rudy v závislosti na ˇcase na volbˇe n, a tuto chybu porovnáme s odhadem (3.5) (tento odhad plat´ı i pro zpˇetné diference). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri. M K(t0 −(−300)) Naprox (−300) − N(−300) ¯ n h e − 1 h 2K 4 −2 7 11 10 3 · 10 9,3 · 10 1,3 · 10 105 3 · 10−3 9,4 · 106 1,3 · 1010 106 3 · 10−4 9,4 · 105 1,3 · 109

ˇ sen´ı ke kapitole 3.4 Reˇ Pˇ r´ıklad 1 Následuj´ıc´ı kód odhadne pomoc´ı Eulerovy metody vývoj populace USA v letech 1790–1950. Za konstantu a dosad´ıme poˇcátek zkoumaného intervalu, tedy rok 1790, a za konstantu b dosad´ıme konec zkoumaného intervalu, tedy rok 1950. Samotné hodnoty odhadované velikosti obyvatelstva USA jsou v promˇenné y.

101

ˇ sen´ı u Reˇ ´loh y0 = 3929000; n = 1e5;24 h = (b-a)/n; t = a : h : b; y(1) = y0; for i = 1 : n y(i+1) = y(i) + h*(0.03134*y(i) - 1.5887e-10*y(i)^ 2); end Pomoc´ı tohoto kódu jsme nalezli aproximace poˇct˚ u obyvatel v letech 1800, 1850, 1900 a 1950, které si m˚ uˇzete prohlédnout v následuj´ıc´ı tabulce. Rok 1790 1800 1850 1900 1950

Poˇcet obyvatel 3 929 000 5 308 000 23 192 000 75 995 000 150 697 000

Odhad poˇctu obyvatel Eulerovou metodou 3 929 000 5 335 985 23 191 645 76 868 359 148 675 121

Tabulka 6.1: Populace USA v letech 1790–1950 a odhady jej´ı velikosti Eulerovou metodou Vývoj populace v letech 1790–1950 odhadnutý Eulerovou metodou si m˚ uˇzete prohlédnout na obr. 6.15. Pro u ´ plnost si uved’me jeˇstˇe analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı zadaných hodnot dostaneme Cauchyovu u ´ lohu ) y ′ (t) = 0,03134y(t) − 1,5887 · 10−10 y 2(t), (6.12) y(1790) = 3 929 000. M˚ uˇzeme si derivován´ım ovˇeˇrit, ˇze ˇreˇsen´ım u ´ lohy (6.12) je funkce 123134,86 . 6,2420023 · + 0,03071579977e−0,03134(t−1790) . . Po dosazen´ı dostaneme y¯(1800) = 5 336 024, y¯(1850) = 23 192 502, . . y¯(1900) = 76 871 300 a y¯(1950) = 148 677 189, coˇz jsou velmi dobré odhady skuteˇcných hodnot. y¯(t) =

24

10−4

Poznamenejme, ˇze n je rozumné volit dostateˇcnˇe velké, protoˇze h je nepˇr´ımo u ´mˇerné n. V kapitole 3.1 jsme si uvedli, ˇze chyba aproximace ˇreˇsen´ı pomoc´ı Eulerovy metody závis´ı pˇr´ımo u ´mˇernˇe na kroku h.

102

poˇcet obyvatel USA (v milionech)

150

100

50

0

1800

1850

1900

1950

rok

Obrázek 6.15: Odhad poˇctu obyvatel USA v letech 1790–1950 Nab´ız´ı se pˇrirozenˇe nápad, pouˇz´ıt výˇse uvedený kód pro odhad poˇctu obyvatel USA i po roce 1950, napˇr. pro rok 2014. Pokud bychom to udˇelali, zjistili bychom, ˇze odhad 188 700 000 obyvatel USA, který bychom z´ıskali, je velmi ˇspatný, a je velmi vzdálený skuteˇcnému poˇctu pˇribliˇznˇe 319 000 000 obyvatel USA v tomto roce. To je zp˚ usobeno nepˇresnost´ı modelu, který nebere v u ´ vahu obrovský nár˚ ust legáln´ı (a pravdˇepodobnˇe i ilegáln´ı) imigrace do USA. Vˇzdyt’ zat´ımco na pˇrelomu 40. a 50. let byla kaˇzdoroˇcn´ı legáln´ı imigrace do USA kolem 200 000 lid´ı, tak v 80. letech to jiˇz bylo mezi 500 000 aˇz 600 000 lidmi, od 90. let je jiˇz legáln´ı imigrace kolem 1 000 000 roˇcnˇe a od roku 2005 jiˇz legáln´ı imigrace neklesla nikdy pod 1 000 000 obyvatel. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace poˇctu obyvatel USA v roce 1950 z´ıskaná Eulerovou metodou na volbˇe n, a tuto chybu porovnáme s odhadem (3.5). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri. M K(1950−t0 ) n h |yaprox (1950) − y¯(1950)| e − 1 h 2K 4 −2 4 5 10 1,6 · 10 2,1 · 10 6,1 · 10 105 1,6 · 10−3 2,1 · 103 6,1 · 104 106 1,6 · 10−4 2,1 · 102 6,1 · 103

103


ˇ sen´ı ke kapitole 3.5 Reˇ Pˇ r´ıklad 1 Následuj´ıc´ı kód nalezne pˇribliˇznˇe pomoc´ı Eulerovy metody ˇcas v minutách, který je potˇrebný pro ochlazen´ı kávy ze 100◦ C na teplotu menˇs´ı neˇz 50◦ C. Za T k0 vol´ıme 100 ◦ C, za konstantu k dosad´ıme 0,04 a za T v vol´ıme 20 ◦ C. Výsledný ˇcas je v promˇenné cas. tN = 60;25 n = 1e5;26 h = tN / n; t = 0 : h : tN; i = 1; T k(1) = T k0; while T k(i) >= 50 T k(i+1) = T k(i) + h*(k*(T v - T k(i))); i = i + 1; end cas = t(i); Pro dostateˇcnˇe velké n, zjist´ıme, ˇze se káva ochlad´ı pod 50◦ C za 24,522 minut. Jinou moˇznost´ı, jak vyˇreˇsit náˇs pˇr´ıklad, je pouˇzit´ı vztahu (3.4) pro nalezen´ı hodnot funkce Tk (t), vykreslen´ı grafu této funkce a následný odhad ˇcasu, který je potˇreba, neˇz se káva ochlad´ı na poˇzadovanou teplotu (viz obr. 6.16). Uved’me si jeˇstˇe analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı uvedených hodnot dostaneme Cauchyovu u ´ lohu ) Tk′ (t) = 0,04(20 − Tk (t)), (6.13) Tk (0) = 100. Derivován´ım lze snadno ovˇeˇrit, ˇze ˇreˇsen´ım rovnice (6.13) je funkce Tk (t) = 80e−0,04t + 20. Z toho plyne 50 = 80e−0,04t + 20. 25 ˇ

Cas tN jsme uhádli“ tak, aby pˇr´ıliˇs nepˇrevyˇsoval potˇrebn´ y ˇcas pro ochlazen´ı kávy. ” Toto hádán´ı“ se dá vyˇreˇsit napˇr. prvotn´ım spuˇstˇen´ım tohoto programu s dostateˇcnˇe ” velk´ ym tN. 26 Poznamenejme, ˇze n je rozumné volit dostateˇcnˇe velké, protoˇze h je nepˇr´ımo u ´mˇerné n. V kapitole 3.1 jsme si uvedli, ˇze chyba aproximace ˇreˇsen´ı pomoc´ı Eulerovy metody závis´ı pˇr´ımo u ´mˇernˇe na kroku h.

104

teplota kávy (ve stupn´ıch Celsia)

100 90 80 70 60 50 40 30 20 10 0

0

5

10

15

20

25

30

35

40

ˇcas (v minutách)

Obrázek 6.16: Pokles teploty kávy v závislosti na ˇcase To znamená, ˇze káva se ochlad´ı pod 50◦ C za pˇribliˇznˇe 24,521 minut. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace teploty kávy v ˇcase T = 24,522 minut z´ıskaná Eulerovou metodou na volbˇe n, a tuto chybu porovnáme s odhadem (3.5). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri. M Tk aprox (T ) − Tk (T ) n h eK(T −t0 ) − 1 h 2K 104 6 · 10−2 3,5 · 10−2 1,6 · 10−1 105 6 · 10−3 3,5 · 10−3 1,6 · 10−2 6 −4 −4 10 6 · 10 3,5 · 10 1,6 · 10−3

ˇ sen´ı ke kapitole 3.6 Reˇ Pˇ r´ıklad 1 Následuj´ıc´ı kód slouˇz´ı k z´ıskán´ı aproximace pr˚ ubˇehu elektrického náboje v elektrickém obvodu na obr. 3.5, který sledujeme v ˇcasovém intervalu h0, 10i s. Za konstantu a dosad´ıme poˇcátek zkoumaného intervalu, tedy 0, a za konstantu b dosad´ıme konec zkoumaného intervalu, tedy 10. Samotné hodnoty odhadovaného elektrického náboje jsou v promˇenné q. 105


q0 = 0; E = 10; R = 2; C = 1; n = 1e5;27 h = (b-a) / n; t = a : h : b; q(1) = q0; for i = 1 : n q(i+1) = q(i) + h * (1/R *(E - 1/C * q(i))); end Pr˚ ubˇeh elektrického náboje odhadnutý Eulerovou metodou si m˚ uˇzete prohlédnout na obr. 6.17.

11

elektrický náboj (v Coulombech)

10 9 8 7 6 5 4 3 2 1 0

0

2

4

6

8

10

ˇcas (v sekundách)

Obrázek 6.17: Odhad pr˚ ubˇehu elektrického náboje v elektrickém obvodu Uved’me si jeˇstˇe analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı zadaných hod27

Poznamenejme, ˇze n je rozumné volit dostateˇcnˇe velké, protoˇze h je nepˇr´ımo u ´mˇerné n. V kapitole 3.1 jsme si uvedli, ˇze chyba aproximace ˇreˇsen´ı pomoc´ı Eulerovy metody závis´ı pˇr´ımo u ´mˇernˇe na kroku h.

106

not dostaneme Cauchyovu u ´ lohu 2q ′ (t) + q(t) = 10, q(0) = 0.

)

(6.14)

Pomoc´ı derivován´ı m˚ uˇzeme ovˇeˇrit, ˇze ˇreˇsen´ım u ´ lohy (6.14) je funkce 1

q¯(t) = −10e− 2 t + 10. Graf této funkce odpov´ıdá pr˚ ubˇehu na obr. 6.17. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace elektrického náboje v ˇcase t = 10 sekund z´ıskaná Eulerovou metodou na volbˇe n, a tuto chybu porovnáme s odhadem (3.5). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri. M h |qaprox (10) − q¯(10)| eK(10−t0 ) − 1 h n 2K 104 10−3 8,4 · 10−5 4,0 · 10−3 5 −4 −6 10 10 8,4 · 10 4,0 · 10−4 106 10−5 8,4 · 10−7 4,0 · 10−5 Pˇ r´ıklad 2 Pomoc´ı stejného kódu jako v pˇr´ıkladu 1 z´ıskáme aproximaci pr˚ ubˇehu elektrického náboje v elektrickém obvodu na obr. 3.5, který sledujeme v ˇcasovém intervalu h0, 10i s. Za konstanty a a b opˇet dosad´ıme poˇcátek a konec zkoumaného intervalu, tedy 0 a 10. Samotné hodnoty odhadovaného elektrického náboje jsou zase v promˇenné q. Jedinou zmˇenu provedeme v nastaven´ı konstanty E a poˇca´teˇcn´ı podm´ınky q0. q0 = 10; E = 0; . . .

Pr˚ ubˇeh elektrického náboje odhadnutý Eulerovou metodou si m˚ uˇzete prohlédnout na obr. 6.18. Uved’me si opˇet jeˇstˇe analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı zadaných hodnot dostaneme Cauchyovu u ´ lohu ) 2q ′ (t) + q(t) = 0, (6.15) q(0) = 10. 107


11

elektrický náboj (v Coulombech)

10 9 8 7 6 5 4 3 2 1 0

0

2

4

6

8

10


Obrázek 6.18: Odhad pr˚ ubˇehu elektrického náboje v elektrickém obvodu Pomoc´ı derivován´ı m˚ uˇzeme ovˇeˇrit, ˇze ˇreˇsen´ım u ´ lohy (6.15) je funkce 1

q¯(t) = 10e− 2 t . Graf této funkce odpov´ıdá pr˚ ubˇehu na obr. 6.18. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace elektrického náboje v ˇcase t = 10 sekund z´ıskaná Eulerovou metodou na volbˇe n, a tuto chybu porovnáme s odhadem (3.5). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri. M n h |qaprox (10) − q¯(10)| eK(10−t0 ) − 1 h 2K 104 10−3 8,4 · 10−5 3,7 · 10−3 5 −4 −6 10 10 8,4 · 10 3,7 · 10−4 106 10−5 8,4 · 10−7 3,7 · 10−5

ˇ sen´ı ke kapitole 3.7 Reˇ Pˇ r´ıklad 1 Následuj´ıc´ı kód slouˇz´ı k z´ıskán´ı aproximace pr˚ ubˇehu výchylek kyvadla v ˇcasovém intervalu h0, 10i s. Za konstantu a dosad´ıme poˇcátek zkoumaného intervalu, tedy 0, a za konstantu b dosad´ıme konec zkoumaného 108

intervalu, tedy 10. Samotné hodnoty odhadované výchylky jsou v promˇenné y. y0 = 0.1; z0 = 0; g = 10; l = 1; n = 1e5;28 h = (b-a) / n; t = a : h : b; y(1) = y0; z(1) = z0; for i = 1 : n y(i+1) = y(i) + h * z(i); z(i+1) = z(i) + h * ((-g/l)*y(i)); end Pr˚ ubˇeh výchylek odhadnutý Eulerovou metodou si m˚ uˇzete prohlédnout na obr. 6.19. Uved’me si jeˇstˇe analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı zadaných hodnot dostaneme Cauchyovu u ´ lohu ) y ′′(t) = −10y(t), (6.16) y(0) = 0,1, y ′(0) = 0. Pomoci derivován´ı m˚ uˇzeme ovˇeˇrit, ˇze ˇreˇsen´ım u ´ lohy (3.18) je funkce y¯(t) = 0,1 cos

√

10t .

Graf této funkce odpov´ıdá pr˚ ubˇehu na obr. 6.19. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace výchylky kyvadla v ˇcase t = 10 sekund z´ıskaná Eulerovou metodou na volbˇe n, a tuto chybu porovnáme s odhadem (3.5). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri. 28

Poznamenejme, ˇze n je rozumné volit dostateˇcnˇe velké, protoˇze h je nepˇr´ımo u ´mˇerné n. V kapitole 3.1 jsme si uvedli, ˇze chyba aproximace ˇreˇsen´ı pomoc´ı Eulerovy metody z´ avis´ı pˇr´ımo u ´mˇernˇe na kroku h.

109


výchylka kyvadla (v radiánech)

0.1

0.05

0

−0.05

−0.1 0

2

4

6

8

10


Obrázek 6.19: Odhad pr˚ ubˇehu výchylek kyvadla n 104 105 106

h 10−3 10−4 10−5

|yaprox (10) − y¯(10)| 5,0 · 10−3 4,5 · 10−4 4,9 · 10−5

M 2K

eK(10−t0 ) − 1 h 3,7 3,5 · 10−1 3,5 · 10−2

Pˇ r´ıklad 2 Pomoc´ı stejného kódu jako v pˇr´ıkladu 1 z´ıskáme aproximaci pr˚ ubˇehu výchylek kyvadla v ˇcasovém intervalu h0, 10i s. Za konstanty a a b opˇet dosad´ıme poˇcátek a konec zkoumaného intervalu, tedy 0 a 10. Samotné hodnoty odhadované výchylky jsou zase v promˇenné y. Jedinou zmˇenu provedeme v nastaven´ı poˇcáteˇcn´ıch podm´ınek y0 a z0. y0 = 0; z0 = 1; . . . Pr˚ ubˇeh výchylek odhadnutý Eulerovou metodou si m˚ uˇzete prohlédnout na obr. 6.20. Uved’me si opˇet jeˇstˇe analytické ˇreˇsen´ı naˇs´ı u ´ lohy. Po dosazen´ı zadaných hodnot dostaneme Cauchyovu u ´ lohu ) y ′′(t) = −10y(t), (6.17) y(0) = 0, y ′(0) = 1. 110

0.4

výchylka kyvadla (v radiánech)

0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4

0

2

4

6

8

10


Obrázek 6.20: Odhad pr˚ ubˇehu výchylek kyvadla

Pomoci derivován´ı m˚ uˇzeme ovˇeˇrit, ˇze ˇreˇsen´ım u ´ lohy (3.19) je funkce

y¯(t) =

√

√ 10 sin 10t . 10

Graf této funkce odpov´ıdá pr˚ ubˇehu na obr. 6.20. V následuj´ıc´ı tabulce si jeˇstˇe ukáˇzeme, jak závis´ı chyba aproximace výchylky kyvadla v ˇcase t = 10 sekund z´ıskaná Eulerovou metodou na volbˇe n, a tuto chybu porovnáme s odhadem (3.5). Vyˇc´ıslen´ı odhadu (3.5) je nároˇcnˇejˇs´ı, a proto v tabulce uvád´ıme jen hodnotu tohoto výrazu a podrobný výpoˇcet pˇrenecháme zv´ıdavému ˇctenáˇri.

n 104 105 106

h 10−3 10−4 10−5

|yaprox (10) − y¯(10)| 3,3 · 10−3 3,3 · 10−4 3,2 · 10−5

M 2K

eK(10−t0 ) − 1 h 11,6 1,11 1,10 · 10−1

111


6.4

´ Ulohy kapitoly 4

ˇ sen´ı ke kapitole 4.1 Reˇ Pˇ r´ıklad 1 Uvedeme si nejprve kód ˇreˇsen´ı: function u = struna1(c,k,l,n) % u = struna1(c,k,l,n) % % Funkce vyresi metodou siti tuto Dirichletovu ulohu % -k*u’’(x) = c na (0,l) % u(0) = u(l) = 0 % % c ... konstantni hustota zatezujici sily % k ... tuhost struny % l ... delka struny % n+1 ... pocet uzlu na strune (pravidelna sit) % u ... vysledny pruhyb struny % vzdalenost mezi sousednimi uzly h = l / n; % matice tuhosti K K = diag(2*ones(n-1,1)) + diag(-ones(n-2,1),1) + diag(-ones(n-2,1),-1); K = k/h^2 * K; % vektor prave strany f f = c * ones(n-1,1); % reseni soustavy - pruhyb struny v uzlech u = K \ f; u = [0; u; 0]; % vykresleni reseni figure x = 0 : h : l; plot(x,u) % srovnani s analytickym resenim u a = c/2/k * x .* (l-x); err abs = max(abs(u - u a’)) % maximalni absolutni chyba 112

Na obr. 6.21 vid´ıme pr˚ uhyb struny pro c = −1, k = 2, ℓ = 2 a r˚ uzná n. M˚ uˇzeme si vˇsimnout, ˇze pr˚ uhyby struny se v uzlových bodech vˇzdy pˇresnˇe 0

ua (x) = x4 (x − 2)

−0.05

n=2 n=4

−0.1

y

n=8

−0.15

−0.2

−0.25 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x

Obrázek 6.21: Analytické ˇreˇsen´ı ua pro c = −1, k = 2 a ℓ = 2 a vypoˇctené pr˚ uhyby struny pro r˚ uzná n shoduj´ı s analytickým ˇreˇsen´ım.29 Je to proto, ˇze metoda s´ıt´ı je pˇresná pro f := c = konst. – nahrazen´ı pomoc´ı diferenc´ı na str. 54 je totiˇz v tomto pˇr´ıpadˇe pˇresné. Pˇ r´ıklad 2 Nejprve si vyjádˇr´ıme analytické ˇreˇsen´ı u ´ lohy pro zatˇeˇzuj´ıc´ı s´ılu 2 30,31 s hustotou f (x) := (x − 1/3) − 1: 4 x − 13 x2 ua (x) = − + + ax + b, 12k 2k " # 4 ℓ − 13 1 1 1 a= − ℓ2 − , b= . 2kℓ 6 486 972k 29

Po dobˇehnut´ı programu s jak´ ymkoli vstupn´ım n dostáváme v promˇenné err abs vˇzdy (poˇc´ıtaˇcovou) nulu. 30 To dostaneme stejn´ ym zp˚ usobem, jak´ y jsme pouˇzili na str. 55 pro konstatn´ı f . 31 Vˇsimnˇeme si, ˇze f (x) < 0 pro x ∈ h0, 4/3 ) a f (x) > 0 pro x ∈ (4/3, 2i, tzn. ˇze strunu zatlaˇcujeme smˇerem dol˚ u v intervalu h0, 4/3 ) a zvedáme nahoru v intervalu (4/3, 2i.

113

ˇ sen´ı u Reˇ ´loh Kód pouˇzijeme stejný jako v pˇr´ıkladˇe 1, jen uprav´ıme pˇredpis pro vektor pravé strany: % x f f

vektor prave strany f = h : h : l-h; = (x - 1/3).^2 - 1; = f’;

a uprav´ıme i ˇcást kódu pro srovnán´ı s analytickým ˇreˇsen´ım: 32 % srovnani s analytickym resenim a = ( (l-1/3)^4/6 - l^2 - 1/486 ) / 2 / k / l; b = 1 / 972 / k; u a = -(x-1/3).^4/12/k + x.^2/2/k + a*x + b; err abs = max(abs(u - u a’)) % maximalni absolutni chyba Na obr. 6.22 vid´ıme pr˚ uhyb struny pro k = 2 a ℓ = 2 a r˚ uzná n. V následuj´ıc´ı 0.02

ua (x) =

0

n=2 n=4 n=8

−0.02

y

−0.04

x (2 72

− x)(3x2 + 2x − 12)

−0.06 −0.08 −0.1 −0.12 −0.14 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x

Obrázek 6.22: Analytické ˇreˇsen´ı ua pro k = 2 a ℓ = 2 a vypoˇctené pr˚ uhyby struny pro r˚ uzná n tabulce pak uvád´ıme výsledky pro k = 2 a ℓ = 2: 32

Takto upravenou funkci struna1 m˚ uˇzeme pˇrejmenovat na struna2, pˇriˇcemˇz vstupn´ımi argumenty funkce struna2 budou pouze promˇenné k, l a n.

114

n 2 4 8 100 300 500 1000

err abs 4,17 · 10−2 1,04 · 10−2 2,6 · 10−3 1,67 · 10−5 1,85 · 10−6 6,67 · 10−7 1,67 · 10−7

Pˇ r´ıklad 3 K ˇreˇsen´ı pouˇzijeme následuj´ıc´ı kód: l = 2; d = 0.2; n = 300; k = 2; % vhodna inicializace while 1 k = k - 0.1; u = struna2(k,l,n); if min(u) <= -d break end end k = k + 0.1; Funkce struna2, která se volá uvnitˇr cyklu, je kód ˇreˇsen´ı pˇr´ıkladu 2 (tedy upravený kód ˇreˇsen´ı pˇr´ıkladu 1). Inicializace tuhosti k se odv´ıj´ı od ˇreˇsen´ı pˇr´ıkladu 2, kde jsme mimo jiné zjistili, ˇze pro k = 2 se struna nacház´ı celá nad pˇrekáˇzkou. Po dobˇehnut´ı programu zjist´ıme, ˇze struna mus´ı m´ıt (s pˇresnost´ı 10−1 ) tuhost nejménˇe k = 1,2, aby se nedotkla zadané pˇrekáˇzky. Na obr. 6.23 jsou znázornˇeny pr˚ uhyby stuny pro postupnˇe klesaj´ıc´ı tuhost.

6.5

´ Ulohy kapitoly 5

ˇ sen´ı ke kapitole 5.2 Reˇ Pˇ r´ıklad 1 Rovnováˇznou polohu volnˇe zavˇeˇsené kuliˇcky dané hmotnosti najdeme jako ˇreˇsen´ı minimalizaˇcn´ı u ´ lohy (5.7). K tomu pouˇzijeme Newtonovu 115


0.05

0

k = 1,9 k = 1,8 .. .

y

−0.05

−0.1

−0.15

k = 1,2 k = 1,1

−0.2

−0.25

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x

Obrázek 6.23: Jednotlivé pr˚ uhyby struny pro klesaj´ıc´ı tuhost metodu pro minimalizaci funkce v´ıce promˇenných. Hledanou rovnováˇznou polohu najdeme pomoc´ı následuj´ıc´ıho kódu. Hmotnost kuliˇcky je v promˇenné m, za kterou postupnˇe dosazujeme hodnoty 1, 2, 3, 4, 5 a 10. function [y,it] = newton(df,ddf,x0,epsilon) %[y,it] = newton(df,ddf,x0,epsilon) % % Zakladni Newtonova metoda pro minimalizaci fce f(x) % % y ... minimum % it ... pocet iteraci % % df ... gradient minimalizovane fce f(x) % ddf ... Hessian minimalizovane fce f(x) % x0 ... pocatecni iterace % epsilon ... ukoncujici podminka global m m = 1;

116

it = 0; x = x0; grad = feval(df,x); Hess = feval(ddf,x); while norm(grad) >= epsilon33 it = it + 1; x = x - Hess \ grad; grad = feval(df,x); Hess = feval(ddf,x); end y = x; Newtonova metoda potˇrebuje v kaˇzdém kroce iteraˇcn´ıho procesu vyˇc´ıslit gradient a Hessián. Gradient funkce f (x1 , x2 ) = x21 +x22 +mx2 vypoˇcteme pomoc´ı následuj´ıc´ıho kódu. function y = df(x) global m y = [2*x(1); 2*x(2) + m]; Hessián této funkce vypoˇcteme pomoc´ı kódu: function y = ddf(x) global m y = [2 0; 0 2]; Pokud ukonˇcuj´ıc´ı podm´ınku epsilon zvol´ıme jako 10−4 , z´ıskáme pro poˇzadované hmotnosti rovnováˇzné polohy kuliˇcky uvedené v tabulce 6.2. Pˇ r´ıklad 2 Pomoc´ı stejného kódu jako v pˇr´ıkladˇe 1 z´ıskáme rovnováˇznou polohu kuliˇcky dané hmotnosti zavˇeˇsené nad pˇrekáˇzkou. Protoˇze polohu kuliˇcky z´ıskáme jako ˇreˇsen´ı minimalizaˇcn´ı u ´ lohy s omezen´ım (5.9), pouˇzijeme nyn´ı metodu kvadratické penalty a na rozd´ıl od pˇr´ıkladu 1 budeme minimali33

Funkce norm vypoˇc´ıtá velikost vektoru grad.

117

ˇ sen´ı u Reˇ ´loh m 1 2 3 4 5 10

souˇradnice x1 0 0 0 0 0 0

souˇradnice x2 -0,5 -1 -1,5 -2 -2,5 -5

Tabulka 6.2: Rovnováˇzné polohy volnˇe zavˇeˇsené kuliˇcky o hmotnosti m zovat funkci f (x1 , x2 ) = x21 + x22 + mx2 + 21 ̺kα(x)k2 , kde α1 (x1 , x2 ) = max{−x1 −x2 −1, 0} a α2 (x1 , x2 ) = max{x1 −x2 −3, 0}. Penalizaˇcn´ı parametr ̺ zvol´ıme jako 105 . Pro Newtonovu metodu mus´ıme opˇet v kaˇzdém kroce iteraˇcn´ıho procesu vyˇc´ıslit gradient a Hessián. Gradient funkce f (x1 , x2 ) = x21 + x22 + mx2 + 1 ̺kα(x)k2 vypoˇcteme pomoc´ı kódu: 2 function y = df(x) global m rho = 1e5; f = x(1)^ 2 + x(2)^ 2 + m*x(2); h(1,1) = -x(1) - x(2) - 1; h(2,1) = x(1) - x(2) - 3; alpha = max(zeros(2,1), h); gf = [2*x(1); 2*x(2) + m]; if alpha(1,1) == 0 && alpha(2,1) == 034 ga = zeros(2,1); elseif alpha(1,1) > 0 && alpha(2,1) == 0 ga = alpha(1,1)*[-1;-1]; elseif alpha(1,1) == 0 && alpha(2,1) > 0 V této ˇcásti kódu vypoˇcteme gradient funkce 12 kα(x)k2 a zap´ıˇseme jej do promˇenné ga. Pokud si ˇctenáˇr chce náˇs v´ ypoˇcet ovˇeˇrit, doporuˇcujeme, aby si pˇrepsal tuto funkci bez pouˇzit´ı funkce max. Tzn. je rozumné napsat si funkˇcn´ı pˇredpis této funkce pro tyto pˇr´ıpady: −x1 − x2 − 1 ≤ 0 ∧ x1 − x2 − 3 ≤ 0, −x1 − x2 − 1 > 0 ∧ x1 − x2 − 3 ≤ 0, −x1 − x2 − 1 ≤ 0 ∧ x1 − x2 − 3 > 0 a −x1 − x2 − 1 > 0 ∧ x1 − x2 − 3 > 0 a pro tyto pˇr´ıpady pak spoˇc´ıtat prvn´ı parciáln´ı derivace této funkce. 34

118

ga = alpha(2,1)*[1;-1]; else ga = alpha(1,1)*[-1;-1] + alpha(2,1)*[1;-1]; end; y = gf + rho*ga; Hessián této funkce vypoˇcteme pomoc´ı kódu: function y = ddf(x) global m rho = 1e5; f = x(1)^ 2 + x(2)^ 2 + m*x(2); h(1,1) = -x(1) - x(2) - 1; h(2,1) = x(1) - x(2) - 3; alpha = max(zeros(2,1), h); Hf = [2 0; 0 2]; if alpha(1,1) == 0 && alpha(2,1) == 035 Ha = zeros(2,2); elseif alpha(1,1) > 0 && alpha(2,1) == 0 Ha = [1 1; 1 1]; elseif alpha(1,1) == 0 && alpha(2,1) > 0 Ha = [1 -1; -1 1]; else Ha = [2 0; 0 2]; end; y = Hf + rho*Ha;

V této ˇcásti kódu vypoˇcteme Hessián funkce 21 kα(x)k2 a zap´ıˇseme jej do promˇenné Ha. Pokud si ˇctenáˇr chce náˇs v´ ypoˇcet ovˇeˇrit, doporuˇcujeme, aby si podobnˇe jako v pˇredchoz´ım kódu pˇrepsal tuto funkci bez pouˇzit´ı funkce max. Tzn. opˇet je rozumné napsat si funkˇcn´ı pˇredpis této funkce pro ˇctyˇri pˇr´ıpady uvedené v pˇredchoz´ı rutinˇe a pro tyto pˇr´ıpady pak spoˇc´ıtat druhé parciáln´ı derivace této funkce. 35

119

ˇ sen´ı u Reˇ ´loh Pokud ukonˇcuj´ıc´ı podm´ınku epsilon zvol´ıme jako 10−4 , z´ıskáme pro poˇzadované hmotnosti rovnováˇzné polohy kuliˇcky uvedené v tabulce 6.3.

m 1 2 3 4 5 10

souˇradnice x1 0 0 0,25 0,5 0,75 1

souˇradnice x2 -0,5 -1 -1,25 -1,5 -1,75 -2

Tabulka 6.3: Rovnováˇzné polohy kuliˇcky o hmotnosti m zavˇeˇsené nad pˇrekáˇzkou

Na následuj´ıc´ıch obrázc´ıch jsou zakresleny rovnováˇzné polohy kuliˇcky zavˇeˇsené nad pˇrekáˇzkou pro vˇsechny uvaˇzované hmotnosti.

−0.5

−0.5

x2

0

x2

0

−1

−1.5

−1

−1.5

−2

−2 −0.5

0

0.5

x1

1

1.5

2

−0.5

0

0.5

x1

1

1.5

2

Obrázek 6.24: Rovnováˇzné polohy kuliˇcky o hmotnosti m=1 (vlevo) a m=2 (vpravo) zavˇeˇsené nad pˇrekáˇzkou

120

−0.5

−0.5

x2

0

x2

0

−1

−1.5

−1

−1.5

−2

−2 −0.5

0

0.5

x1

1

1.5

2

−0.5

0

0.5

x1

1

1.5

2


−0.5

−0.5

x2

0

x2

0

−1

−1.5

−1

−1.5

−2

−2 −0.5

0

0.5

x1

1

1.5

2

−0.5

0

0.5

x1

1

1.5

2


121

122

Apendix Apendix A: Inverzn´ı funkce Definice A.1 Bud’ f : R 7→ R. Funkci f −1 , pro niˇz plat´ı souˇcasnˇe: i) definiˇcn´ı obor Df −1 funkce f −1 je roven oboru hodnot funkce f , ii) pro kaˇzdé x ∈ Df −1 plat´ı, ˇze f −1 (x) = y ⇔ f (y) = x, nazveme funkc´ ı inverzn´ ı k funkci f . Pozn´ amka A.1 Lze ukázat, ˇze f −1 existuje právˇe tehdy, je-li f prostá. Graf f −1 je pˇritom osovˇe soumˇerný s grafem f dle pˇr´ımky y = x.

123


Apendix B: Funkce arkussinus

Definice B.1 Funkci inverzn´ı k funkci sinus z´ uˇzené na interval − π2 , π2 nazveme arkussinus a oznaˇcujeme jako arcsin, tj. −1 . arcsin := sin|h− π , π i 2 2

Pozorov´ an´ı B.1 Funkce arkussinus má tyto vlastnosti (viz také n´ıˇze uvedený obrázek): • definiˇcn´ı obor je roven h−1, 1i,

• oborem hodnot je interval − π2 , π2 ,

• funkce arcsin je lichá, tj. arcsin x = arcsin (−x) pro kaˇzdé x ∈ h−1, 1i.

π/2

arcsin x 1

y

sin x

0

−1 −π/2 −π/2

124

−1

0

x

1

π/2

Apendix C: Funkce arkustangens Definice C.1 Funkci inverzn´ı k funkci tangens z´ uˇzené na interval (− π2 , π2 ) nazveme arkustangens a oznaˇcujeme jako arctg, tj. −1 . arctg := tg|(− π , π ) 2 2

Pozorov´ an´ı C.1 Funkce arkustangens má tyto vlastnosti: • definiˇcn´ı obor je roven R, • oborem hodnot je interval (− π2 , π2 ),

• funkce arctg je lichá, tj. arctg x = arctg(−x) pro kaˇzdé x ∈ R.

tg x

π/2

arctg x

y

π/4 0

−π/4 −π/2

− π2

−1

0

x

π

1 2

125


Apendix D: Soustavy line´ arn´ıch rovnic a jejich ˇ reˇsen´ı Definice D.1 Soustavou n line´ arn´ ıch rovnic o n nezn´ am´ ych x1 , . . . , xn nazýváme systém rovnic ve tvaru:  a11 x1 + · · · + a1n xn = b1   .. .. .. .. (7.1) . . . .   an1 x1 + · · · + ann xn = bn

ˇ ısla aij , i = 1, . . . , n, j = 1, . . . , n, nazýváme koeficienty soustavy a ˇc´ısla C´ bi , i = 1, . . . , n, nazýváme pravé strany.

Soustavu (7.1) budeme u ´ spornˇe zapisovat pomoc´ı matice A a vektoru b, kde     a11 · · · a1n b1    ..  , .. A =  ... b =  ...  . . .  an1 · · · ann bn

Matici A nazýváme matic´ ı soustavy (7.1) a vektor b nazýváme pravou stranou soustavy (7.1). Má-li vektor x za sloˇzky neznámé x1 , . . . , xn , m˚ uˇzeme soustavu (7.1) zapsat v maticové podobˇe Ax = b. ˇ sen´ı x soustavy (7.1) m˚ Reˇ uˇzeme naj´ıt pomoc´ı inverzn´ı matice k matici ˇ sen´ı x m˚ A (pojem inverzn´ı matice jsme definovali v kapitole 2.1). Reˇ uˇzeme −1 nalézt tak, ˇze soustavu Ax = b vynásob´ıme zleva matic´ı A a po u ´ pravˇe z´ıskáme x = A−1 b. Protoˇze výpoˇcet inverzn´ı matice je numericky znaˇcnˇe nároˇcný, pouˇz´ıvaj´ı se pro ˇreˇsen´ı velkých soustav jiné metody. Bˇeˇznˇe pouˇz´ıvaná je Gaussova eliminaˇ cn´ ı metoda, která spoˇc´ıvá v nahrazen´ı p˚ uvodn´ı soustavy soustavou, která má stejné ˇreˇsen´ı jako soustava p˚ uvodn´ı, ale toto ˇreˇsen´ı z n´ı snadno vyˇcteme. Pro nalezen´ı nové jednoduˇsˇs´ı soustavy pouˇz´ıváme tzv. ekvivalentn´ı u ´ pravy. Gaussovu eliminaˇcn´ı metodu má ˇrada matematických softwar˚ u jiˇz implementovánu, napˇr. v programu Matlab m˚ uˇzeme ˇreˇsen´ı soustavy Ax = b zapsat pˇr´ıkazem x = A \ b.

126

Literatura [1] J. Bouchala: Matematick´ a analýza 1. ké ˇskoly bán ˇ ské - Technické univerzity http://homel.vsb.cz/∼bou10/MA1/ma1.html

Skripta Ostrava

Vyso(2005).

[2] M. Braun: Differential Equations and Their Applications. Springer Verlag (1993). [3] M. Barnsley: Fractals Everywhere. Academic Press, Boston (1993). [4] Z. Dostál, P. Beremlijski: Metody optimalizace. Text vytvoˇrený pˇri realizaci projektu Matematika pro inˇzenýry 21. stolet´ı“, ” Vysoká ˇskola bán ˇ ská - Technická univerzita Ostrava (2012). http://mi21.vsb.cz/modul/metody-optimalizace [5] Z. Dostál, V. Vondrák: Line´ arn´ı algebra. Text vytvoˇrený pˇri realizaci projektu Matematika pro inˇzenýry 21. stolet´ı“, Vy” soká ˇskola bán ˇ ská - Technická univerzita Ostrava (2011). http://mi21.vsb.cz/modul/linearni-algebra [6] B.G. Katzung: Základn´ı a klinick´ a farmakologie. H & H (1994). [7] T. Kozubek, T. Brzobohatý, V. Hapla, M. Jaroˇsová, A. Markopoulos: Lineárn´ı algebra s Matlabem. Text vytvoˇrený pˇri realizaci projektu Matematika pro inˇzenýry 21. stolet´ı“, Vysoká ˇskola bán ˇ ská - Technická ” univerzita Ostrava (2012). http://mi21.vsb.cz/modul/linearni-algebras-matlabem [8] B. Krajc, P. Beremlijski: Obyˇcejné diferenci´ aln´ı rovnice. Text vytvoˇrený pˇri realizaci projektu Matematika pro inˇzenýry 21. sto” let´ı“, Vysoká ˇskola bán ˇ ská - Technická univerzita Ostrava (2012). http://mi21.vsb.cz/modul/obycejne-diferencialni-rovnice [9] M. Tennenbaum, H. Pollard: Ordinary Differential Equations: An Elementary Textbook for Students of Mathematics, Engineering, and the Sciences. Dover Publications (1985). 127

LITERATURA [10] P. Vodstrˇcil, J. Bouchala: Drobn´ a pˇrekvapen´ı spojen´ a s numerickou integrac´ı. Pokroky matematiky, fyziky a astronomie, svazek 55, ˇc´ıslo 4 (2010). http://dml.cz/dmlcz/141970

128

Fakulta elektrotechniky a informatiky Katedra aplikované matematiky. P. Beremlijski, M. Sadowská. Matematika s radostí

Recommend Documents