ˇ´ıtac ˇova ´ cvic ˇen´ı Poc
Petr Beremlijski, Marie Sadowsk´a
Katedra aplikovan´e matematiky Fakulta elektrotechniky a informatiky ˇ - Technick´a univerzita Ostrava VSB
ˇen´ı 1 - Matlab - na ´ stroj pro matematicke ´ modelova ´ n´ı Cvic Abychom se mohli vˇenovat numerick´emu ˇreˇsen´ı matematick´ ych u ´loh, potˇrebujeme vhodn´e prostˇred´ı, kter´e n´am to umoˇzn´ı. A tak jako fyzik ˇci chemik maj´ı svou laboratoˇr nebo patolog pitevnu, maj´ı i numeriˇct´ı matematici svoj´ı Maticovou laboratoˇr1 - Matlab. Podrobnˇe se tomuto pracovn´ımu prostˇred´ı a jeho pˇr´ıkaz˚ um vˇenuje pˇriloˇzen´ y Matlabovsk´ y slabik´aˇr2 . My si v tomto textu uvedeme pouze struˇcn´ y pˇrehled matlabovsk´ ych promˇenn´ ych a pˇr´ıkaz˚ u, kter´ ym se budeme vˇenovat. Prostˇ red´ı help, demos, intro, who, whos, clear, size, length
Promˇ enn´ e • Skal´ary • Vektory • Matice
Pˇ r´ıkazy • Skal´arn´ı funkce - sin, cos, tan, exp, log, abs, sqrt, round • Vektorov´e funkce a generov´an´ı vektor˚ u - max, min, sort • Maticov´e funkce a generov´an´ı matic - det, rand, ones, zeros, eye • Skal´arn´ı operace - +, −, ∗, /, b
• Maticov´e a vektorov´e operace - +, −, ∗, ´(transponov´an´ı), \ (A\v = x ⇔ Ax = v) Operace “po prvc´ıch” - .∗ , .b, ./
• 2D grafika (vykreslen´ı graf˚ u funkc´ı jedn´e promˇenn´e) - plot, hold on, hold off, figure • 3D grafika (vykreslen´ı graf˚ u funkc´ı dvou promˇenn´ ych) - meshgrid, mesh, contour, hold on, hold off, figure ˇ ıd´ıc´ı pˇr´ıkazy - if (podm´ınˇen´ • R´ y pˇr´ıkaz) for, while (pˇr´ıkazy cyklu se zn´am´ ym poˇctem opakov´an´ı a podm´ınkou na zaˇca´tku) 1 2
MATrix LABoratory K. Sigmon - MATLAB Primer
1
• Relace a logick´e operace - <, >, <=, >=, ==, ˜=, &, |, ˜ • Skripty a funkce - function
Vˇse si vyzkouˇs´ıme pˇri ˇreˇsen´ı n´asleduj´ıc´ıch u ´loh. Pˇ r´ıklad 1 Kolik ˇclen˚ u harmonick´e ˇrady3 mus´ıte nejm´enˇe seˇc´ıst, aby tento ˇca´steˇcn´ y souˇcet ˇrady mˇel hodnotu alespoˇ n 10 (15, 20)? Pˇ r´ıklad 2 Zkuste odhadnout s vyuˇzit´ım Matlabu souˇcet ˇrady
P∞
n=1
1 n
·
1 n+1
.
Pˇ r´ıklad 3 Sestrojte grafy n´asleduj´ıc´ıch funkc´ı: • f (x) = x2 √ • f (x) = 1 − x2 • f (x) = x2 · sin • f (x) = |x|
1 x2
Pˇ r´ıklad 4 Sestrojte grafy n´asleduj´ıc´ıch funkc´ı: • f (x, y) = x2 + y 2 p • f (x, y) = x2 + y 2 • f (x, y) = (x2 + y 2 ) · sin • f (x, y) =
p
1 x2 +y 2
|xy|
3ˇ
Radou (re´aln´ ych ˇc´ısel) rozum´ıme v´ yraz a1 + a2 + · · · + an + · · · = P ∞ an ∈ R. Harmonickou naz´ yv´ame ˇradu n=1 n1 .
2
P∞
n=1
an , kde pro kaˇzd´e n ∈ N je
ˇen´ı 2 - Co doka ´ˇ ´ˇ Cvic ze a nedoka ze Matlab ˇeˇ ´u ´ lohy s permutac´ı -r sen´ı jedne Programy pro numerick´e v´ ypoˇcty, mezi nˇeˇz patˇr´ı i Matlab, um´ı vyˇreˇsit mnoho i velmi komplikovan´ ych u ´loh v “rozumn´em” ˇcase. Nejsou vˇsak vˇsemocn´e a maj´ı sv´e limity. Vhodnˇe to ilustruje tato u ´loha: Naleznˇete nejmenˇs´ı pˇrirozen´e ˇc´ıslo n, pro kter´e nejvˇetˇs´ı spoleˇcn´ y dˇelitel ˇc´ısel (n17 + 9) a ((n + 1)17 + 9) nen´ı 1. Pokud byste chtˇeli tuto u ´lohu ˇreˇsit poˇc´ıtaˇcem, ˇcekali byste opravdu dlouho, protoˇze ˇreˇsen´ım t´eto u ´lohy je n = 84244329255928893292881973223089006724594204607924334 . My si “omezen´ı” Matlabu uvˇedom´ıme d´ıky n´asleduj´ıc´ımu pˇr´ıkladu. V tomto pˇr´ıkladˇe budeme m´ıt za u ´kol zjistit pravdˇepodobnost v´ yskytu dan´eho jevu. Nejjednoduˇsˇs´ı, ale ne nejm´enˇe pracnou, moˇznost´ı, jak vypoˇc´ıst pravdˇepodobnost urˇcit´eho jevu, je zjistit poˇcet pˇr´ızniv´ ych moˇznost´ı a podˇelit jej poˇctem vˇsech moˇznost´ı. Pˇ r´ıklad 1 Mˇejme posloupnost pˇrirozen´ ych ˇc´ısel 1, 2, 3, . . . , n. Pot´e ji n´ahodnˇe prom´ıchejte. Jak´a je pravdˇepodobnost, ˇze ani jedno ˇc´ıslo nebude na sv´e p˚ uvodn´ı pozici?5 Vyˇreˇste pomoc´ı Matlabu. Pro jak velk´e n jste schopni tuto pravdˇepodobnost zjistit? N´apovˇeda: Pro algoritmizaci u ´lohy pouˇzijte rekurzi6 .
4
Odhadnˇete, jak dlouho budete ˇcekat na nalezen´ı ˇreˇsen´ı, pokud m´ate poˇc´ıtaˇc, kter´ y zvl´adne miliardu operac´ı za sekundu a pˇredpokl´ad´ate, ˇze nalezen´ı nejvˇetˇs´ıho spoleˇcn´eho dˇelitele pro jedno konkr´etn´ı n se d´a dos´ahnout jednou operac´ı. 5 Vˇsimnˇete si, jak´e ˇc´ıslo z´ısk´ate, pokud vypoˇctete pˇrevr´acenou hodnotu zjiˇstˇen´e pravdˇepodobnosti. 6 V programov´an´ı rekurze pˇredstavuje opakovan´e vnoˇren´e vol´an´ı stejn´e funkce (podprogramu), takovou funkci pak nazveme rekurzivn´ı. Souˇca´st´ı rekurzivn´ı funkce mus´ı b´ yt ukonˇcuj´ıc´ı podm´ınka, kter´a urˇcuje, kdy se m´a vnoˇrov´an´ı zastavit. Pro pouˇzit´ı rekurze je zapotˇreb´ı, aby programovac´ı jazyk umoˇzn ˇoval vol´an´ı podprogramu jeˇstˇe pˇred ukonˇcen´ım jeho pˇredchoz´ıho vol´an´ı. Po kaˇzd´em kroku vol´an´ı sebe sama doch´az´ı ke zjednoduˇsen´ı probl´emu, pokud nen´ı splnˇena ukonˇcuj´ıc´ı podm´ınka, provede se rekurzivn´ı krok.
3
ˇen´ı 3 - Pra ´ ce s vektory a maticemi Cvic ´ va ´ n´ı v databa ´ z´ıch a trocha geometrie navrch - vyhleda Pod´ıvejme se na dvˇe aplikace aritmetick´ ych vektor˚ u a matic. ´ va ´ n´ı v databa ´ z´ıch Vyhleda Nejprve si zavedeme pojem aritmetick´ y vektor. N-rozmˇ ern´ y aritmetick´ y vektor je uspoˇra´dan´a n-tice ˇc´ısel, jej´ıˇz prvky se naz´ yvaj´ı sloˇzky. Tyto uspoˇra´dan´e n-tice budeme zapisovat do hranat´ ych z´avorek do ˇra´dk˚ u nebo sloupc˚ u. Pˇripomeˇ nme si skal´ arn´ı souˇ cin dvou aritmetick´ ych vektor˚ u o n sloˇzk´ach, kter´ y je definov´an vztahem
kde ||u|| =
p Pn
i=1
u · v = ||u|| ||v|| cos ϕ, u2i a ϕ je u ´hel, kter´ y sv´ıraj´ı vektory u a v. u
PSfrag replacements v
ϕ
O vektorech ˇrekneme, ˇze jsou si “bl´ızk´e”, pokud u ´hel ϕ je bl´ızk´ y 0. Napˇr. na n´asleduj´ıc´ım obr´azku je vektor u “bl´ızk´ y” v. u
PSfrag replacements
v
ϕ1 ϕ2
w
“Bl´ızkosti” vektor˚ u m˚ uˇzeme vyuˇz´ıt pˇri vyhled´av´an´ı v datab´az´ıch. Ilustrujme si to na n´asleduj´ıc´ım pˇr´ıkladˇe. Mˇejme kuchaˇrku, kter´a obsahuje n recept˚ u, kde kaˇzd´ y recept obsahuje nˇekterou z k surovin. Takovou datab´azi recept˚ u m˚ uˇzeme popsat n vektory r1 , r2 , . . . , rn o k sloˇzk´ach. Pokud chceme vybrat recept se surovinami, kter´e nejl´epe splˇ nuj´ı n´aˇs poˇzadovan´ y v´ ybˇer, pak m˚ uˇzeme n´aˇs poˇzadavek zapsat jako vektor p o k sloˇzk´ach a naˇs´ım u ´kolem je naj´ıt vektor ri z vektor˚ u r1 , r2 , . . . , rn , kter´ y je “nejbliˇzˇs´ı” vektoru p (tj. plat´ı, ˇze ri ·p rl ·p = cos ϕi je nejvˇetˇs´ı ze vˇsech ||rl || ||p|| = cos ϕl ∀l = 1, . . . , n). ||ri || ||p||
4
Pˇ r´ıklad 1 Popiˇste n´asleduj´ıc´ı recepty pomoc´ı vektor˚ u jako v pˇredchoz´ım textu: • 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´ • 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´ • 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, 10g m´ asla, 2 lˇ z´ıce plnotuˇ cn´ e hoˇrˇ cice, 2 vejce, petrˇ zelka • 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 • Ev´ıkova mˇ namka - ingredience: 4 kuˇrec´ı ˇr´ızky, 8 pl´ atk˚ u anglick´ e slaniny, s´ yr podle chuti (eidam, blat’a ´ck´ 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 • 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´ı • Chleb´ıˇ ckov´ a pomaz´ anka - ingredience: 100 g brambor, 1 cibule, 1 lˇ z´ıce tatarsk´ e om´ aˇ cky, s˚ ul, pepˇr • 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 • 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 • 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 • 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 • 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´ı • 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 • 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´ı • 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) • 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´ı • Vepˇrov´ a k´ yta - ingredience: 6 ˇr´ızk˚ u z 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
Zkuste na z´akladˇe “bl´ızkosti” vyhledat, kter´e recepty obsahuj´ı ingredience nejl´epe odpov´ıdaj´ıc´ı tˇemto tˇrem poˇzadavk˚ um: • 100 gram˚ u cukru, 20 gram˚ u droˇzd´ı, 100 gram˚ u Hery, 2 lˇz´ıce kakaa, 500 gram˚ u mouky a skoˇrice 5
• s˚ ul, s´ yr a worchester • 2 strouˇzky ˇcesneku, 1 lˇz´ıce hoˇrˇcice, 3 kuˇrec´ı ˇr´ızky, major´anka, pepˇr, 2 pl´atky slaniny a s˚ ul
ˇ i za ´ pisu geometricky ´ ch zobrazen´ı Pouˇ zit´ı matic pr Zaˇcneme zaveden´ım pojmu matice. Necht’ jsou d´any prvky a11 , a12 , . . . , amn z dan´e mnoˇziny F. Matice typu (m, n) je obd´eln´ıkov´a tabulka a11 . . . a1n .. , ... A = ... . am1 . . . amn kter´a m´a m · n prvk˚ u aij uspoˇra´dan´ ych do m ˇra´dk˚ u riA a n sloupc˚ u sA ze j , takˇ r1A A A = ... = sA 1 , . . . , sn , A rm
a1j .. riA = [ai1 , . . . , ain ] , sA . . j = amj
Struˇcnˇe p´ıˇseme t´eˇz A = [aij ]. Nyn´ı si zavedeme operaci n´asoben´ı matic. K tomu ale potˇrebujeme nejdˇr´ıve definovat n´asoben´ı matice a vektoru. Souˇcinem matice A = [aij ] typu (m, n) a sloupcov´eho vektoru x = [xi ] o n sloˇzk´ach naz´ yv´ame vektor y o m sloˇzk´ach definovan´ y pˇredpisem A y = Ax = x1 sA 1 + · · · + x n sn .
Rozeps´an´ım definice po sloˇzk´ach dostaneme [y]i = [Ax]i = ai1 x1 + · · · + ain xn = riA x. Ted’ m˚ uˇzeme pˇrej´ıt k n´asoben´ı matic. Jestliˇze A je matice typu (m, p) a B je matice typu (p, n), pak souˇ cin matic A a B je matice AB typu (m, n) definovan´a pˇredpisem B AB = AsB 1 , . . . , Asn .
Rozep´ıˇseme-li si definici n´asoben´ı matic po sloˇzk´ach, dostaneme
[AB]ij = ai1 b1j + · · · + aip bpj = riA · sB j . 6
Neˇz pouˇzijeme v´ yˇse zavedenou operaci, pˇripomeneme si nˇekter´e pojmy z oblasti geometrick´ ych zobrazen´ı. Geometrick´ ym zobrazen´ım nazveme zobrazen´ı, kter´e kaˇzd´emu bodu A u ´tvaru U pˇriˇrazuje pr´avˇe jeden bod A0 u ´tvaru U 0 . Bod A je tzv. vzor a bod A0 se oznaˇcuje jako obraz. My se budeme zab´ yvat tˇremi zobrazen´ımi, a to stejnolehlost´ı, rotac´ı a posunut´ım. Pro pˇripomenut´ı: Stejnolehlost: Mˇejme v rovinˇe ˇci prostoru ρ bod S. Geometrick´e zobrazen´ı, pˇri nˇemˇz obrazem bodu S je bod S a obrazem kaˇzd´eho A ∈ ρ, A 6= S, je takov´ y bod A0 ∈ ρ, ˇze pro vektor SA0 plat´ı SA0 = κSA, kde κ ∈ R\{0, 1} je pevnˇe zvolen´e, se naz´ yv´a stejnolehlost´ı (homoteti´ı). Bod S naz´ yv´ame stˇredem stejnolehosti a κ koeficientem stejnolehlosti. Stejnolehlost s κ = −1 je stˇredovou soumˇernost´ı. y
PSfrag replacements y
y
stejnolehlost s κ = 2, S = [0, 0] x
PSfrag replacements
vzor zobrazen´ı
vzor zobrazen´ı x
x
Rotace: Otoˇcen´ı (rotace) v rovinˇe je geometrick´e zobrazen´ı, kter´e je charakterizov´ano t´ım, ˇze spojnice vˇsech bod˚ u s pevnˇe zvolen´ ym bodem S, tzv. stˇredem otoˇcen´ı, se zmˇen´ı o stejn´ yu ´hel ϕ, tzv. u ´hel otoˇcen´ı, a vzd´alenost bod˚ u od stˇredu otoˇcen´ı z˚ ust´av´a nezmˇenˇena. PSfrag replacements y
y
PSfrag replacements
y
x vzor zobrazen´ı
rotace s ϕ = π/4, S = [0, 0]
vzor zobrazen´ı x
x
Posunut´ı: Posunut´ı (translace) je geometrick´e zobrazen´ı, kter´e je charakterizov´ano t´ım, ˇze vˇsechny body transformovan´e mnoˇziny bod˚ u zmˇen´ı sv´e kart´ezsk´e souˇradnice o stejnou hodnotu, tj. 7
PSfrag replacements y
y
PSfrag replacements
y
x vzor zobrazen´ı
posunut´ı s vektorem posunut´ı [1, 1]
vzor zobrazen´ı x
x
ke kaˇzd´emu bodu pˇriˇcteme stejn´ y vektor posunut´ı. Nyn´ı si uk´aˇzeme, jak lze pomoc´ı element´arn´ıch maticov´ ych operac´ı zapsat v´ yˇse uveden´a zobrazen´ı. Zapiˇsme souˇradnice vzoru geometrick´eho zobrazen´ı, tj. bodu z R2 ˇci R3 , do sloupcov´eho vektoru. Pokud vzorem zobrazen´ı je v´ıce bod˚ u, zap´ıˇseme je jako sloupcov´e vektory do matice P . Obraz bod˚ u ve stejnolehlosti s koeficientem κ a stˇredem v poˇca´tku zap´ıˇseme jako souˇcin transformaˇcn´ı matice T typu (n, n), kde n je rozmˇer prostoru, ve kter´em zobrazujeme (2 nebo 3), a matice P . Matice T m´a vˇsechny prvky na hlavn´ı diagon´ale rovny koeficientu stejnolehlosti κ. Vˇsechny dalˇs´ı prvky jsou nulov´e. Obraz bod˚ u v rotaci s u ´hlem otoˇcen´ı ϕ a stˇredem otoˇcen´ı v poˇca´tku zap´ıˇseme jako souˇcin transformaˇcn´ı matice R typu (2, 2) a matice P , kde cos ϕ − sin ϕ R= . sin ϕ cos ϕ Posunut´ı s vektorem posunut´ı pos realizujeme tak, ˇze ke kaˇzd´emu sloupci matice P pˇriˇcteme sloupcov´ y vektor pos. Pˇ r´ıklad 2 Implementujte geometrick´a zobrazen´ı - stejnolehlost, rotaci a posunut´ı. Naleznˇete obraz troj´ uheln´ıku s vrcholy [1, −1], [2, 0] a [1, 1] ve stejnolehlosti se stˇredem v poˇca´tku a s koeficientem stejnolehlosti 2, v rotaci se stˇredem otoˇcen´ı v poˇca´tku a u ´hlem otoˇcen´ı π/2 a v posunut´ı s vektorem posunut´ı [1, −1]. Zobrazte vzor a jednotliv´e obrazy (pouˇzijte funkce plot a patch).
8
ˇen´ı 4 - Monte Carlo Cvic Hledejme obsah kruhu o polomˇeru r (pˇredstavme si, ˇze nezn´ame pˇr´ısluˇsn´ y vzorec a nic nev´ıme o ˇc´ıslu π). Prvn´ım n´apadem by mohlo b´ yt zhotoven´ı v´alcov´ ych n´adob s r˚ uzn´ ymi polomˇery podstav (napˇr. s polomˇery o d´elce 1 a 2 jednotek) a jednotkovou v´ yˇskou.
v=1
v=1
r=1
r=2
Objem vody, kter´ y se do takov´ ych n´adob vejde, je roven obsahu podstavy v´alce, tj. obsahu kruhu s polomˇerem r. Rychle si vˇsimneme, ˇze pokud zvˇetˇs´ıme polomˇer podstavy dvakr´at, zvˇetˇs´ı se objem ˇctyˇrikr´at a n´aslednˇe odvod´ıme, ˇze obsah kruhu je pˇr´ımo u ´mˇern´ y druh´e mocninˇe polomˇeru. Tak´e zjist´ıme, ˇze druhou mocninu polomˇeru kruhu mus´ıme vyn´asobit vhodnou konstantou, abychom dostali spr´avnou hodnotu obsahu dan´eho kruhu. Tuto konstantu oznaˇc´ıme π. Existuje mnoho moˇznost´ı, jak tuto konstantu odhadnout. Snad nejjednoduˇsˇs´ı zp˚ usob, jak stanovit meze pro π, je vepsat do kruhu o polomˇeru r ˇctverec a stejn´emu kruhu opsat jin´ y ˇctverec.
r
r
Protoˇze je snadn´e √ spoˇc´ıtat obsahy dan´ ych ˇctverc˚ u, zjist´ıme, ˇze π ∈ (2, 4) (d´elka strany menˇs´ıho ˇctverce je 2r, d´elka strany vˇetˇs´ıho ˇctverce je 2r). Mnohem rozumnˇejˇs´ı v´ ysledek z´ısk´ame, pokud budeme dan´emu kruhu vepisovat n-´ uheln´ıky a poˇc´ıtat jejich obsahy. Tuto metodu naz´ yv´ame vyˇcerp´avac´ı (exhaustn´ı) a pravdˇepodobnˇe prvn´ı ji pouˇzil Eudoxos 7 . Neˇz se budeme vˇenovat odhad˚ um ˇc´ısla π, pod´ıvejme se kr´atce na v´ ypoˇcet obvodu kruhu. Jistˇe v´ıme, ˇze obvod kruhu je pˇr´ımo u ´mˇern´ y dvojn´asobku jeho polomˇeru, ale abychom dostali spr´avnou hodnotu, je nutno 2r vyn´asobit vhodnou konstantou. Na n´asleduj´ıc´ım obr´azku provedeme pˇreuspoˇra´d´an´ı kruhu na u ´tvar, kter´ y se pro zjemˇ nuj´ıc´ı se dˇelen´ı kruhu bl´ıˇz´ı obd´eln´ıku. Porovn´an´ım obsahu kruhu a vznikl´eho obd´eln´ıku je n´azornˇe vidˇet, ˇze tato konstanta je opˇet π. 7
Eudoxos (410 nebo 408 pˇr. n. l. – 355 nebo 347 pˇr. n. l.) – ˇreck´ y astronom, matematik a fyzik, student Plat´ona
9
o = 2kr r
PSfrag replacements r r
r
k·r
k·r
zjemnˇen´ı dˇelen´ı
r
r k·r k=π
Nyn´ı si uk´aˇzeme nˇekolik zp˚ usob˚ u, jak nal´ezt pˇribliˇznou hodnotu ˇc´ısla π. Buffonova metoda ˇ sen´ım tzv. Buffonova probl´emu s jehlou8 je aproximace ˇc´ısla π. Uloha ´ Reˇ spoˇc´ıv´a v opakovan´em h´azen´ı jehly o d´elce l na rovinu, na kter´e m´ame vyznaˇcenu s´ıt’ rovnobˇeˇzek o vzd´alenosti s ≥ l. Jestliˇze jehlu hod´ıme n-kr´at a x-kr´at n´am bˇehem tˇechto pokus˚ u po dopadu zkˇr´ıˇz´ı nˇekterou z rovnobˇeˇzek, pak v pˇr´ıpadˇe, ˇze s = 2l, ˇc´ıslo n x aproximuje s danou pˇresnost´ı ˇc´ıslo π. V roce 1975 p´anov´e Perlman a Wichera publikovali tento v´ ysledek t´ ykaj´ıc´ı se pˇresnosti Buffonovy metody: S pravdˇ e podobnost´ ı 95 procent nem´ a chyba aproximace hodnotu vˇetˇs´ı √ neˇz 5/ n. Tzn. napˇr´ıklad pro 10000 pokus˚ u n´am s pravdˇepodobnost´ı 95 procent chyba nepˇrekroˇc´ı hodnotu 0, 05. 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.
8
G. L. Buffon (1707–1788) – francouzsk´ y pˇr´ırodovˇedec
10
Metoda Monte Carlo Monte Carlo je tˇr´ıda v´ ypoˇcetn´ıch algoritm˚ u zaloˇzen´a na prov´adˇen´ı n´ahodn´ ych experiment˚ u. T´eto metody se ˇcasto pouˇz´ıv´a pro simulaci fyzik´aln´ıch a matematick´ ych syst´em˚ u. V´ ysledkem proveden´ı velk´eho mnoˇzstv´ı experiment˚ u je obvykle pravdˇepodobnost urˇcit´eho jevu. Na z´akladˇe z´ıskan´e pravdˇepodobnosti a zn´am´ ych vztah˚ u pak spoˇc´ıt´ame potˇrebn´e v´ ysledky. Protoˇze metoda vyˇzaduje generov´an´ı velk´eho souboru n´ahodn´ ych dat, je vhodn´e pro jej´ı implementaci pouˇzit´ı poˇc´ıtaˇce. Metod Monte Carlo se pouˇz´ıv´a v pˇr´ıpadˇe, kdy je pˇr´ıliˇs pracn´e nebo nemoˇzn´e nal´ezt pˇresn´ y v´ ysledek jin´ ym zp˚ usobem. Jej´ı v´ yhodou je jednoduch´a implementace, nev´ yhodou relativnˇe mal´a pˇresnost. Metoda byla vytvoˇrena skupinou fyzik˚ u pracuj´ıc´ıch na projektu jadern´e pumy v Los Alamos, jm´eno metody bylo navrˇzeno v roce 1940 von Neumannem9 . V matematice se Monte Carlo pouˇz´ıv´a zejm´ena pro v´ ypoˇcet urˇcit´ ych integr´al˚ u (zejm´ena v´ıcen´asobn´ ych urˇcit´ ych integr´al˚ u), kter´e je obt´ıˇzn´e ˇci nemoˇzn´e vyˇc´ıslit analyticky nebo jinou vhodnou numerickou metodou. R 1 2 Napˇr. obsah plochy ohraniˇcen´e shora grafem funkce 2 y = x na intervalu h0, 1i (tj. 0 x dx) je moˇzn´e metodou Monte Carlo vypoˇc´ıst n´asleduj´ıc´ım zp˚ usobem. Necht’ n´aˇs program generuje n´ahodnˇe dvojice ˇc´ısel [x, y], pˇriˇcemˇz kaˇzd´e z ˇc´ısel x a y je vybr´ano nez´avisle z intervalu h0, 1i. Tuto dvojici budeme ch´apat jako souˇradnice bodu, kter´ y je n´ahodnˇe zvolen ve ˇctverci h0, 1i × h0, 1i. Pravdˇepodobnost toho, ˇze bod leˇz´ı uvnitˇr zadan´eho ˇctverce, je 1. Pravdˇepodobnost toho, ˇze bod leˇz´ı uvnitˇr podmnoˇziny E jednotkov´eho ˇctverce (tj. ˇctverce, jehoˇz strana m´a jednotkovou velikost), je rovna obsahu plochy E. Takˇze obsah plochy, kter´a je podmnoˇzinou jednotkov´eho ˇctverce, m˚ uˇzeme odhadnout jako pravdˇepodobnost, ˇze n´ahodnˇe zvolen´ y bod leˇz´ı v t´eto podmnoˇzinˇe. y y = x2
PSfrag replacements E x
Pˇ r´ıklad 2 Implementujte metodu Monte Carlo pro v´ ypoˇcet
R1 0
x2 dx.
Pro odhad pˇresnoti metody Monte √ Carlo plat´ı: S pravdˇepodobnost´ı 95 procent nem´a u n´am s pravchyba aproximace hodnotu vˇetˇs´ı neˇz 1/ n. Tzn. napˇr´ıklad pro 10000 pokus˚ dˇepodobnost´ı 95 procent chyba nepˇrekroˇc´ı hodnotu 0, 01. R1√ Pokud chceme pouˇz´ıt metodu Monte Carlo k aproximaci ˇc´ısla π, vypoˇcteme 0 1 − x2 dx. Snadno si uvˇedom´ıme, ˇze t´ımto zp˚ usobem z´ısk´ame aproximaci hodnoty π/4. 9
John von Neumann (1903–1957) – mad’arsk´ y matematik
11
y
y=
PSfrag replacements
√ 1 − x2
E x
R1√ Pˇ r´ıklad 3 Implementujte metodu Monte Carlo pro v´ ypoˇcet 0 1 − x2 dx a pouˇzijte ji k aproximaci ˇc´ısla π. Porovnejte vaˇsi aproximaci se skuteˇcnou hodnotou ˇc´ısla π a urˇcete chybu aproximace. ˇ Rady Posledn´ı moˇznost´ı k aproximaci ˇc´ısla π, kterou si v tomto pˇrehledu uk´aˇzeme, je vyuˇzit´ı ˇrady10 . K t´eto aproximaci pouˇzijeme Gregoryho ˇradu11 , kter´a rozv´ıj´ı funkci arctg x: x3 x5 x7 arctg x = x − + − + ... 3 5 7 Tato ˇrada m´a koneˇcn´ y souˇcet pro x ∈ h−1, 1i (ˇr´ık´ame, ˇze konverguje), nav´ıc plat´ı, ˇze ˇc´ım je |x| menˇs´ı neˇz 1, t´ım m´enˇe ˇclen˚ u ˇrady potˇrebujeme pouˇz´ıt k nahrazen´ı arctg x s “uspokojivou” pˇresnost´ı. Prvn´ı moˇznost´ı, jak aproximovat π, je tud´ıˇz nahradit arctg 1 Gregoryho ˇradou (arctg 1 = π/4). Aproximaci s rychlejˇs´ı konvergenc´ı z´ısk´ame, pokud pouˇzijeme rovnost arctg 1 = arctg (1/2) + arctg (1/3).
10 ˇ
Radou (re´aln´ ych ˇc´ısel) rozum´ıme v´ yraz a1 + a2 + · · · + an + · · · = an ∈ R. 11 James Gregory (1638 – 1675) – skotsk´ y matematik a astronom
12
P∞
n=1
an , kde pro kaˇzd´e n ∈ N je
