Deterministické modely Lenka Přibylová 15. února 2019 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Obsah Model a jeho tvorba 6 Statické modely a komparativní statika 33 Statické modely interakcí a teorie her 46 Dynamické modely 59 Rovnovážná dynamika 64 Základní spojité modely růstu 73 Nerovnovážná dynamika 95 Strukturovaný spojitý dynamický model 100 Spojitá a diskrétní dynamika v Rm . 112 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Lineární diskrétní model v rovině 121 Strukturovaný diskrétní dynamický lineární model populace 127 Nelineární dynamika a linearizace 134 Dynamické modely v rovině 140 Dynamika chemických reakcí 154 Dynamické modely interakcí 162 Evoluční hry 186 Teorie her a dynamika 196 Dynamický model difúze a šíření 202 Model difúze s advekcí 234 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Reakčně-difúzní model 241 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model a jeho tvorba Model je zjednodušená reprezentace reálného objektu nebo systému reálných objektů zapsaná rovnicemi nebo počítačovým programem. Deterministickým modelem rozumíme model, ve kterém popisované veličiny spontánně nemění svůj stav a jsou vázány pevně danými vztahy. Stav tedy není náhodná veličina, ale veličina deterministicky určená vztahy, počátečními podmínkami, okrajovými podmínkami apod. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Determinismus je přesvědčení, že vývoj světa je předem dán jeho současným stavem (případně jeho stavem v kterémkoliv bodě v minulosti či na počátku) a absolutně platnými přírodními zákony. Dle tohoto přesvědčení neexistují skutečně náhodné (stochastické) jevy, pocit náhodnosti je dán pouze naší neznalostí příčin. Determinismus dle některých interpretací vylučuje existenci svobodné vůle (inkompatibilismus), jejich slučitelnost je ale možná v podobě dualismu (kompatibilismus). Deterministické přesvědčení bylo silné v 18. a 19. století po objevech mnohých přírodních, zvláště fyzikálních, zákonů. Po objevu kvantové fyziky vliv determinismu mezi vědci zeslábl, přestože ve vědě zesláblo i přesvědčení o svobodné vůli. Tolik z Wikipedie... Mnoho lidí determinismus chápe právě tímto způsobem. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Mou snahou bude předložit poněkud komplexnější pohled. Determinismus v moderním pojetí není v rozporu, ale často jde ruku v ruce se stochastickými jevy a může je dokonce vysvětlovat. Myšlenky tohoto pojetí světa vyslovil poprvé Ilya Prigogine v 70. letech minulého století a ovlivnil tak celou moderní vědu, zvláště oblasti chemické a biochemické, fyzikální, např. právě kvantovou mechaniku, ale ovlivnil i sociální vědy. V pozadí jeho úvah stojí nelineární dynamické jevy, bifurkace a nerovnovážná dynamika. Jedním z úžasných důsledků takového pojetí světa je vysvětlení vzniku řádu z chaosu, vzniku složitých struktur v případě, že je systém vzdálen od své rovnováhy. Takový systém je možný pouze v případě, že si vyměňuje energii nebo informace s okolím, tedy není izolovaný. Izolované systémy spějí nenávratně k rovnováze, stavu s maximální entropií. Věci se rozpadají, káva chladne. Interakce s okolím a výměna energie způsobuje vznik složitých struktur, nerovnovážných avšak organizovaných dějů. Možná i život. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Tvorba modelu: Účel modelu Je realizovatelný? Konstrukce modelu Zhodnocení modelu Akceptovat model Revidovat model Zamítnout model a začít znovu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Účel modelu Deskriptivní a prediktivní mo- dely: Explikativní modely: Hlavní účel: management, tvorba plánu, predikce. Hlavní účel: porozumění principům, rozvoj teorie. Důležitá je numerická přesnost i na úkor jednoduchosti. Numerická přesnost není podstatná, model popisuje princip a má být co nejjednodušší. Některé procesy můžeme ignorovat, pokud nejsou numericky podstatné. Některé procesy můžeme ignorovat, pokud nejsou principiálně podstatné. Předpoklady jsou kvantitativní. Předpoklady jsou kvalitativní. Model je tvořen ”na míru”. Model je aplikovatelný na širokou oblast. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Je model realizovatelný? Nejčastější omezující podmínky jsou • čas - náročnost odhadujte spíše pesimisticky, je lépe začít s jednoduchým modelem a ten pak rozšiřovat • data - zde naopak uvažujte spíše optimisticky, mnohdy nejsou některé parametry modelu třeba, lze je obejít nebo nejsou podstatné • kapacita a výkon počítače - pokud nezpracováváte zrovna kvantitativní model počasí nebo množství hmoty ve Vesmíru můžete být klidní ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Konstrukce modelu: koncepce → diagram → rovnice → počítačová realizace ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Koncepce modelu • Které proměnné jsou pro model podstatné? • Které z nich budou stavové proměnné a které exogenníproměnné a parametry. Stavovou proměnnou je proměnná, která určuje stav popisovaného systému, exogenní proměnnou je obecně funkce nezávislá na stavu systému, parametrem je konstanta nezávislá na stavu systému. Často je lépe začít s více stavovými proměnnými, které během tvory modelu přesouváme mezi exogenní proměnné a parametry. • Jak detailní bude model? Je třeba rozhodnout, které veličiny budeme považovat za identické. Při volbě velké agregace může dojít k chybám, pokud se popisované veličiny chovají odlišným způsobem, pak je třeba je rozdělit, tzv. strukturovat (druhově, věkově apod.), mluvíme pak o strukturovaném modelu. Mnohdy i přes velkou agregaci je nestrukturovaný model vzhledem jeho účelu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × vhodný. Stejně tak přílišný detail vede k přílišné složitosti modelu a mnoha parametrům, které je třeba odhadovat z mnoha dat - a ta nemusí být k dispozici. Musíme vhodně volit mezi chybou danou modelem a chybou danou parametry. Diagram • Zvolené (pro model podstatné) proměnné ”uložíme do krabiček”. • Zakreslíme vzájemné vztahy, které nám pomohou rozhodnout, zda je daná proměnná stavová nebo exogenní, nebo ji můžeme považovat za parametr. • Zakreslení vztahů je první kontrolou vhodné volby agregace. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice • V prvé řadě volíme mezi statickým a dynamickým modelem. Pokud je účelem modelu najít rovnováhu systému bez ohledu na to, jakým způsobem (a zda vůbec!) se tato rovnováha ustanoví, je možné volit model statický. V opačném případě je nutné použít dynamické rovnice. • Je třeba rozhodnout o typu dynamických rovnic. Základním vodítkem je diskrétní resp. spojitý běh času. V diskrétním případě je vhodné použít diferenční rovnice, ve spojitém diferenciální rovnice. Můžeme použít ODR, PDR, rovnice se zpožděním apod. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × • Je třeba rozhodnout, zda budou procesy mezi stavovými proměnnými záviset na jedné nebo více proměnných a parametrech a jak, zda můžeme např. míry těchto procesů považovat za parametry a exogenní proměnné (tedy konstanty nebo funkce nezávislé na stavových proměnných) nebo zda závisí také na stavových proměnných. Složitost v popisu vztahů mezi stavovými veličinami v modelu může být dána jednak samotnými principy nebo také např. nelineárními odhady z naměřených dat. • Rovnice v modelu musí ”sedět” jednotkově. V okamžiku, kdy máme sestaveny rovnice, můžeme je zjednodušit co se týče počtu parametrů vhodnou transformací času a stavových proměnných (nondimensionalization - zbavení se jednotek). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Nondimenzionalizace: N′ = rN 1 − N K , diferenciální rovnice popisující populaci bakterií Zavedením nové stavové proměnné x = N K dx dt = dN dt · 1 K = rN 1 − N K · 1 K = rx(1 − x) a nového času τ = rt dostaneme dx dτ = dx dt · dt dτ = rx(1 − x) · 1 r = x(1 − x). Nová stavová proměnná x je bezrozměrná a představuje míru zaplnění Petriho misky, x = 1 je 100% zaplnění Petriho misky do její kapacity. Časová jednotka τ je vůči jednotce t zkrácena nebo prodloužena tak, aby za ni došlo ke zdvojnásobení počtu bakterií v misce. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Nondimenzionalizace: N′ = rN 1 − N K , diferenciální rovnice popisující populaci bakterií • N hustota populace (např. počet miliónů bakterií v mililitru) • r > 0 je specifická míra růstu (veličina daná poměrem nově vzniklých bakterií ku stávajícím za časovou jednotku na počátku expe- rimentu) • K kapacita prostředí (maximální množství miliónů bakterií, které prostředí uživí - např. Petriho miska) • N′ = dN dt je změna počtu bakterií za časovou jednotku Jednotky odpovídají. Zavedením nové stavové proměnné x = N K dx dt = dN dt · 1 K = rN 1 − N K · 1 K = rx(1 − x) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Nondimenzionalizace: N′ = rN 1 − N K , diferenciální rovnice popisující populaci bakterií Zavedením nové stavové proměnné x = N K dx dt = dN dt · 1 K = rN 1 − N K · 1 K = rx(1 − x) a nového času τ = rt dostaneme dx dτ = dx dt · dt dτ = rx(1 − x) · 1 r = x(1 − x). Nová stavová proměnná x je bezrozměrná a představuje míru zaplnění Petriho misky, x = 1 je 100% zaplnění Petriho misky do její kapacity. Časová jednotka τ je vůči jednotce t zkrácena nebo prodloužena tak, aby za ni došlo ke zdvojnásobení počtu bakterií v misce. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Nondimenzionalizace: N′ = rN 1 − N K , diferenciální rovnice popisující populaci bakterií Zavedením nové stavové proměnné x = N K dx dt = dN dt · 1 K = rN 1 − N K · 1 K = rx(1 − x) a nového času τ = rt dostaneme dx dτ = dx dt · dt dτ = rx(1 − x) · 1 r = x(1 − x). Nová stavová proměnná x je bezrozměrná a představuje míru zaplnění Petriho misky, x = 1 je 100% zaplnění Petriho misky do její kapacity. Časová jednotka τ je vůči jednotce t zkrácena nebo prodloužena tak, aby za ni došlo ke zdvojnásobení počtu bakterií v misce. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Nondimenzionalizace: N′ = rN 1 − N K , diferenciální rovnice popisující populaci bakterií Zavedením nové stavové proměnné x = N K dx dt = dN dt · 1 K = rN 1 − N K · 1 K = rx(1 − x) a nového času τ = rt dostaneme dx dτ = dx dt · dt dτ = rx(1 − x) · 1 r = x(1 − x). Nová stavová proměnná x je bezrozměrná a představuje míru zaplnění Petriho misky, x = 1 je 100% zaplnění Petriho misky do její kapacity. Časová jednotka τ je vůči jednotce t zkrácena nebo prodloužena tak, aby za ni došlo ke zdvojnásobení počtu bakterií v misce. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Počítačová realizace • Maple - vhodný spíše pro teoretické modely • Matlab - vhodný pro maticové zápisy • R - freeware ;-) • Matcont - kontinuační balík pod placený Matlab • XppAut - freeware, vhodný pro parametrickou analýzu • Tabulkové procesory - vhodné pro diskrétní modely ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Zhodnocení modelu: Bylo by jednoduché říct, že zhodnotíme vhodnost modelu nakreslením reálných a simulovaných dat do jednoho grafu a porovnáme je. Není to tak, protože záleží na účelu modelu, krátkodobosti nebo dlouhodobosti predikce, možnostech dobrého odhadu parametrů apod. Žádný model nemůže být realitou, proto zhodnocení modelu nutně v některém okamžiku selže. Je na nás rozhodnout, zda je model už ”dostatečně blízko ”. Daleko jednodušší je porovnávat více modelů mezi sebou. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Proč tak složitě, když známe regresní modely? • regresní model prokládá naměřenými daty křivku a slouží k predikci, JENŽE • je použitelný většinou jen pro krátkodobou predikci • neřekne nic o principu chování systému a vztazích v popisovaném systému • je použitelný jen na konkrétní situaci, výsledky nelze zobecnit • popisuje pouze trend nebo naopak detail, ne obojí • nemůže odhalit, které parametry jsou pro systém podstatné a použitelné např. pro jeho řízení NAROZDÍL OD DETERMINISTICKÉHO MODELU!!! ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Tři základní rady nebojte se LHÁT PODVÁDĚT a KRÁST ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Tři základní rady nebojte se LHÁT PODVÁDĚT a KRÁST ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Tři základní rady nebojte se LHÁT PODVÁDĚT a KRÁST ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Tři základní rady nebojte se LHÁT PODVÁDĚT a KRÁST ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Tři základní rady nebojte se LHÁT PODVÁDĚT a KRÁST To se samozřejmě nemá..., ale zde je to myšleno následovně. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Lhát Každý model obsahuje nekorektní předpoklady. Modely musí být tak zjednodušené, aby množství jejich parametrů a stavových veličin nepřesáhlo dostupná data a možnosti. Zvlášť explikativní modely musí být tak jednoduché, aby bylo vidět co dělají a proč. Reálný svět takový, bohužel a bohudík, není. Proto musí modely ignorovat některá fakta nebo procesy a nahradit je jednoduššími, jistojistě nepravdivými . . . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Podvádět Přesněji, dělejte věci, které budou statistiky znervózňovat, jako například použijte data závislá na jedné proměnné k odhadu parametrů rovnice závislé na mnoha proměnných, používejte znalosti z jiných oborů a používejte intuici. Data jsou pouze jeden z faktorů, které ovlivňují tvorbu modelu, další jsou zkušenost a znalost modelované problematiky. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Krást Nápady si berte odkudkoliv, nezáleží na vědním oboru. Nové vědecké objevy jsou mnohdy výsledkem konvenčních modelů s použitím konvenčních funkčních tvarů v rovnicích - jen v jiném vědeckém oboru. Jestliže již někdo vytvořil rozumný model pro proces, který se objevuje ve vašem modelu, vyzkoušejte ho. Když už někdo věnoval čas a úsilí k odhadu parametrů, použijte ho. Buďte však kritičtí a neváhejte zahodit, co jste si ukradli, pokud to nebude fungovat. Zkuste to spravit a přizpůsobit, třeba to fungovat bude . . . Pokud kradete, citujte odkud. Samozřejmě. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Statické modely a komparativní statika Definice: Statickým modelem rozumíme model nezávislý na běhu času. Popisuje strukturu reálného objektu v rovnovážném stavu. Definice: Komparativní analýzou statického modelu rozumíme analýzu stavových proměnných statického modelu v závisloti na exogenních proměnných a parametrech modelu. Poznámka 1. Rovnice statického modelu nezávisí na čase. Rovnovážný stav je jejich řešením, tedy nalezením stavových proměnných jako funkcí proměnných exogenních a parametrů. Komparativní statiku popisují parciální derivace stavových proměnných podle exogenních proměnných a parametrů. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Třísektorový model uzavřené ekonomiky. Chceme zjistit zda a jak ovlivňují vládní výdaje hrubý národní produkt. Budeme tedy vytvářet teoretický model, jistě realizovatelný. Koncepce: Proměnnými budou jistě vládní výdaje G a hrubý národní produkt (důchod) Y. Vládní výdaje jsou hrazeny z daní T, to bude další proměnná. Na celkovou národní produkci můžeme nahlížet také jako na celkovou sumu peněz za tento produkt, ty jsou rozděleny na tři základní části - investice I, spotřebu C a vládní výdaje G. Uvědomme si jak hrubou agregaci jsme provedli a také jaké předpoklady jsou v pozadí. Tím nejpodstatnějším je, že vše, co je vytvořeno danou ekonomikou, zde je také koupeno. Neexistuje import a export, jde o uzavřenou ekonomiku. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: G Y I α C T γ δ β β Do diagramu doplníme další vztahy. Spotřeba C je závisí zčásti na disponibilním důchodu (Y − T) a zčásti ne - tzv. autonomní výdaje α. Daně T jsou podobně tvořeny daněmi z příjmu Y a jinými typy daní γ. Vzhledem k účelu modelu předpokládáme, že jsou vztahy mezi proměnnými lineární. Vidíme, že proměnné G a I můžeme přesunout mezi exogenní proměnné. Stavovými proměnnými budou Y, C a T. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: Y = I + C + G, C = α + β(Y − T), T = γ + δY, kde α > 0, 0 < β < 1, γ > 0 a 0 < δ < 1 jsou parametry, G a I exogenní proměnné. Rovnice jsou v jednotkách peněz, parametry β a δ jsou bezrozměrné (množství parametrů bychom mohli ještě snížit). Jde o statický model. Vyřešením rovnic dostáváme rovnováhu Y = α − βγ + I + G 1 − β(1 − δ) . Jmenovatel 1 − β(1 − δ) je podle předpokladů kladný, čitatel může být záporný jen v případě vysokého γ, tedy v případě nepřiměřené výše nepříjmových daní dojde ke kolapsu ekonomiky. Účelem modelu bylo zjistit jak ovlivňují vládní výdaje hrubý národní produkt, odtud je zřejmé, že jej zvyšují. Tzv. vládní výdajový multiplikátor je definován ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × jako číslo, které udává o kolik se zvedne GNP, jestliže se zvednou vládní výdaje o jednotku ∂Y ∂G = 1 1 − β(1 − δ) > 0. Relevantními parametry jsou tedy β a δ, které lze odhadnout z naměřených dat. Vyhodnocení: Model je teoretický a těžko porovnatelný s reálnými daty, vzhledem k předpokladu uzavřené ekonomiky. Slouží k pochopení principů a tento účel splnil. Vhodnou revizí by bylo zavedení nových proměnných: exportu a importu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Věta: Nechť x = (x1, . . . , xn) je vektor stavových proměnných, α = (α1, . . . , αm) je vektor exogenních proměnných a parametrů a funkce F = (F1 , . . . , Fn ) : Rn+m → Rn je hladká. Nechť x(α) : Rm → Rn je řešení rovnic modelu F(x, α) = 0, tedy F(x(α), α) = 0. Je-li jacobián F v x nenulový, tj. |DF(x)| = F1 x1 · · · F1 xn ... ... ... Fn x1 · · · Fn xn |x=x = 0, pak je řešení úlohy závislosti rovnovážné stavové proměnné x na některé exogenní proměnné αi řešením soustavy DF(x) · ∂x ∂αi = −Fαi . (1) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Důkaz. Vzhledem k předpokladu hladkosti funkce F, můžeme rovnice modelu F1 (x(α), α) = 0, ... Fn (x(α), α) = 0 derivovat podle proměnné αi, dostáváme tedy F1 x1 ∂x1 ∂αi + · · · + F1 xn ∂xn ∂αi + F1 αi = 0, ... Fn x1 ∂x1 ∂αi + · · · + Fn xn ∂xn ∂αi + Fn αi = 0, což je maticově    F1 x1 · · · F1 xn ... ... ... Fn x1 · · · Fn xn    ·     ∂x1 ∂αi ... ∂xn ∂αi     =    −F1 αi ... −Fn αi    . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Poznámka 2. Parciální derivace vektoru stavových proměnných podle zvoleného parametru ∂x ∂αi = ∂x1 ∂αi , · · · , ∂xn ∂αi T často nemusíme hledat všechny. Vzhledem k tomu, že |DF(x)| = 0, je úloha (1) jednoznačně řešitelná a řešení můžeme hledat pomocí Cramerova pravidla: ∂xj ∂αi = |Dji| |DF(x)| , kde Dji je Jacobiho matice F v x s j-tým sloupcem nahrazeným vektorem −Fαi = (−F1 αi , . . . , −Fn αi )T . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model trhu. Chceme zjistit zda a jak ovlivňuje spotřebitelský důchod cenu výrobku a jeho množství. Koncepce: Proměnnými budou spotřebitelský důchod Y, cena P a množství Q výrobku. Rovnováha trhu se ustanoví, pokud se nabídka S vyrovná poptávce D. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: S Q D P Y = = Do diagramu doplníme další vztahy. Cena P ovlivňuje nabídku S i poptávku D, poptávku ovlivňuje také spotřebitelský důchod Y. Vidíme, že proměnná Y je exogenní. Stavovými proměnnými budou P a Q. Navíc předpokládejme, že platí ∂S ∂P > 0, ∂D ∂P < 0 a ∂D ∂Y > 0. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: F1(P, Q, Y) = S(P) − Q = 0, F2(P, Q, Y) = D(P, Y) − Q = 0. Předpokládejme, že (P, Q) je tržní rovnováha. Derivací rovnic modelu podle Y dostaneme ∂F1(P,Q,Y) ∂Y = ∂S(P) ∂P ∂P ∂Y − ∂Q ∂Y = 0, ∂F2(P,Q,Y) ∂Y = ∂D(P,Y) ∂P ∂P ∂Y − ∂Q ∂Y + ∂D(P,Y) ∂Y = 0. Maticově můžeme rovnice zapsat takto:   ∂S(P) ∂P −1 ∂D(P,Y) ∂P −1   · ∂P ∂Y ∂Q ∂Y = 0 −∂D(P,Y) ∂Y . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Jacobiho matice DF(P, Q, Y) je regulární, protože ∂S(P) ∂P −1 ∂D(P,Y) ∂P −1 = ∂D(P,Y) ∂P − ∂S(P) ∂P < 0. Parciální derivace můžeme tedy jednoznačně vyjádřit pomocí Cramerova pravidla jako ∂P ∂Y = 0 −1 −∂D(P,Y) ∂Y −1 ∂S(P) ∂P −1 ∂D(P,Y) ∂P −1 = ∂D(P,Y) ∂Y ∂S(P) ∂P − ∂D(P,Y) ∂P > 0 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × ∂Q ∂Y = ∂S(P) ∂P 0 ∂D(P,Y) ∂P −∂D(P,Y) ∂Y ∂S(P) ∂P −1 ∂D(P,Y) ∂P −1 = ∂S(P) ∂P ∂D(P,Y) ∂Y ∂S(P) ∂P − ∂D(P,Y) ∂P > 0 Zvýšení příjmů tedy vede ke zvýšení ceny i množství. Relevantními parametry jsou tedy ∂S ∂P, ∂D ∂P a ∂D ∂Y , které lze odhadnout z naměřených dat. Vyhodnocení: Model je teoretický, odpovídá očekávanému výsledku, můžeme jej srovnat s reálnými daty. Simulovat ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Statické modely interakcí a teorie her Pokud model slouží k popisu interakcí subjektů, zvlášť pokud jde o model rozhodovacího procesu v situaci, kdy dochází ke střetu zájmů, je modelem tzv. hra. Teorie her se zabývá analýzou širokého spektra konfliktních i kooperativních rozhodovacích procesů, od aukcí, přes tržní konkurenci, volby, rodinné konflikty až po evoluci a chování zvířat. Teorie her slouží především pro nalezení svým způsobem optimálního řešení konfliktu. Teorie her se zabývá jak statickými, tak dynamickými modely, modely s úplnou informací (deterministickými), tak s neúplnou informací (stochastickými). V této kapitole uvedeme některé statické modely s úplnou informací. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Výběrové řízení. Dvě firmy se zajímají o dva trhy zakázek brněnského magistrátu za 18 a 12 mil. korun. Každá z firem má finanční prostředky buď na velký úplatek jednoho úředníka, nebo na menší úplatky obou úředníků rozhodujících o přidělení zakázek. Předpokládejme, že účinnost úplatků obou firem je stejná a úředníci rozdělují podle těchto pravidel: • Dá-li úplatek jen jedna firma, dostane všechny zakázky trhu. • Dají-li úplatky téhož typu obě firmy, dělí se zakázky na polovinu. • Dá-li jedna firma velký a druhá malý úplatek získá prvně jmenovaná 2/3 zakázek a druhá 1/3 zakázek. Jaké jsou optimální strategie firem? ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × strategie V1 V2 M V1 ( 15, 15) (18, 12) (12, 18) V2 (12, 18) (15, 15) (8, 22) M (18, 12) (22, 8) (15, 15) Strategie obou firem jsou buď velký úplatek V1 prvnímu úředníkovi, velký úplatek V2 druhému úředníkovi nebo dva malé úplatky oběma M (neuplácení ponechme prozatím stranou zájmu). Představme si, že jsme v pozici modré firmy. Pokud by hrála strategii V2 (2. řádek), při jakékoliv volbě strategie červené firmy, získala by méně než při volbě strategie M. Takovouto strategii nazýváme striktně dominovanou jinou strategií, V2 ≺ M. Stejně tak V1 ≺ M. Striktně dominované strategie modrá firma nebude hrát, stejně se zachová i červená firma. Obě zvolí strategii dvou malých úplatků. Toto řešení má tu vlastnost, že při jednostranném odchýlení od této strategie si ani jedna firma nepolepší, říkáme mu rovnovážné řešení. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Definice: Hra v normálním tvaru pro n hráčů je tvořena prostory strategií jednotlivých hráčů S1, . . . , Sn a jejich výplatními funkcemi u1, . . . un , kde každé ui zobrazuje S1 × · · · × Sn do R. Označením ui(si, s−i) budeme rozumět ui(s1, . . . , sn), kde sj ∈ Sj. Definice: Nechť s′ i, s′′ i ∈ Si jsou dvě možné strategie i-tého hráče. Řekneme, že strategie s′ i je striktně dominovaná strategií s′′ i , s′ i ≺ s′′ i , jestliže pro každou kombinaci strategií ostatních hráčů je výplata i-tého hráče při strategii s′ i menší než při strategii s′′ i , tj. ui(s′ i, s−i) < ui(s′′ i , s−i), pro libovolné strategie protihráčů s−i. 1. příklad: Ukažte, že strategie nedávat úplatek nebo dát jen jeden malý úplatek je striktně dominovaná strategií uplatit oba. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Definice: Strategie s∗ 1, . . . , s∗ n tvoří Nashovu rovnováhu, jestliže pro každého hráče je s∗ i nejlepší odpovědí na strategie s∗ −i ostatních, tedy ui(s∗ i , s∗ −i) ≥ ui(si, s∗ −i) pro libovolné si ∈ Si. Jinak řečeno, s∗ i je řešením extremální úlohy max si∈Si ui(si, s∗ −i). Pokud při eliminaci striktně dominovaných strategií zůstane jediná kombinace strategií, je jedinou Nashovou rovnováhou. Eliminací striktně dominovaných strategií obecně zmenšíme hru a pokud existuje Nashova rovnováha, zůstává mezi zbylými strategiemi menší hry. Obecně Nashova rovnováha nemusí existovat (v takto zavedených, tzv. ryzích strategiích). Navíc pokud existuje nemusí být pareto-optimální, tj. může existovat strategie s lepší výplatou pro daného hráče přičemž ostatní si nepohorší. Tato strategie ale není rovnovážná, protože vychýlení z této strategie by bylo pro některého hráče výhodnější. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Cournotův model trhu. Chceme nalézt optimální množství výrobků, jež budou ochotny na trh dodávat firmy. Koncepce: Exogenními proměnnými budou poptávané množství M a mezní náklady c na výrobu jednoho výrobku, endogenní proměnné jsou množství qi výrobků od jednotlivých firem. Model bude statický firmy se v daném okamžiku rozhodnou a nezávisle na sobě volí optimální strategii. Volme tyto zjednodušující předpoklady: • poptávková funkce je lineární tvaru P(Q) = M − Q , kde Q je celkové množství dodávané na trh pro Q ≥ M je P(Q) = 0 • postavení firem je rovnocenné a jejich produkt je homogenní, tj. Q = q1 + · · · + qn. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × • mezní náklady c jsou konstantní, tj. nákladová funkce Ci(qi) = cqi, c < M • výstup je libovolně dělitelný, prostor strategií tak můžeme označit Si = [0, ∞) Rovnice: Výplatní funkcí je zisková funkce firem: πi(qi, q−i) = qi P(Q) − c = qi M − (q1 + · · · + qn) − c Abychom našli Nashovu rovnováhu tohoto problému, musí každá firma řešit optimalizační problém max 0≤qi<∞ πi(qi, q∗ −i) = max 0≤qi<∞ qi M − (qi + ∑q∗ −i) − c Úkolem je tedy najít maximum kvadratické funkce v proměnné qi. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Podmínky prvního řádu tedy jsou: 2q1 + q2 + · · · + qn = M − c q1 + 2q2 + · · · + qn = M − c ... = ... q1 + q2 + · · · + 2qn = M − c Řešením jsou hodnoty q∗ 1 = q∗ 2 = . . . = q∗ n = M − c n + 1 , které jsou vzhledem ke konkávnosti funkcí πi(qi, q∗ −i) maximem. Cena odpovídá poptávkové funkci P∗ = M − Q = M − n M − c n + 1 = 1 n + 1 M + n n + 1 c a zisk firmy je π∗ = q∗ P∗ − c = M − c n + 1 1 n + 1 M + n n + 1 c − c = M − c n + 1 2 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Q∗ P∗ π∗ monopol 1 2 (M − c) 1 2 M + 1 2 c 1 4 (M − c)2 duopol 2 3 (M − c) 1 3 M + 2 3 c 1 9 (M − c)2 oligopol n n + 1 (M − c) 1 n + 1 M + n n + 1 c M − c n + 1 2 dok. konkurence M − c c 0 V případě monopolu, je na trhu jediná firma nabízející množství q∗ = M − c 2 v případě duopolu nabízí dvě firmy q∗ = M − c 3 . Limitním případem je pro n → ∞ dokonalá konkurence, jež na trh dodává celkové množství lim n→∞ n n + 1 (M − c) = M − c, při němž firmy dosahují nulového zisku. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Monopolista dodává na trh menší množství výrobků než duopolisté a prodává ho za vyšší cenu. Pokud porovnáme celkový zisk monopolisty a obou duopolistů, je vidět, že by pro duopolisty bylo výhodnější vyrábět polovinu produkce monopolisty. Tato strategie ale není rovnovážná, každá firma by v takové situaci musela odolávat pokušení nevyrábět více, z čehož by získala výhodu, pokud by druhá firma udržovala produkci na dané hladině. V takové situaci mohou firmy uzavřít dohodu, problém však spočívá v tom, že vzhledem k antimonopolním opatřením jsou takové dohody protizákonné a tajná dohoda je legálními prostředky nevymahatelná. Pro firmy je ale natolik výhodné vytvářet velké kartely a chovat se jako monopolista, že dokonalá soutěž zůstává jen v myšlenkách idealistů. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 2. příklad: V Cournotově modelu duopolu nakreslete tzv. reakční křivky obou firem, tj. funkce qi = R(q−i), které udávají nejlepší odpovědi (best response functions). Ukažte zde Nashovu rovnováhu. 3. příklad: Ukažte, že v případě monopolu jsou všechny ostatní strategie striktně dominované strategií q∗ = M − c 2 . Nápověda: ostatní strategie jsou q = q∗ + x, kde x = 0. 4. příklad: Najděte Nashovu rovnováhu v Bertrandově modelu duopolu, kde se dva výrobci rozhodují o optimální ceně za poptávané množství qi(pi, p−i) = a − pi + bp−i, které závisí na ceně obou výrobků. Přitom b ∈ 0, 2) je tzv. elasticita nebo míra substituce. Řešení v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × V ryzích strategiích ne vždy existuje Nashova rovnováha. Proto zavádíme pojem smíšené strategie, která udává pravděpodobnost volby ryzí strategie. Smíšená strategie je tedy pro každého hráče vektor, jehož k-tá složka udává pravděpodobnost, s níž hráč volí k-tou strategii ze svého prostoru strategií. Zde už existenci rovnováhy zaručuje Nashova věta. Definice: Uvažujme konečnou hru n hráčů v normálním tvaru, kde počet prvků prostoru strategií Si libovolného hrače i označíme symbolem mi. Smíšenou strategií hráče i se rozumí vektor pravděpodobností xi = (xi 1, . . . , xi mi )T , kde xi j ≥ 0 a ∑ j xi j = 1. Výplatní funkcí pro smíšenou strategii xi i-tého hráče je pak vážený průměr výplatních funkcí ui(s j i, s−i) všech možných strategií s1 i , . . . s mi i s vahami danými pravděpodobnostmi hrát danou strategii. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Věta (Nashova): Konečná hra n hráčů má v prostoru smíšených strategií alespoň jednu Nashovu rovnováhu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Dynamické modely Definice: Dynamickým modelem rozumíme model závislý na běhu času. Popisuje chování reálného objektu v průběhu času. Dynamický model je popsán dynamickým systémem (rovnicí, soustavou rovnic, formulí). Definice: Dynamickým systémem rozumíme trojici {T, X, ϕt }, kde T je číselná množina (čas), X je metrický prostor, který nazýváme fázovým prostorem, a ϕt je parametrický systém evolučních operátorů s parametrem t ∈ T definovaných jako zobrazení ϕt : X → X, které zobrazuje počáteční stav x0 ∈ X na nějaký stav xt = ϕt x0 ∈ X. Poznámka 3. V případě, že T = Z mluvíme o diskrétním dynamickém systému, je-li T = R mluvíme o spojitém dynamickém systému. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Definice: Deterministickým dynamickým systémem rozumíme systém {T, X, ϕt } splňující podmínku ϕ0 = id, kde id je identita na X, tj. ∀x ∈ X : id x = x. Tato vlastnost říká, že systém spontánně nemění svůj stav. Definice: Autonomním dynamickým systémem rozumíme deterministický systém {T, X, ϕt } splňující podmínku ϕt+s = ϕt ◦ ϕs , tj. ∀x ∈ X : ϕt+s x = ϕt (ϕs x), pokud jsou definovány obě strany rovnosti. Tato vlastnost říká, že se „zákony evoluce“ nemění během času. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Definice: Trajektorie s počátečním bodem x0 ∈ X je uspořádaná podmnožina fázového prostoru X {x ∈ X : x = ϕt x0, ∀t ∈ T, pro které je ϕt x0 definováno} V případě spojitého systému jde o orientované křivky v X, v případě diskrétního systému jsou to posloupnosti bodů v X. Fázovým portrétem dynamického systému rozumíme rozmístění trajektorií ve fázovém prostoru X. Definice: Bod x0 ∈ X nazýváme rovnovážným bodem (nebo též singulárním, stacionárním, pevným bodem) dynamického systému, jestliže pro všechna t ∈ T platí ϕt x0 = x0. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Definice: Cyklem rozumíme periodickou trajektorii L0, která není rovnovážným bodem, splňující ∀x0 ∈ L0 ϕt+T0 x0 = ϕt x0, pro nějaké T0 > 0, ∀t ∈ T. Nejmenší takové T0 nazýváme periodou cyklu L0. Poznámka 4. V systému s cyklem vznikají periodické oscilace. Cyklus spojitého systému je uzavřená křivka v X. Limitním cyklem rozumíme cyklus, v jehož okolí nejsou jiné cykly. Definice: Invariantní množinou S rozumíme podmnožinu X splňující x0 ∈ S ⇒ ϕt x0 ∈ S ∀t ∈ T. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Poznámka 5. Rovnovážný bod i cyklus jsou invariantní množiny. Definice: Invariantní množina S se nazývá stabilní, jestliže • ∀U ⊃ S libovolně malé okolí invariantní množiny existuje okolí V ⊃ S takové, že ∀x ∈ V a ∀t > 0 platí ϕt x ∈ U (tento typ stability nazýváme ljapunovskou stabilitou), • existuje okolí U0 ⊃ S takové, že ϕt x → S pro x ∈ U0 a t → ∞ (tento typ stability nazýváme asymptotickou stabilitou). V opačném případě je S nestabilní. Poznámka 6. Existují další typy stability, my se budeme většinou setkávat s rovnovážnými body a cykly, které jsou jak ljapunovsky, tak asymptoticky stabilní. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnovážná dynamika Dynamické modely oproti statickým modelů zachycují vývoj stavových veličin v čase. Až do druhé poloviny minulého století se v aplikovaných vědách objevovaly většinou dynamické modely, které směřovaly k rovnovážnému stavu. Implicitně se tedy předpokládalo, že dynamický systém z libovolné relevantní počáteční hodnoty směřuje k rovnováze, což je přesně pojem stability, dokonce asymptotické. Uveďme jako příklad neviditelnou ruku trhu, která má za každých ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × okolností přivést ekonomický systém k makroekonomické rovnováze. V biochemii uveďme například Michaelisův-Mentenové model enzymatické reakce, který si později podrobně rozebereme. Celá klasická termodynamika předpokládá postupné směřování systému k rovnováze (vyrovnání teplot a postupné dosažení maximální entropie). Takováto stabilní dynamická rovnováha odpovídá právě statické rovnováze, kterou jsme studovali v předchozích modelech. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Základním principem takovýchto modelů je následující úvaha. Čím více se systém odchýlí od své rovnováhy, tím větší má tendenci k ní směřovat. Tato úvaha je v mnohých případech velmi racionální a aplikovatelná na velké množství situací. Tato úvaha v sobě ale implictně zahrnuje existenci dynamické rovnováhy a její asymptotickou stabilitu. Takovýto předpoklad nutně vede k rovnovážné dynamice. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × • Uveďme jako základní příklad Newtonův zákon ochlazování, kdy teplota tělesa se mění tím rychleji, čím větší je rozdíl teplot tělesa a jeho okolí. • Stejně tak bychom ale mohli použít lineární makroekonomický model nabídky a poptávky, kdy růst nabídky je tím větší, čím větší je převis poptávky nad nabídkou atd. • Stáda antilop migrují společně a pokud se některá dostane mimo stádo, má tím větší tendenci se k němu připojit, čím dál od něj je. • Dokonce i chování lidí je možné tímto způsobem modelovat. Většina lidí má tendenci nevybočovat z davu a své chování měnit tím více, čím větší je jeho odlišnost od běžné normy. Můžeme tím vysvětlit např. to, že i ateisté zmlknou v kostele nebo že si i přísný abstinent dá na Silvestra skleničku sektu, byť ji nevypije. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Newtonův model ochlazování. Koncepce: Představme si kelímek kávy právě vytažené z automatu (o teplotě T0) a postavené do místnosti s teplotou T∗ . Stavová proměnná bude teplota kávy T, parametrem bude k ∈ R, které bude záviset ostatních fyzikálních veličinách (měrná tepelná kapacita kávy, tvar a plast kelímku apod.). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: Rovnice: dT(t) dt = k(T∗ − T(t)). (2) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 5. příklad: Vyřešte rovnici s počáteční podmínkou T(0) = T0. Odhadněte k pro konkrétní hrnek kafe. 6. příklad: Najděte rovnovážný bod rovnice a určete jeho stabilitu. Odhadněte, za jak dlouho káva ”vystydne”. Vyhodnocení: Teoretické výsledky srovnejte s měřením. Pokud neodpovídají, vysvětlete a navrhněte revizi. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Věta: Uvažujme dostatečně hladké zobrazení f : Rm → Rm a rovnici x′ = f (x). (3) Rovnovážný bod x∗ spojitého systému (3) splňuje f (x∗ ) = 0. Uvažujme nejprve případ m = 1 a Taylorův rozvoj f v rovnovážném bodě x∗ . Pro x ≈ x∗ platí f (x) ≈ f (x∗ ) + D f (x∗ )(x − x∗ ) + · · · = D f (x∗ )(x − x∗ ) + . . . . V dostatečně blízkém okolí x∗ tedy platí x′ ≈ D f (x∗ )(x − x∗ ). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Věta: Mějme rovnici (3) pro m = 1 a f hladkou v okolí rovnovážného bodu x∗ . Jestliže D f (x∗ ) < 0, pak je rovnovážný bod x∗ stabilní (atraktor). V opačném případě, když D f (x∗ ) > 0, je x∗ nestabilní (repeler). x∗ x∗ Df(x∗ ) < 0 Df(x∗ ) > 0 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Základní spojité modely růstu Spojitý exponenciální růst. Uvažujme populaci, kterou můžeme modelovat spojitě - např. množství sinic na přehradě budeme měřit v g/m3 , nebudeme je počítat. Stejně tak budeme spojitý přístup používat u populace, která nemá daná období rozmnožování (jako má mnoho druhů zvířat narozdíl od člověka). Zajímá nás, jak se bude populace vyvíjet v čase. Označíme-li x(t) velikost populace v čase t, b okamžitou míru reprodukce a d míru vymírání (zde předpokládáme, že jsou míry b a d konstantní), pak můžeme populaci popsat diferenciální rovnicí x′ = bx − dx = rx, kde r je konstantní míra růstu populace a x′ představuje okamžitou ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × změnu velikosti populace. Připomeňme zde definici derivace: x(t + ∆t) − x(t) ∆t = rx(t), pro ∆t → 0. Řešením je exponenciální funkce x(t) = x0ert . Pokud je r < 0, tj. d > b, populace vymře (rovnovážný stav x(t) ≡ 0 je stabilní), pokud je r > 0, tj. d < b, populace bude růst nade všechny meze (rovnovážný stav x(t) ≡ 0 je nestabilní). 7. příklad: Kdy je třeba vyhlásit zákaz koupání v přehradě, jestliže jsme minulé 4 dny v 8:00 ráno naměřili v odběrné nádobě hodnoty 2, 3, 5 a 7 mikrogramů. Hranice toxicity je 30 µg. Simulovat v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Model lze použít v případech, kdy nám stačí krátkodobá předpověď, nebo je-li dynamika populace vzhledem k jiné modelované proměnné daleko pomalejší. Takovým příkladem může být například dynamika trhu práce a kapitálu v ekonomii, kdy dynamiku trhu práce můžeme popsat rovnicí s konstantní mírou růstu práce (odpovídající míře růstu populace). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Spojitý logistický růst. Uvažujme nyní tuto modifikaci předchozího modelu: x′ = r(x)x, kde míra růstu populace r(x) závisí na velikosti populace. Volbou r(x) dostáváme následující rovnice populačního růstu: r(t, x) = r0 1 − x K logistická Verhulstova rovnice, r(t, x) = r0 1 − x K β , β > 0 Richardsova rovnice, r(t, x) = r0 1 − x K 1 + c x K , c > 0 Smithova rovnice, r(t, x) = r0 ln K x Gompertzova rovnice atd. Všechny uvedené rovnice jsou autonomní, r nezávisí na čase explicitně, ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × pouze v závislosti na velikosti populace. K > 0 je tzv. kapacita prostředí. 8. příklad: S pomocí Maplu pro uvedené rovnice nakreslete řešení počáteční úlohy x0 = 3, pro r0 = 2, K = 100 a vhodně volené případné další parametry, nakreslete také funkce r(t, x). Najděte obecná řešení rovnic a zkoumejte jejich tvar. Najděte inflexní body řešení a vysvětlete, co znamenají. Řešení v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 9. příklad: S pomocí vhodného programu (nebo bez něj) najděte rovnovážné body výše uvedených rovnic a vyšetřete jejich stabilitu. x′ = f (x) := r0 1 − x K x, logistická Verhulstova rovnice Předpokládejme, že r0 > 0 (typický případ). s.b.: x = 0, x = K, D f (x) = r0 1 − x K − r0 K x D f (0) = r0 > 0, x = 0 je nestabilní s.b. D f (K) = −r0 < 0, x = K je stabilní s.b. Řešení v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 9. příklad: S pomocí vhodného programu (nebo bez něj) najděte rovnovážné body výše uvedených rovnic a vyšetřete jejich stabilitu. x′ = f (x) := r0 1 − x K x, logistická Verhulstova rovnice Předpokládejme, že r0 > 0 (typický případ). s.b.: x = 0, x = K, D f (x) = r0 1 − x K − r0 K x D f (0) = r0 > 0, x = 0 je nestabilní s.b. D f (K) = −r0 < 0, x = K je stabilní s.b. Řešení v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 9. příklad: S pomocí vhodného programu (nebo bez něj) najděte rovnovážné body výše uvedených rovnic a vyšetřete jejich stabilitu. x′ = f (x) := r0 1 − x K x, logistická Verhulstova rovnice Předpokládejme, že r0 > 0 (typický případ). s.b.: x = 0, x = K, D f (x) = r0 1 − x K − r0 K x D f (0) = r0 > 0, x = 0 je nestabilní s.b. D f (K) = −r0 < 0, x = K je stabilní s.b. Řešení v Maplu Pro rovnovážný bod x∗ platí x′ = f (x) = 0. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 9. příklad: S pomocí vhodného programu (nebo bez něj) najděte rovnovážné body výše uvedených rovnic a vyšetřete jejich stabilitu. x′ = f (x) := r0 1 − x K x, logistická Verhulstova rovnice Předpokládejme, že r0 > 0 (typický případ). s.b.: x = 0, x = K, D f (x) = r0 1 − x K − r0 K x D f (0) = r0 > 0, x = 0 je nestabilní s.b. D f (K) = −r0 < 0, x = K je stabilní s.b. Řešení v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 9. příklad: S pomocí vhodného programu (nebo bez něj) najděte rovnovážné body výše uvedených rovnic a vyšetřete jejich stabilitu. x′ = f (x) := r0 1 − x K x, logistická Verhulstova rovnice Předpokládejme, že r0 > 0 (typický případ). s.b.: x = 0, x = K, D f (x) = r0 1 − x K − r0 K x D f (0) = r0 > 0, x = 0 je nestabilní s.b. D f (K) = −r0 < 0, x = K je stabilní s.b. Řešení v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 9. příklad: S pomocí vhodného programu (nebo bez něj) najděte rovnovážné body výše uvedených rovnic a vyšetřete jejich stabilitu. x′ = f (x) := r0 1 − x K x, logistická Verhulstova rovnice Předpokládejme, že r0 > 0 (typický případ). s.b.: x = 0, x = K, D f (x) = r0 1 − x K − r0 K x D f (0) = r0 > 0, x = 0 je nestabilní s.b. D f (K) = −r0 < 0, x = K je stabilní s.b. Řešení v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 9. příklad: S pomocí vhodného programu (nebo bez něj) najděte rovnovážné body výše uvedených rovnic a vyšetřete jejich stabilitu. x′ = f (x) := r0 1 − x K x, logistická Verhulstova rovnice Předpokládejme, že r0 > 0 (typický případ). s.b.: x = 0, x = K, D f (x) = r0 1 − x K − r0 K x D f (0) = r0 > 0, x = 0 je nestabilní s.b. D f (K) = −r0 < 0, x = K je stabilní s.b. Řešení v Maplu Podobně pro další rovnice. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Harrodův-Domarův model ekonomického růstu Modelujme růst hrubého domácího produktu. Koncepce: Uvažujme uzavřenou ekonomiku a předpokládejme, že jsou splněny následující podmínky. Kapitál K vzniká investicemi I, přitom dochází k jeho amortizaci. Spotřeba a úspory S jsou pevným podílem produktu Y, zbytek produktu investujeme do tvorby kapitálu. Relativní přírůstek kapitálu se projevuje relativním přírůstkem produkce. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: Y S amortizace K I 1 − s s δ Je zřejmé, že amortizaci, spotřebu a úspory lze odvodit z produktu, máme tedy pouze tři stavové proměnné: produkt Y(t) > 0, kapitál K(t) > 0 a investice I(t) > 0. Míra úspor a spotřeby je označena s (mezní sklon k úsporám a spotřebě), míra amortizace δ. Zřejmě s ∈ (0, 1) a δ ∈ (0, 1). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: K′ = I − δK, I = (1 − s)Y, K′ K = Y′ Y . Všimněme si nyní, že platí 0 = K′ K − Y′ Y = K′Y−Y′K Y2 Y K = K Y ′ Y K . Odtud K Y ′ = 0. Existuje tedy konstanta r ∈ R, taková že K Y = r. Toto číslo můžeme interpretovat jako kapitálovou náročnost jednotky produkce. 10. příklad: Odvoďte diferenciální rovnici pro růst produktu a vyřešte ji. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Závěr analýzy modelu nyní můžeme přeformulovat: je-li r < 1−s δ pak produkce roste, je-li r = 1−s δ pak produkce stagnuje, je-li r > 1−s δ pak produkce klesá. To odpovídá zkušenosti: je-li kapitálová náročnost jednotky produkce příliš velká, pak produkce nemůže růst. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diskrétní exponenciální růst - malthusovský model. Uvažujme populaci, která se rozmnožuje a vymírá v pevně daných intervalech. Může jít o jakoukoliv populaci - ryb, rostlin nebo peněz. Může jít také o populaci, která je v pevných časových intervalech kontrolována a jiné informace o ní nemáme. Zajímá nás, jak se bude populace vyvíjet v čase. Označíme-li xn velikost populace v čase n, b míru reprodukce a d míru vymírání (zde předpokládáme, že jsou míry b a d konstantní), pak můžeme populaci v následujícím čase n + 1 popsat diferenční rovnicí xn+1 − xn = bxn − dxn = rxn, kde r je konstantní míra růstu populace, neboli xn+1 = µxn, kde µ = 1 + r = 1 + b − d. Simulovat v Matlabu exponencialnirust.m Řešením je geometrická posloupnost xn = x0µn . Pokud je µ < 1, tj. d > b, populace vymře (rovnovážný stav xn ≡ 0 je stabilní), pokud je ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × µ > 1, tj. d < b, populace bude růst nade všechny meze (rovnovážný stav xn ≡ 0 je nestabilní). Věta: Uvažujme dostatečně hladké zobrazení f : Rm → Rm a rovnici xn+1 = f (xn). (4) Pevný bod x∗ diskrétního systému (4) splňuje f (x∗ ) = x∗ . Uvažujme nejprve případ m = 1 a Taylorův rozvoj f v pevném bodě x∗ . Pro x ≈ x∗ platí f (x) ≈ f (x∗ ) + D f (x∗ )(x − x∗ ) + · · · = x∗ + D f (x∗ )(x − x∗ ) + . . . . V dostatečně blízkém okolí x∗ tedy platí xn+1 − x∗ ≈ D f (x∗ )(xn − x∗ ). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Věta: Mějme zobrazení(4) pro m = 1, hladké v okolí pevného bodu x∗ . Jestliže |D f (x∗ )| < 1, pak |xn+1 − x∗ | < |xn − x∗ |, a pevný bod x∗ je stabilní (atraktor). V opačném případě, když |D f (x∗ )| > 1, je x∗ nestabilní (repeler). Pavučinový diagram: Vhodným zobrazením dynamiky zobrazení (4) pro m = 1 je následující graf: f(x) x xn xn+1 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diskrétní logistický růst - Verhulstův model. Uvažujme nyní takovou modifikaci předchozího modelu, že míra růstu r bude lineárně klesat v závislosti na velikosti populace yn. Pokud dosáhne populace určité velikosti K, kterou nazýváme kapacita prostředí, bude míra růstu nulová, pokud tuto kapacitu překročí, bude velikost populace klesat, tj. yn+1 − yn = r 1 − yn K yn. Pokud je r = 0 (triviální případ), můžeme provést transformaci yn = 1 + r r Kxn, kterou zmenšíme počet parametrů: xn+1 = µxn(1 − xn), kde µ = 1 + r. (5) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 11. příklad: Tato rovnice má dva pevné body. Najděte je, určete pro ně podmínky stability za předpokladu, že r ∈ (0, 2). Co se děje, pro o něco vyšší hodnoty r? Vzniká stabilní cyklus periody 2. Připomeňme, že cyklus periody 2 je uspořádaná dvojice [x1, x2], kde x1 = f (x2) = f ( f (x1)) = f (2) (x1), x1 je tedy pevným bodem zobrazení f (2) (x) = µ2 x(1 − x)(1 − µx(1 − x)). 12. příklad: Najděte všechna řešení rovnice f (2) (x) = x pro r = 2.1 a ukažte, že cyklus periody 2 je stabilní. Rada: vyřaďte ta řešení, která jsou zároveň pevným bodem f (x) (proč?), spočtěte v nich D f (2) . Výpočet v Maplu V programu XppAut spusťte soubor cobweb.ode Simulovat v Matlabu logistickyrust.m ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Pokud zakreslíme závislost pevných bodů na parametru µ, dostaneme tzv. bifurkační diagram. Postupné zdvojování periody přechází v deterministický chaos. V programu XppAut spusťte soubor logbif.ode Co je to chaos? Slovo chaos je řeckého původu a znamená nepředvídatelnost. Deterministický chaos je neperiodické deterministické chování, které je • velice citlivé na počáteční podmínky, • topologicky transitivní - což znamená, že libovolný interval transformuje na libovolný další interval • má husté trajektorie DETERMINISTICKÝ NEZNAMENÁ PŘEDVÍDATELNÝ!!! ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Nerovnovážná dynamika Jak je vidět, předpoklad samovolného asymptotického směřování systému k jeho rovnováze lze docela jednoduše narušit. Vznik chaotického nepředvídatelného chování trajektorie diskrétní logistické rovnice je toho důkazem. Od 70. let 20. století začíná získávat nerovnovážná dynamika ve většině aplikovaných věd své místo a nelineární dynamika otvírá cestu pro propojení deterministického a stochastického modelování. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Řízení (kontrola) chaosu metodou OGY V roce 1990 Ott, Grebogi a Yorke uvedli praktickou metodu (úspěšnou i v aplikacích) stabilizace nestabilních chaotických cyklů. Metoda je založena na faktu, že chaotický atraktor obsahuje nekonečné husté množství nestabilních cyklů. Ty jsou stabilizovány malými perturbacemi kontrolního parametru a. Uvažujme zobrazení xn+1 = f (xn, a), (6) kde a je dostupný parametr, který můžeme změnit v nějakém okolí své ”nominální”hodnoty a0. Označme x∗ (a) nestabilní pevný bod zobrazení (6). V malém okolí a0 můžeme aproximovat xn+1 − x∗ (a0) = D f (x∗ (a0), a0)(xn − x∗ (a0)) + c(a − a0), (7) kde c = ∂ f ∂a (x∗ (a0), a0). Vzhledem k transitivnosti a hustotě chaotické trajektorie musí v nějakém malém okolí x∗ (a0) pro nějaké xn platit a − a0 = −k(xn − x∗ (a0)). (8) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Substitucí (8) do (7) dostaneme xn+1 − x∗ (a0) = (D f (x∗ (a0), a0) − ck)(xn − x∗ (a0)). Volbou k můžeme dosáhnout stability regulovaného pevného bodu, tj. najdeme k tak, aby |D f (x∗ (a0), a0) − ck| < 1. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Řízení (kontrola) chaosu v logistickém zobrazení Uvažujme logistickou rovnici (5), ve které ovlivňujeme dynamiku neustálými pulzy xi = kxi po p iteracích. Definujme zobrazení F(x) = k f (p) (x). Pevný bod x∗ regulovaného zobrazení F(x) tedy bude splňovat k f (p) (x∗ ) = x∗ a bude stabilní, pokud |kD f (p) (x∗ )| < 1. Označíme-li Cp (x) = x f (p)(x) D f (p) (x), dostáváme podmínku pro oblast hodnot, pro které jsme schopni chaos změnit ve stabilní dynamiku: |Cp (x)| < 1. Výpočet Cp v Maplu Simulace v Matlabu chaoscontrol.m ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × CESTA DO VÍCE STAVOVÝCH PROMĚNNÝCH... ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Strukturovaný spojitý dynamický model Strukturované modely se používají v případě, že je potřeba rozlišovat složky stavové proměnné podle nějakého kritéria, které ovlivňuje dynamiku. Typickým příkladem jsou epidemiologické modely, kdy v populaci rozlišujeme jedince v různých stádiích nemoci. Účelem modelu je porozumět průběhu epidemie a předpovědět, kdy epidemie odezní. Modely použitelné např. na reálné chřipkové epidemie jsou samozřejmě komplikovanější, než v této přednášce uvedené základní epidemiologické modely, princip je však stejný. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model SI. Chceme modelovat epidemii infekční nemoci, kterou neumíme léčit, která však není smrtelná, např. herpes labialis, opar rtu. Koncepce: Stavovou proměnnou budou infikovaní jedinci I a náchylní jedinci S. Předpokládáme nulovou úmrtnost způsobenou nemocí a také rovnováhu mezi počtem nově narozených a přirozeně zemřelých jedinců. Toto hrubé zjednodušení můžeme použít, pokud rychlost šíření infekční nemoci je podstatně větší než růst populace. Parametrem bude samozřejmě rychlost šíření infekce β > 0. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: S I β V čase t = 0 existuje S0 > 0 náchylných jedinců a I0 > 0 nakažlivých. Můžeme předpokládat, že počet nově infikovaných je přímo úměrný počtu náchylných a nakažlivých jedinců. Koeficient β bude závislý na četnosti kontaktů v populaci a pravděpodobnosti nákazy při kontaktu náchylného a nakažlivého jedince. Rovnice: Model je popsán následujícím systémem diferenciálních rovnic: S′ = −βSI, I′ = βSI. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Řešením počáteční úlohy S(0) = S0, I(0) = I0, S(t) + I(t) ≡ N. je funkce I(t) = N 1 + N I0 − 1 e−βNt Grafem je logistická křivka, která má inflexní bod   ln N I0 − 1 βN , N 2   . Z hlediska dynamiky je zajímavý graf funkce I′ (t) = βN2 N I0 − 1 eβNt N I0 − 1 + eβNt 2 který ukazuje přírůstky infikovaných. Výpočet a simulace v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Je vidět, že celá populace se nakonec nakazí, což jsme očekávali. U oparu např. je v dospělosti nakaženo asi 90% populace. 10. příklad: Najděte rovnovážné body rovnice infikovaných jedinců a vyšetřete jejich stabilitu. I′ = f (I) := β(N − I)I, rovnovážné body: I = 0, I = N, D f (I) = β(N − 2I), D f (0) = βN > 0, I = 0 je nestabilní s.b. D f (N) = −βN < 0, I = N je stabilní s.b. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Je vidět, že celá populace se nakonec nakazí, což jsme očekávali. U oparu např. je v dospělosti nakaženo asi 90% populace. 10. příklad: Najděte rovnovážné body rovnice infikovaných jedinců a vyšetřete jejich stabilitu. I′ = f (I) := β(N − I)I, rovnovážné body: I = 0, I = N, D f (I) = β(N − 2I), D f (0) = βN > 0, I = 0 je nestabilní s.b. D f (N) = −βN < 0, I = N je stabilní s.b. V každém okamžiku platí S(t) = N − I(t). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Je vidět, že celá populace se nakonec nakazí, což jsme očekávali. U oparu např. je v dospělosti nakaženo asi 90% populace. 10. příklad: Najděte rovnovážné body rovnice infikovaných jedinců a vyšetřete jejich stabilitu. I′ = f (I) := β(N − I)I, rovnovážné body: I = 0, I = N, D f (I) = β(N − 2I), D f (0) = βN > 0, I = 0 je nestabilní s.b. D f (N) = −βN < 0, I = N je stabilní s.b. Pro rovnovážný bod platí I′ = f (I) = 0. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Je vidět, že celá populace se nakonec nakazí, což jsme očekávali. U oparu např. je v dospělosti nakaženo asi 90% populace. 10. příklad: Najděte rovnovážné body rovnice infikovaných jedinců a vyšetřete jejich stabilitu. I′ = f (I) := β(N − I)I, rovnovážné body: I = 0, I = N, D f (I) = β(N − 2I), D f (0) = βN > 0, I = 0 je nestabilní s.b. D f (N) = −βN < 0, I = N je stabilní s.b. Jacobiho ”matice”, v jednorozměrném případě derivace pravé strany ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Je vidět, že celá populace se nakonec nakazí, což jsme očekávali. U oparu např. je v dospělosti nakaženo asi 90% populace. 10. příklad: Najděte rovnovážné body rovnice infikovaných jedinců a vyšetřete jejich stabilitu. I′ = f (I) := β(N − I)I, rovnovážné body: I = 0, I = N, D f (I) = β(N − 2I), D f (0) = βN > 0, I = 0 je nestabilní s.b. D f (N) = −βN < 0, I = N je stabilní s.b. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Je vidět, že celá populace se nakonec nakazí, což jsme očekávali. U oparu např. je v dospělosti nakaženo asi 90% populace. 10. příklad: Najděte rovnovážné body rovnice infikovaných jedinců a vyšetřete jejich stabilitu. I′ = f (I) := β(N − I)I, rovnovážné body: I = 0, I = N, D f (I) = β(N − 2I), D f (0) = βN > 0, I = 0 je nestabilní s.b. D f (N) = −βN < 0, I = N je stabilní s.b. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Poznámka 7. Je evidentní, že pro použití modelu bude nejpodstatnější odhad parametru β. Zkusme najít průměrný počet nakažených za jednotku času. Aby se někdo nakazil, musí se setkat infikovaný jedinec s náchylným a musí dojít k nákaze. Jaká je pravděpodobnost, že vybereme z N lidí jednoho infikovaného a jednoho náchylného? p = I N S N − 1 + S N I N − 1 ≈ 2 N2 SI. Tuto aproximaci můžeme provést ve velké skupině lidí, kde N2 >> N, jinak je třeba použít prvně uvedený vztah. Pokud γ > 0 označíme průměrný počet interakcí za jednotku času a c průměrný počet nakažení při SI interakci, tj. 0 < c ≤ 1, je počet nově nakažených za jednotku času I(t + ∆t) − I(t) ∆t = 2cγ N2 SI. Provedením limitního přechodu t → 0 dostáváme β = 2cγ N2 . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Později se podíváme na složitější epidemiologické modely, např. model SIR a SIRS. K tomu ale budeme potřebovat něco málo další teorie, protože vstupujeme do fázového prostoru o více než jednom rozměru. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Spojitá a diskrétní dynamika v Rm . Lineární algebra - připomenutí Pro vlastní číslo (vlastní hodnotu) matice A ∈ Rm×m příslušné vlastnímu vektoru v ∈ Rm platí Av = λv, tj. vlastní čísla hledáme jako kořeny charakteristického polynomu det(A − λI) = 0. Matice A má v komplexním oboru m vlastních hodnot {λ1, . . . , λm} a příslušné vlastní vektory {vλ1 , . . . , vλm } tvoří bázi Cm . Matice T tvořená vlastními vektory (po sloupcích) pak splňuje A · T = T ·    λ1 · · · 0 0 ... 0 0 · · · λm    . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × V případě násobných vlastních hodnot může obsahovat bloky tvaru λ 1 0 λ , přičemž sloupce matice T v tomto případě tvoří tzv. zobecněné vlastní vektory. Jde o vektor splňující Av = λv a další vektor w, který splňuje Aw = λw + v. Pokud je násobnost vlastní hodnoty vyšší než dva, bude se takto vytvářet kaskáda zobecněných vlastních vektorů wi+1 splňující Awi+1 = λwi+1 + wi, která bude spolu s vektorem v tvořit bázi prostoru zobecněných vlastních vektorů. Lineární regulární transformace A → T−1 AT převádí na komplexní Jordanův kanonický tvar. Reálný tvar s reálným blokem α β −β α dostaneme, pokud použijeme místo komplexně sdružených vektorů v a v reálnou a imaginární část u a w vektoru v = u + iw. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Lineární diferenciální systém - opakování Uvažujme lineární diferenciální autonomní systém x′ = Ax, (9) kde x ∈ Rm a A ∈ Rm×m s počáteční podmínkou x(0) = x0. Nechť λ ∈ C je vlastní číslo matice A a v příslušný vlastní vektor. • V případě λ ∈ R je t → eλt v reálným řešením rovnice (9). • V případě λ ∈ R, které je k-násobným kořenem charakteristického polynomu jsou t → eλt i ∑ j=1 ti−jvj (i−j)! , i = 1, . . . k reálnými řešeními rovnice (9), kde vi je systém k zobecněných vlastních vektorů (Av1 = λv1 a Avi = λvi + vi−1 pro i > 1). • V případě λ = α ± iβ je vlastní vektor v = u ± iw a reálnými řešeními rovnice (9) jsou pak t → eαt (cos βt · u − sin βt · w), t → eαt (sin βt · u + cos βt · w). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Protože x0 můžeme zapsat jako lineární kombinaci vlastních vektorů: x0 = k1vλ1 + k2vλ2 + · · · + kmvλm , můžeme řešení x(t) (v případě jednonásobných vlastních čísel, obecně komplexních) zapsat jako x(t) = k1eλ1t vλ1 + k2eλ2t vλ2 + · · · + kmeλ3t vλm . V případě násobných vlastních čísel přibývají k exponenciálním funkcím polynomy. Uvedená řešení jsou lineárně nezávislá a tvoří bázi prostoru řešení. Jejich lineární kombinace je také řešením (9). Maticové zobrazení t → Φ(t) těchto řešení se nazývá fundamentální matice řešení příslušného homogenního lineárního systému (9). Je zřejmé, že rovnovážným bodem systému (9) je počátek, který je stabilní, pokud Re λi < 0 pro všechna i ∈ {1, . . . , m}. Oscilace způsobují komplexní vlastní hodnoty. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Lineární diferenční systém - opakování Uvažujme lineární diferenční autonomní systém xn+1 = Axn, (10) kde xn ∈ Rm , A ∈ Rm×m , n ∈ N0 s počáteční podmínkou x = x0. Odtud xn = An x0. Podobně jako ve spojitém případě má matice A obecně m vlastních hodnot λi, která jsou řešením charakteristické rovnice det(A − λI) = 0. Označme je sestupně |λ1| ≥ |λ2| ≥ · · · ≥ |λm|. Protože x0 můžeme zapsat jako lineární kombinaci vlastních vektorů: x0 = k1xλ1 + k2xλ2 + · · · + kmxλm , ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × můžeme řešení xn zapsat jako xn = An (k1xλ1 + k2xλ2 + · · · + kmxλm ) = k1λn 1 xλ1 + k2λn 2 xλ2 + · · · + kmλn mxλm = λn 1 k1xλ1 + k2 λ2 λ1 n xλ2 + · · · + km λm λ1 n xλm Pevným bodem systému (10) je počátek, který je stabilní, pokud |λ1| < 1. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Lineární diferenciální a diferenční rovnice bývají často zapsány ve tvaru 0 = amx(m) (t) + am−1x(m−1) (t) + · · · + a0x(t), resp. 0 = amxn+m + am−1xn+m−1 + · · · + a0xn. V takovém případě hledáme vlastní čísla jako kořeny charakteristického polynomu p(λ) = amλm + am−1λm−1 + · · · + a0. Poznámka 8. Podkud je levá strana rovnice nenulová, tj. ve tvaru f (t) = . . . (nehomogenní rovnice), pak obecné řešení nehomogenní rovnice je součtem libovolného partikulárního řešení nehomogenní rovnice a obecného řešení příslušné lineární homogenní rovnice (s nulovou levou stranou). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Polynom p(λ) je ve skutečnosti charakteristickým polynomem det(A − λI) = 0 lineárního systému y′ 1(t) = y2(t), ... y′ m−1(t) = ym(t), y′ m(t) = − 1 am (am−1ym(t) + · · · + a0y1(t)), kde y1(t) = x(t) ve spojitém případě. Podobně pro diskrétní případ y1 n+1 = y2 n, ... ym−1 n+1 = ym n , ym n+1 = − 1 am (am−1ym n + · · · + a0y1 n), kde y1 n = xn. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 11. příklad: Dokažte uvedené tvrzení pro 0 = ax′′ + bx′ + cx, resp. 0 = axn+2 + bxn+1 + cxn , tj. ukažte, že kořeny p(λ) jsou vlastní čísla Jacobiho matice jistého dvourozměrného lineárního systému. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Lineární diskrétní model v rovině Samuelsonův model interakce multiplikátoru a akcelerátoru Chceme zjistit jak ovlivňuje GNP multiplikační a akcelerační princip. Multiplikačním efektem rozumíme to, že růst vládních výdajů vede k růstu GNP. Akcelerační efekt je růst investic díky růstu GNP. Koncepce: Proměnnými budou jistě vládní výdaje G a hrubý národní produkt Y, který je součem investic I, spotřeby C a vládních výdajů G. Uvažujeme uzavřenou ekonomiku. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: G Y I C α β ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: Yt = It + Ct + G, Ct = αYt−1, It = β(Ct − Ct−1), kde α ∈ (0, 1) je sklon ke spotřebě, β > 0 je míra růstu investic. GNP Yt, spotřeba Ct a investice It jsou stavové proměnné, G je exogenní proměnná, α a β jsou parametry. Jde o dynamický diskrétní model. Sloučením rovnic dostáváme lineární nehomogenní diferenční rovnici 2. řádu pro Y: Yt+2 − α(β + 1)Yt+1 + αβYt = G (11) Pevný bod Y∗ (rovnováha) splňuje Y∗ − α(β + 1)Y∗ + αβY∗ = G, tj. Y∗ (1 − α) = G ⇒ Y∗ = G 1−α . Dostáváme multiplikační efekt, růst vládních výdajů vede k růstu rovnováhy Y∗ , multiplikátor je 1 1−α . To je jednoduchá komparativní statika. Nás ale bude tentokrát zajímat ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × dynamika systému. Dynamika je dána lineární nehomogenní diferenční rovnici 2. řádu (11). Příslušná homogenní rovnice má charakteristický polynom λ2 − α(β + 1)λ + αβ = 0 s vlastními hodnotami λ1,2 = α(β + 1) ± α2(β + 1)2 − 4αβ 2 . Podle věty o stabilitě diskrétního systému je rovnováha Y∗ stabilní, pokud platí |λ1,2| < 1, tj. α(β + 1) ± α2(β + 1)2 − 4αβ 2 < 1. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 12. příklad: Ukažte, že postačující podmínkou stability Y∗ je αβ < 1. 13. příklad: Napište obecné řešení rovnice (11). Ukažte, že osciluje pro α < 4β (β + 1)2 . 14. příklad: Vyšetřete průběh funkce α = 4β (β + 1)2 . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Následující graf ukazuje oblasti stability a nestability, resp. oscilací rovnováhy Y∗ . Vyhodnocení: Samuelsonův model multiplikátoru-akcelerátoru je prvním modelem, který vysvětluje princip vzniku oscilací GNP. Taky za něj (nejen za něj :-) ) dostal Paul Samuelson v roce 1970 Nobelovu cenu za ekonomii. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Strukturovaný diskrétní dynamický lineární model populace Leslieho model. Model věkově strukturované populace. Můžeme jej použít např. pro modelování populace víceleté rostliny, populace ryb nebo i lidí. Obecně je tedy účelem modelu znát (diskrétní) vývoj struktury populace. Koncepce: Proměnnými budou jistě jednotlivé věkové třídy populace: x1 , . . . xm . Populace se kontroluje po určitých pevných intervalech. Některé skupiny produkují nové jedince, a to s různou mírou reprodukce bi > 0 (dospělí jedinci), jiné mají míru reprodukce nulovou, bi = 0 (nedospělí jedinci). Po nějakém čase přechází určitá část dané třídy xi do následující třídy xi+1 (tyto míry přežití označíme pro každou třídu ci.) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: x1 x2 · · · xm c1 c2 cm−1 b1 b2 bmRovnice:      x1 n+1 x2 n+1 ... xm n+1      =        b1 b2 b3 · · · bm−1 bm c1 0 0 · · · 0 0 0 c2 0 · · · 0 0 ... ... ... ... ... ... 0 0 0 · · · cm−1 0             x1 n x2 n ... xm n      Dostáváme lineární systém diferenčních rovnic s Leslieho maticí L a vektorem iterací struktury populace xn = (x1 n, x2 n, . . . , xm n ), tj. xn+1 = Lxn. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Definice: Nechť A ∈ Rm×m je matice a λ1, . . . λm její vlastní čísla. Striktně dominantní vlastní hodnotou λi matice rozumíme kladnou reálnou vlastní hodnotu jednoduché násobnosti, pro kterou platí |λj| < λi, i = j. Věta (Speciální případ Perronovy - Frobeniovy věty): Předpokládejme, že pro matici L a 1 ≤ i ≤ m platí: bi ≥ 0, existuje nějaké i tak, že bi > 0 a bi+1 > 0, a 0 ≤ ci ≤ 1. Pak má matice L tzv. striktně dominantní vlastní hodnotu λ > 0 a jí příslušný vlastní vektor xλ má všechny složky kladné. Poznámka 9. Protože je λ striktně dominantní, bude pro velká n xn ≈ kλn xλ. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Věková struktura populace se tedy stabilizuje proporcionálně vlastnímu vektoru xλ. Procentní vyjádření je tedy dáno normalizovaným vektorem P = xλ |xλ| , kde výrazem |xλ| rozumíme součet (kladných) složek vektoru xλ. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Model lze prakticky ověřit a je používán nejen pro projekci budoucí struktury populace, ale také například pro kontrolu dynamického systému (odhad trvale udržitelného rybolovu, kácení lesního porostu, pěstování víceletých rostlin apod.). 15. příklad: Uvažujte populaci žen ve věkovém rozmezí 0-14, 15-29, 30-44 a více let. Vysvětlete následující diagram, zvolte stavové proměnné, predikujte situaci za 30 let s počátečními podmínkami danými tabulkou a odhadněte dlouhodobou strukturu populace. 0-14 15-29 30-44 45 a více 1200 1500 1000 1300 Výpočet v Maplu x1 x2 x3 x4 0.99 0.9 0.7 0.5 2 1 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Řízení systému, model těžby: Modifikujme nyní předchozí model a uvažujme nyní řízený systém, kdy populaci částečně vytěžujeme. Může jít o pěstování rostlin, lov ryb, těžbu dřeva nebo o kontrolu populace škůdců apod. Buď D =      d1 0 · · · 0 0 d2 · · · 0 ... ... ... ... 0 0 · · · dm      matice vytěžování, 0 ≤ di ≤ 1. Rovnice modelu má tedy nyní tvar xn+1 = (I − D)Lxn. Naší snahou je udržitelná těžba a stabilizace populace na úrovni x, tj. x = (I − D)Lx, ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × kde x odpovídá vlastnímu vektoru matice (I − D)L příslušnému vlastní hodnotě λ1 = 1. 16. příklad: Najděte podmínku pro udržitelnou těžbu d v případě, že di = d pro všechna i, tj. těžba je věkově nezávislá - rovnoměrná, a je-li λ1 striktně dominantní vlastní hodnota Leslieho matice L. 17. příklad: Uvažujme populaci ryb s Leslieho maticí L =   0 4 3 0.5 0 0 0 0.25 0   . Můžeme zvolit rovnoměrný výlov nebo výlov některé věkové skupiny. Je některá z variant výlovu dlouhodobě udržitelná? Výpočet v Maplu ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Nelineární dynamika a linearizace Uvažujme nyní znovu rovnici (3) resp. (4) x′ = f(x) resp. xn+1 = f(xn) a hyperbolický rovnovážný bod x∗ ∈ Rm (nemá vlastní číslo s nulovou reálnou částí pro spojitý, resp. na jednotkovém kruhu pro diskrétní případ). Podobně jako v jednorozměrném případě můžeme v okolí x∗ funkci f aproximovat Taylorovým rozvojem f(x) ≈ f(x∗ ) + Df(x∗ )(x − x∗ ) + . . . . V dostatečně blízkém okolí x∗ tedy platí x′ ≈ Df(x∗ )(x − x∗ ) resp. xn+1 − x∗ ≈ Df(x∗ )(xn − x∗ ) a nelineální systém (3) resp. (4) se chová v okolí x∗ ”stejně”, jako jeho linearizace. Slovem stejně rozumíme topologickou ekvivalenci (nebudeme dále rozebírat), v prvé řadě jde o lokální stabilitu nebo nestabilitu rovnováhy. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Věta (Věta o linearizaci): Mějme systém (3) resp. (4) s f hladkou v okolí hyperbolického rovnovážného bodu x∗ a jeho linearizaci. Pak v okolí x∗ jsou tyto systémy topologicky ekvivalentní, zejména platí: • Jestliže mají ve spojitém případě všechny vlastní hodnoty matice Df(x∗ ) záporné reálné části, v diskrétním případě jsou-li všechny vlastní hodnoty v absolutní hodnotě menší než 1, pak je x∗ asymptoticky stabilní. • Jestliže ve spojitém případě má alespoň jedna vlastní hodnota matice Df(x∗ ) kladnou reálnou část, v diskrétním je-li alespoň jedna vlastní hodnota v absolutní hodnotě větší než 1, pak je x∗ nestabilní. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Poznámka 10. Charakteristický polynom v rovině. Uvažujme dvourozměrný systém (3) resp. (4), tj. x = (x1, x2)T ∈ R2 . Označme J = Df(x∗ ) Jacobiho matici. Matice J má pak dvě vlastní hodnoty λ1, λ2, které jsou kořeny charakteristické rovnice det(J − λI) = λ2 − σλ + ∆ = 0, kde σ = tr J = λ1 + λ2 je stopa Jacobiho matice a ∆ = det J = λ1λ2 je její determinant. Věta: Postačujícími podmínkami asymptotické stability rovnováhy x∗ spojitého systému (3) v rovině jsou podmínky ∆ = det J > 0 a σ = tr J < 0, kde J = Df(x∗ ) je Jacobiho matice f v rovnovážném bodě. 18. příklad: Dokažte! ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Topologická klasifikace hyperbolického rovnovážného bodu v rovině: stabiln´ı (n+, n−) nestabiln´ı nestabiln´ı (0, 2) (1, 1) (2, 0) F´azov´y portr´et StabilitaVlastn´ı hodnoty uzel ohnisko uzel ohnisko sedlo ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Věta: Postačujícími podmínkami asymptotické stability rovnováhy x∗ diskrétního systému (4) v rovině jsou podmínky |∆| = | det J| < 1, 1 − σ + ∆ = 1 − tr J + det J > 0 1 + σ + ∆ = 1 + tr J + det J > 0, kde J = Df(x∗ ) je Jacobiho matice f v rovnovážném bodě. 19. příklad: Dokažte! ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Stabilita hyperbolického rovnovážného bodu v rovině: ∆ σ σ ∆ 1 10 nestabiln´ı nestabiln´ı nestabiln´ı sedlo sedlo sedlo stabiln´ı stabiln´ı spojit´y syst´em diskr´etn´ı syst´em ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Dynamické modely v rovině V této části použijeme poznatky z kapitoly o systémech diferenciálních a diferenčních rovnic v rovině a aplikujeme je na některé jednoduché modely. Navážeme na statický herní model Cournotova duopolu přidáním dynamiky (spojité i diskrétní), lineární model modifikujeme na nelineární. Uvedeme slavný Samuelsonův model multiplikátoru-akcelerátoru a strukturovaný epidemiologický model SI rozšíříme o další vztahy a přechody mezi skupinami. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Dynamika Cournotova modelu duopolu - spojitý přístup Modelujme nyní dynamiku dříve uvedeného statického herního modelu duopolu. Jde o revizi modelu, kdy si uvědomujeme, že změnit množství výroby směrem k optimu zabere určitý čas a výroba bude klesat nebo růst postupně. Koncepce: Připomeňme Cournotův model. Exogenními proměnnými jsou poptávané množství M a mezní náklady c na výrobu jednoho výrobku, endogenní proměnné jsou množství qi výrobků od jednotlivých firem. Model je dynamický, a proto qi = qi(t) jsou funkcí času t. Poptávková funkce je tvaru P(Q) = M − Q , kde Q = Q(t) = q1(t) + q2(t) je celkové množství dodávané na trh. Produkt je homogenní, mezní náklady c jsou konstantní, tj. nákladová funkce Ci(qi) = cqi(t), c < M. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Výplatní funkcí je zisková funkce firem: πi(qi, qj) = qi P(Q) − c = qi M − (q1 + q2) − c . Nashova rovnováha řeší optimalizační problém max 0≤qi<∞ πi(qi, q−i) = max 0≤qi<∞ qi M − (q1 + q2) − c , tj. platí ∂πi ∂qi = M − q−i − c − 2qi = 0. Optimem je tedy qi = M − c − q−i 2 . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: Změnit množství výroby směrem k optimu qi zabere určitý čas. Budeme předpokládat, že rychlost změny bude přímo úměrná rozdílu mezi optimem a skutečnou výrobou, tj. q′ i = βi(qi − qi), kde βi > 0 je koeficient změny. q′ 1 = β1 M − c − q2 2 − q1 , q′ 2 = β2 M − c − q1 2 − q2 . (12) 20. příklad: Najděte stacionární bod systému (12) (a ukažte, že skutečně existuje), spočtěte pro něj Jacobiho matici a určete jeho typ. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 21. příklad: Nakreslete fázový portrét systému (12) a ukažte, že rovnováha statického Cournotova modelu duopolu je globálně asymptoticky stabilní. Vyhodnocení: Dynamický model oproti statickému modelu ukazuje navíc jakým způsobem se systém rovnováze přibližuje a že k tomu skutečně dochází při jakémkoliv počátečním stavu. Při znalosti odhadu parametrů může sloužit také k odhadu doby, za kterou dojde k dosažení vhodně blízkého okolí této rovnováhy. 22. příklad: V některém z dříve používaných programů vytvořte simulaci a zkoumejte vliv exogenních proměnných a parametrů na dynamiku modelu. V programu XppAut spusťte soubor cournot.ode ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × To, že jsme získali dynamickou stabilní rovnováhu, není nic překvapujícího. Vzpomeneme-li si na chladnoucí kávu, musíme připustit, že jsme použili model přesně kopírující tuto klasickou ukázku implicitně předpokládané stabilní rovnováhy. 23. příklad: Uvažujte revizi tohoto dynamického Cournotova modelu. Předpokládajte racionální chování firem tak, že budou měnit výrobu v závislosti na změně zisku. Čím větší je z navýšení výroby profit, tím ochotněji budou výrobu navyšovat a naopak. Spusťte cournotspojity.mw ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Dynamika Cournotova modelu duopolu - diskrétní přístup Modelujme nyní znovu dynamiku statického herního modelu duopolu, tentokrát diskrétně. Rovnice: Změnit množství výroby směrem k optimu qi zabere určitý čas. Budeme předpokládat, že firma bude měnit množství výroby směrem k optimu, tj. qi(t + 1) = (1 − αi)qi(t) + αiqi(t), kde αi ∈ (0, 1 je rychlost adaptace. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × q1(t + 1) = (1 − α1)q1(t) − α1 2 q2(t) + α1(M−c) 2 , q2(t + 1) = −α2 2 q1(t) + (1 − α2)q2(t) + α2(M−c) 2 . (13) 24. příklad: Najděte pevný bod systému (13) a diskutujte jeho stabilitu. 25. příklad: Porovnejte rovnice spojitý a diskrétní případ a pokuste se z diskrétního přejít ke spojitému limitním přechodem. Vysvětlete, co v obou případech znamenají parametry αi a βi. 26. příklad: Předpokládejme, že se firma řídí mezním ziskem s rychlostí adaptace βi: qi(t + 1) = qi(t) + βiqi(t)∂πi ∂qi (q1(t), q2(t)). Najděte pevné body a diskutujte jejich stabilitu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Další epidemiologické modely Modelujme dynamiku SIR a SIRS modelu. Model SIR. Koncepce: Chceme modelovat epidemii infekční nemoci, kdy jedinci infikovaní přecházejí do skupiny již uzdravených (recovered) a imunních (případně umírají). Jde o nejjednodušší model SIR, Kermack-McKendrickův model s následujícím diagramem: Diagram: S I R β ν ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: S′ = −βSI I′ = βSI − νI R′ = νI, (14) kde β, ν > 0 jsou parametry a S(t), I(t), R(t) stavové proměnné reprezentující okamžitý počet náchylných, infekčních a odolných jedinců v čase. Předpokládáme, že populace se v čase nemění S(t) + I(t) + R(t) = N > 0. (15) a S(0) = S0 > 0, I(0) = I0 > 0, R(0) = 0, S0 + I0 = N. 27. příklad: Z (15) vyjádřete R(t) a zjednodušte model (14) na dvourozměrný se stavovými proměnnými S a I. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 28. příklad: Ukažte, že pokud βS(0) < ν, infekce se vytratí. Zavádíme proto prahovou hodnotu ν β , kterou musí počáteční populace náchylných překročit, aby se epidemie začala šířit. 29. příklad: Najděte stacionární body dvojrozměrného systému. Pokuste se nakreslit fázový portrét s pomocí dI dS . 30. příklad: Spočteme druhou derivaci d2I dS2 a ukažte, že trajektorie jsou konkávní a I nabývá své maximální hodnoty Imax = ν β (ln ν βS0 − 1) + S0, pro Smax = ν β . 31. příklad: Ukažte, že platí S(t) = S(0)e− βR(t) ν a proveďte limitní přechod t → ∞, abyste nalezli rozsah infekce daný mírou R(∞) N = 1 − S(∞) N = 1 − ρ. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 32. příklad: Ve vhodném programu simulujte model SIR. Vyhodnocení: Výstupy z modelu jsou v souladu s realitou. Chceme-li omezit rozsah epidemie, je třeba zvětšit ρ, je tedy potřeba zvýšit rychlost izolace infikovaných jedinců (snížit koeficient β) a zvýšit odolnost jedinců vůči nakažení infekcí při kontaktu s infikovaným jedincem. Navíc získáváme další důležité epidemiologické informace: prahovou hodnotu, maximum infikovaných jedinců apod. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model SIRS. Koncepce: Chceme modelovat epidemii infekční nemoci, kdy infikovaní jedinci přecházejí do skupiny uzdravených (recovered), oproti předchozímu modelu však nezůstávají imunní a mohou znovu onemocnět. Diagram: S I R β ν γ 33. příklad: Sestavte rovnice modelu SIRS pro konstatní populaci (tj. při splnění podmínek (15)). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 34. příklad: Podobně jako v modelu SIR přejděte na dvourozměrný se stavovými proměnnými S a I. 35. příklad: Najděte stacionární body dvojrozměrného systému, spočtěte zde Jacobiho matici a určete jejich stabilitu. Pokuste se nakreslit fázový portrét. Vyhodnocení: Model můžeme samozřejmě dále rozšiřovat: 36. příklad: Nalezněte na internetu nějaké informace o modelu SIRS. Pokuste se najít nějaký vědecký článek, který je jeho rozšířením. Vytvořte ve vhodném programu simulaci. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Dynamika chemických reakcí Chemické a biochemické reakce je vhodné popisovat pomocí diferenciálních rovnic. Elementární reakce podléhají kinetické rovnici, která popisuje rychlost, se kterou interagují dvě látky a vytvářejí třetí: A + B k → C Koncentrace látek se značí v hranatých závorkách a uvedenou reakci můžeme popsat rovnicí d[C] dt = k[A][B], kde derivace koncentrace [C] je okamžitá změna koncentrace [C], tedy rychlost, s jakou je tvořen produkt reakce. Konstanta k je rychlostní konstanta, která vlastně konstantou není – závisí např. na teplotě nebo homogenitě směsi. Budeme ale předpokládat, že se teplota nemění a látky jsou dobře promíchané. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Většina biochemických reakcí probíhá oběma směry: A + B k+ ⇄ k− C Změna koncentrace [A] pak splňuje d[A] dt = −k+[A][B] + k−[C]. Ve skutečnosti je většina reakcí složitější a bude tedy popsána systémem diferenciálních rovnic. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model Michaelise-Mentenové Koncepce: Enzymy E jsou katalyzátory chemických reakcí, při kterých pomáhají ze substrátu S vytvořit produkt P, přičemž z reakce vycházejí samy v nezměněné formě. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: Substr´at Enzym Komplex enzym-substr´at Produkt E + S k1 ⇄ k−1 C k2 → E + P ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: Kinetické rovnice reakcí tedy můžeme popsat následujícími diferenciálními rovnicemi: d[S] dt = k−1[C] − k1[S][E], d[E] dt = (k−1 + k2)[C] − k1[S][E], d[C] dt = k1[S][E] − (k2 + k−1)[C], d[P] dt = k2[C]. Navíc předpokládáme, že produkt P okamžitě odebíráme, aby nešel do zpětné reakce. Je evidentní, že platí d[E] dt + d[C] dt = 0, tj. [E] + [C] = e0 je počáteční koncentrace enzymu, [E] tedy můžeme eliminovat. Rovnici produktu můžeme oddělit a integrovat zvlášť. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Označme [S] = s a [C] = c. Úpravou tedy dostáváme dvě diferenciální rovnice: ˙s = k−1c − k1s(e0 − c), ˙c = k1s(e0 − c) − (k2 + k−1)c s počátečními podmínkami c(0) = 0 a s(0) = s0 >> e0. 37. příklad: Dokažte, že počátek je asymptoticky stabilní rovnovážný bod. 38. příklad: Nakreslete fázový portrét a graficky analyzujte systém a nakreslete přibližně tvar řešení s uvedenou počáteční podmínkou. 39. příklad: Simulujte řešení ve vhodném programu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Z výsledků je zřejmé, že koncentrace komplexu c nejprve roste ke své maximální hodnotě a pak monotónně klesá k nule. Tato maximální hodnota je cmax = e0s K + s , kde K = k2+k−1 k1 je tzv. Michaelisova konstanta. Vzhledem k tomu, že pro počáteční koncentrace enzymu a substrátu platí e0 << s0, je trajektorie řešení velmi rychle přitahována k ˙c = 0 nulklině, kterou následně ”kopíruje”, tj. změna koncentrace komplexu je téměř stálá. Chemici tomuto říkají kvazi-stacionární stav nebo kvazi-rovnováha, kdy platí ˙c = k1s(e0 − c) − (k2 + k−1)c = 0 Jde o jeden ze základních bichemických dynamických modelů. Tento model vznikl na počátku minulého století a je dodnes hojně využíván. Ve složitějších modelech bichemických reakcí v buňkách jsou právě ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Michaelisovy konstanty různých dílčích katalytických reakcí vstupujících do dynamiky systému parametry. Kvazi-rovnováh se pak využívá pro popis složitějších enzymatických reakcí (replikace DNA, dělení buněk apod.), přičemž se předpokládá, že komplex splňuje výše uvedenou podmínku, tj. c(t) = e0s(t) K + s(t) . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Dynamické modely interakcí Snad nejznámějším deterministickým dynamickým modelem je model interakce dravec-kořist. Tím nejjednodušším je Lotkův-Volterrův model, který stojí u základů vědní disciplíny zvané matematická ekologie. Model dravec-kořist. Koncepce: Modelujme dvě vzájemně provázané populace - populaci kořisti x a dravce y. Je zřejmé, že velikost a dynamika populace kořisti bude ovlivňovat dynamiku a velikost populace dravce a naopak. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: x y f(x,y) g(x,y) Rovnice: x′ = x f (x, y) y′ = yg(x, y), (16) kde stavové proměnné x a y reprezentují populace kořisti a dravce a f (x, y) = r − λy a g(x, y) = eλx − d pro parametry r, λ, e, d > 0. 40. příklad: Interpretujte parametry modelu (16). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 41. příklad: Najděte stacionární body systému (16), spočtěte zde Jacobiho matici a určete jejich stabilitu. Lze použít větu o linearizaci? Pokuste se nakreslit fázový portrét. 42. příklad: Najděte předpis netriviální trajektorie (16) ve fázovém prostoru s počáteční podmínkou y(x0) = y0, použijte znalost toho, že y′ x′ = dy dx . Řešte samozřejmě v 1. kvadrantu, který je pro model smysluplný. Může trajektorie tento kvadrant opustit? Jak trajektorie vypadá? 43. příklad: Srovnejte Lotkův-Volterrův model s Kermack-McKendrickovým modelem SIR. 44. příklad: Podívejte se na Scholarpedii ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Model můžeme samozřejmě dále rozšiřovat: 45. příklad: Nalezněte na internetu nějaké informace o modelu dravec-kořist. Pokuste se najít nějaký vědecký článek, který je jeho rozšířením. Vytvořte ve vhodném programu simulaci. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Revize Lotkova-Volterrova modelu úpravou dynamiky populace kořisti: Míra růstu kořisti je v Lotkově-Volterrově modelu konstantní, bez přítomnosti predátora se kořist bude množit exponenciálně. Revidovat můžeme např. zavedením kapacity prostředí nebo prahu přežití. 1. r(x) = r 2. r(x) = r(1 − x K ), kde K je kapacita prostředí 3. r(x) = r(1 − x K )( x A − 1), kde K je kapacita prostředí a A práh přežití (tzv. silný Alleeho efekt) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × A K r 1. exponenci´aln´ı 2. logistick´y 3. slab´y a siln´y Alleeho efekt r(x) x ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Revize Lotkova-Volterrova modelu úpravou funkce predace: Predace je v Lotkově-Volterrově modelu přímo úměrná velikosti (hustotě) populace kořisti. Jeden predátor vyhledá (a uloví) za čas ∆ts ∆x = λx∆ts jedinců kořisti. Okamžitá změna množství kořisti ulovená jedním predátorem, tedy jakási schopnost lovu, se nazývá funkční odpověď predátora. V případě Lotkova-Volterova modelu je tato funkční odpověď Φ(x) = λx. Mluvíme o lineární funkční odpovědi. Pokud lineární funkční odpověď ohraničíme hladinou nasycení, dostaneme funkční odpověď Hollingova I. typu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Pokud uvažujeme, že predátor po nalezení kořisti potřebuje k jejímu ulovení a strávení nějaký další čas h (handling time), pak čas na vyhledání kořisti je zkrácen o tuto dobu, tj. ∆ts = ∆tt − h∆x. Dosazením pak ∆x = λx(∆tt − h∆x), ∆x + λxh∆x = λx∆tt, ∆x = λx 1+hλx ∆tt. Funkční odpověď predátora je pak Hollingova II. typu Φ(x) = λx 1+hλx, která pro h = 0 odpovídá predaci v předchozím Lotkově-Volterrově modelu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Další možností revize je predace více než přímo úměrně závislá na populaci kořisti. Odtud pak podobně vyplývá funkční odpověď Hollingova III. typu Φ(x) = λxk 1+hλxk , pro k > 1. Typu I. odpovídají např. dravci, kteří se krmí filtrací (planktonu, bakterií, hmyzu apod), typu II. odpovídá většinou hmyz a paraziti, typu III. pak obratlovci. Jako důvod je nejčastěji uváděn proces učení lovu, který je v populaci o malé hustotě daleko pomalejší. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Φ(x) x I. II. III. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Goodwinův model hospodářského cyklu Představíme si nyní ekonomický model interakcí, který vede na model dravec-kořist. Koncepce: Budeme vycházet z Harrodova-Domarova modelu, přičemž budeme předpokládat, že veškerá čistá produkce, tj. produkce bez vyplacených mezd, je investována. Označme L množství zaměstnaného obyvatelstva, které za svou práci dostává mzdu W (jde vlastně o střední hodnotu mzdy). N bude množství práceschopného (nebo práceochotného) obyvatelstva. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Pro zjednodušení zavedeme následující veličiny: • produktivita práce = střední množství produktu vytvořeného jedním pracujícím člověkem a = Y L , • relativní zaměstnanost v = L N , • podíl mzdy na produkci u = W a = WL Y . Dále předpokládejme, že míra růstu obyvatel β je konstantní, projevuje se stálý technický pokrok, tj. konstantní relativní růst produktivity práce α a relativní změna mzdové sazby závisí na relativní zaměstnanosti. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Phillipsova křivka: Závislost relativní změny mzdové sazby na relativní zaměstnanosti (nebo nezaměstnanosti) popisuje Phillipsova křivka, jejíž vlastnosti byly zjištěny empiricky. Funkce ϕ : [0, 1) → R je diferencovatelná funkce, která je rostoucí a konvexní a splňuje nerovnosti ϕ(0) < 0, (17) tj. při malé zaměstnanosti (velké nezaměstnanosti) mzdy klesají (je-li práce vzácná, lidé jsou ochotni pracovat za nízkou mzdu), lim v→1− ϕ(v) > 0, (18) tj. při velké zaměstnanosti mzdy rostou (chceme-li při téměř plné zaměstnanosti získat nového pracovníka, musíme ho přeplatit). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Phillipsova křivka jako závislost ˙W W na relativní nezaměstnanosti, tedy jako funkce 1 − v, je klesající konvexní funkce (otočení okolo osy v = 1 2 ). Někdy se místo relativní změny mzdové sazby analogicky vyjadřuje inflace. 0 1 0 1 v ˙W W ˙W W ϕ(v) m´ıra nezamˇestnanosti m´ıra zamˇestnanosti ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: Oproti Harrodovu-Domarovu modelu nyní platí I = Y − LW, tedy při původním označení kapitálové náročnosti jednotky produkce r = K Y můžeme změnu kapitálu psát jako K′ = rY′ = Y − LW − δrY. Odtud Y′ Y = 1 r (1 − u) − δ. 46. příklad: Ukažte, že za daných předpokladů platí v′ v = 1 r (1 − u) − δ − α − β. 47. příklad: Za předpokladu W′ W = ϕ(v) ukažte, že platí u′ u = ϕ(v) − α. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice v′ = v(1 r (1 − u) − δ − α − β) u′ = u(ϕ(v) − α) odpovídají modelu dravec-kořist, kde v je kořist a u dravec. Vyhodnocení: Goodwinův model je významný model vysvětlující endogenní fluktuace ekonomiky. Stejně jako model dravec-kořist přispěl k pochopení, že oscilace nemusí být vyprovokovány vnějšími periodickými vlivy a že oscilace nejsou o patologickým jevem, ať už jde o interakce v biosystémech, ekosystémech, chemických reakcích či v ekonomii. Pochopení vzniku oscilací a jiných nerovnovážných stavů v dynamických systémech dává tušit mnoho nového. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Například moderní pojetí termodynamiky jako obecně nerovnovážné dynamiky otevřených systémů vysvětluje na základě modelů s limitními cykly (nebo jinými i chaotickými atraktory) vznik složitých struktur a jejich uspořádání. Klasická termodynamika popisuje především izolované systémy, které po určitém čase dosáhnou rovnováhy (teplota nerovnoměrně zahřátého tělesa se časem vyrovná, sůl ve skleničce vody se časem rozpustí a koncentrace solného roztoku se vyrovná). Mluvíme o dosažení maximální entropie. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Pro otevřené systémy je situace odlišná. Otevřené systémy, kterými jsou i živé organismy, využívají okolní energii na udržení a zvýšení své vlastní uspořádanosti a snižují takto entropii systému. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × U vzniku teorie nerovnovážné termodynamiky, která vysvětluje tento jev samoorganizace, stojí především dva laureáti Nobelovy ceny za chemii Lars Onsager (1968) a Ilya Prigogine (1977). Výsledkem je vznik nových vědních oblastí čerpajících z nerovnovážné teorie. Často je najdeme jako aplikace matematické disciplíny nazývané nelineární dynamika nebo teorie bifurkací a teorie chaosu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyberme namátkou biochemii (pochopení např. biochemických přepínačů u enzymatických reakcí v buňkách), ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × ekonomii (neviditelná ruka trhu vedoucí k rovnovážnému stavu již není dogmatem), ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × biologii (vznik vzorů u zvířat), Watanabe M, Kondo S.: Changing clothes easily: connexin41.8 regulates skin pattern variation, Pigment Cell Melanoma Res. 2012 Feb 7 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × neurovědu (např. popis chování neuronů), ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × meteorologii (dynamika počasí), fyziku (např. jevy hydrodynamického proudění), psychologii, sociologii, dopravu až po kosmologii a další. Těm, kteří mají zájem o tuto oblast matematiky doporučuji předmět M6201 Nelineární dynamika a její aplikace vyučovaný příští rok. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Evoluční hry Teorie evolučních her je spojením teorie her a dynamických systémů v biologických aplikacích. Strategie v teorii her byly strategiemi rozumných hráčů, ti jednali na základě optimalizace (Nashova rovnováha). Jistě se nedá očekávat, že se podobně racionálně budou chovat zvířata nebo rostliny (ani lidé to často nedělají). Přesto jisté strategie v přírodě ”vyhrávají”ve smyslu přežití, říkáme jim evolučně stabilní strategie. Teorie evolučních her je vcelku nová disciplína rozvíjející se každým dnem, proto si ukážeme jen některé její hlavní principy na nejznámějším modelu ”hawk and dove”a na tomto jednoduchém modelu si odvodíme některé obecnější věty. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Co je v evoluci hrou a co je evolučně stabilní strategie? V evoluční hře budou hráči subpopulace druhu s možnou strategií, tj. fenotypem chování. Výplatní funkcí je fitness (biologická, reprodukční zdatnost), která představuje zdatnost zachovat své geny a rozšířit je v genotypu populace (genotypem rozumíme soubor všech genů, které má organismus k dispozici). Klasická darwinistická a neodarwinistická teorie evoluce předpokládá, že kritériem evolučního úspěchu jedince je jeho fitness. Podle klasických představ by se v důsledku přirozeného výběru měly v populaci zachovat ty geny, které svému nositeli poskytují největší fitness. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × S nástupem evoluční teorie her se však ukázalo, že z dlouhodobého hlediska není důležité, jak příslušný gen mění fitness nositele, ale to, jestli podmiňuje evolučně stabilní strategii, tj. takovou strategii, která pokud jednou v populaci převládne, nemůže být potlačena žádnou jinou minoritní strategií. Od statického modelu tedy musíme nutně přejít k dynamickému. V evoluci nevyhrává vždy zdatnější, ale stabilnější... ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model jestřáb a holubice Koncepce: Uvažujme dvě populace jednoho druhu H a D bojující mezi sebou o zdroj potravy. H představuje fenotyp chování hawk - jestřáb, který bojuje tvrdě, chladnokrevně a vzdává se jen tehdy, když je vážně zraněný, dove D - holubice se uchyluje jen k symbolické hrozbě a při přímém útoku utíká nezraněná. Cílem evoluční teorie her je určit zastoupení těchto dvou strategií v populaci a která z nich převládne. Diagram: strategie H D H 1 2 (G − C) G D 0 1 2 G ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Uvažujme nejprve populaci fenotypu D. Vstoupí-li do ní jedinec fenotypu H, bude se v ní s jistotou šířit, protože fitness - výplatní funkce u(H, D) = G > 1 2 G = u(D, D). D tedy není evolučně stabilní strategie, protože není odolná vůči vstupu mutantního fenotypu. Uvažujme naopak populaci fenotypu H. Vstoupí-li do ní jedinec fenotypu D, bude pro fitness platit u(D, H) = 0 a u(H, H) = 1 2 (G − C) a logicky bude záviset na tom, zda zisk z boje G bude větší nebo menší než náklady na boj C. V případě, že G > C, pak u(D, H) < u(H, H) a mutantní fenotyp se v populaci nebude moci šířit. V této situaci bude fenotyp jestřába H evolučně stabilní strategií. V opačném případě, kdy náklady C převýší zisk G, bude se fenotyp holubice D moci šířit v populaci jestřábů. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Uvažujeme-li smíšenou strategii, kdy se jedinec chová jako jestřáb s pravděpodobností x a jako holubice s pravděpodobností 1 − x, pak fitness jednotlivých fenotypů je daná vztahy u(H, xH + (1 − x)D) = xu(H, H) + (1 − x)u(H, D) = x 1 2 (G − C) + (1 − x)G u(D, xH + (1 − x)D) = xu(D, H) + (1 − x)u(D, D) = x · 0 + (1 − x)1 2 G 48. příklad: Ukažte, že pokud x < G C , dochází k šíření fenotypu H a pokud x > G C , dochází k šíření fenotypu D. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: Uvažujme populaci o H jedincích fenotypu jestřába a D holubice, velikost celé populace je N = H + D. Předpokládáme, že každý fenotyp se rozmnožuje úměrně svému fitness uH = rH(H, D) a uD = rD(H, D), který závisí samozřejmě na zastoupení jedinců fenotypu H a D v populaci: H′ = rH(H, D)H D′ = rD(H, D)D Dynamika růstu celé populace je pak dána rovnicí N′ = rH(H, D)H + rD(H, D)D = rH(H, D) H N N + rD(H, D) D N N = ¯rN, kde ¯r = rH x + rD(1 − x), přičemž x představuje podíl jestřába v celé populaci. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 49. příklad: Ukažte, že platí tzv. replikátorová rovnice x′ = x(rH − ¯r). Pro výplatní matici A = 1 2 (G − C) G 0 1 2 G pak fitness fenotypu jestřába bude rH = xu(H, H) + (1 − x)u(H, D) = 1 2 (G − C)x + G(1 − x), fitness fenotypu holubice bude rD = xu(D, H) + (1 − x)u(D, D) = 1 2 G(1 − x) a ¯r = rH x + rD(1 − x) = x(1 2 (G − C)x + (1 − x)G) + (1 − x)1 2 G(1 − x). Replikátorová rovnice pro fenotyp jestřába je proto x′ = C 2 x(x − 1)(x − G C ). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 50. příklad: Odvoďte :-) 51. příklad: Najděte stacionární body a určete jejich stabilitu. 52. příklad: Určete zastoupení strategií jestřába a holubice v populaci v dlouhodobém horizontu. 53. příklad: Ve vhodném programu simulujte populaci jestřábů a holubice. Předpokládejte náklady na boj ve výši C = 4, zisk G = 1 a počáteční populaci jestřábů a holubic v poměru 1:100. 54. příklad: Porovnejte výsledek dynamického modelu s herním modelem se smíšenými strategiemi. Smíšenou strategii (x, 1 − x)T můžeme v populaci jestřábů a holubic vnímat jako pravděpodobnost chování náhodného jedince (samozřejmě předpokládáme, že každého jedince potkáme se stejnou pravděpodobností, tj. např. holbice se ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × neshlukují). Tato smíšená strategie náhodného jedince je tedy dána právě poměrem fenotypů v populaci. Vyhodnocení: Model je samozřejmě velmi jednoduchý, právě pro svou přehlednost je jedním ze základních modelů biologie, vysvětluje sice dynamiku evoluce pouze dvou fenotypů, ale jeho princip lze použít obecně. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Teorie her a dynamika V této části použijeme předchozí model hawk-dove pro odvození obecnějšího principu. Půjde nám o vyjádření replikátorových rovnic pro populaci složenou z n fenotypů. Definice: Maticovou symetrickou hrou dvou hráčů rozumíme hru se stejnými konečnými n-rozměrnými prostory strategií S1 = S2 = S a symetrickými výplatními funkcemi u1(i, j) = u2(j, i) = (aij), i, j ∈ S. Výplatní funkci pro smíšené strategie x, y ∈ (S) pak zapisujeme pomocí výplatní matice A, tj. u1(x, y) = xT Ay = u2(y, x) pro x, y ∈ (S). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Poznámka 11. Smíšená strategie ¯x je rovnovážnou strategií maticové symetrické hry s maticí A právě tehdy, když pro všechny smíšené strategie x ∈ (S) platí ¯xT A ¯x ≥ xT A ¯x, tj. ¯xT A ¯x = max x∈(S) xT A ¯x. Označme ryzí strategie ei = (0, . . . , 0, 1, 0, . . ., 0)T (vektor s i-tou nenulovou složkou) a pro smíšenou strategii x = (x1, . . . xn)T definujme množinu indexů nenulových pravděpodobností C(x) = {k : xk > 0}. Pak platí následující věta: ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Věta: Smíšená strategie ¯x symetrické maticové hry s maticí A je rovnovážná právě tehdy, když ¯xT A ¯x ≥ eT i A ¯x pro všechna i /∈ C( ¯x) a ¯xT A ¯x = eT i A ¯x pro všechna i ∈ C( ¯x). Poznámka 12. Ryzí strategie ei je tedy rovnovážnou strategií právě tehdy, když pro všechna j platí aii ≥ aji . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Replikátorové rovnice: Uvažujme nyní populaci n fenotypů o velikosti Ni, i = 1 . . . n, rozložení populace je tedy x = (x1, . . . , xn)T , kde xi = Ni N . Výplatní matice A ∈ Rn×n určuje fitness (a tedy růst populace) i-tého fenotypu takto: N′ i = riNi, kde ri = eT i Ax. Růst celé populace je určen průměrnou mírou růstu N′ = ¯rN, kde ¯r = n ∑ i=1 xiri = n ∑ i=1 xieT i Ax = xT Ax. Pro jednotlivé fenotypy tedy platí x′ i = xi(ri − ¯r) = xi(eT i Ax − xT Ax). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 55. příklad: Napište repikátorové rovnice pro symetrickou maticovou hru s maticí A =   0 1 0 0 0 2 0 0 1   . 56. příklad: Ukažte, že simplex x1 + x2 + x3 = 1 je invariantní množinou dynamického systému daného replikátorovými rovnicemi. (Návod: uvažujte dynamiku nové proměnné s = x1 + x2 + x3 pro s = 1.) 57. příklad: Ve vhodném programu nakreslete fázový portrét (2D i 3D), v dvojrozměrném nakreslete nulkliny, spočtěte stacionární body a na základě Jacobiho matice určete jeho stabilitu (pokud to lze). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 58. příklad: Uvažujte interakci dvou populací - prodejců a kupujících. Prodejce se může řídit dvěma strategiemi - buď být čestný, nebo podvádět. Kupující může buď prověřit nebo neprověřit, co kupuje. Jde o tzv. bimaticovou hru (prodávající a kupující mají obecně nesymetrické výplatní matice) P = 2 1 3 4 , K = 3 2 4 1 . Předpokládáme, že prodávající a kupující budou používat danou strategii tím více, čím úspěšnější bude (máme tu jakousi fitness daného fenotypu). Odvoďte replikátorové rovnice a nakreslete jejich fázový portrét. Návod v článku. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Dynamický model difúze a šíření Difúze je proces, při kterém se částice pohybují proti směru gradientu koncentrace. Je výsledkem pohybu mnoha malých částeček v náhodných směrech (Brownův pohyb). Jedním ze způsobů modelování je proto agregační přístup pomocí náhodné procházky. Uvažujme nejjednodušší případ, pohyb jedné částečky po přímce. Za jednotku času se posune náhodně vpravo nebo vlevo o ξ, za čas t urazí vzdálenost x(t) od místa, ve kterém se nacházela v čase t = 0. Pokud budeme sledovat mnoho takových částeček, urazí za čas t průměrně jakousi střední vzdálenost. Platí pro ni, že její čtverec s časem lineárně roste. Ukázal to v roce 1905 Einstein při modelování Brownova pohybu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Intuitivní představa je tato: V počátečním okamžiku umístíme i-tou částečku do počátku a budeme sledovat její vzdálenost xi(t) od počátku po t krocích (v čase t). Po t − 1 krocích se nacházíme ve vzdálenosti xi(t − 1) a po dalším kroku bude naše vzdálenost rovna jedné z následujících možností xi(t) = xi(t − 1) − ξ nebo xi(t) = xi(t − 1) + ξ, tedy (xi(t))2 = (xi(t − 1))2 ± 2xi(t − 1)ξ + ξ2 . Při sledování dostatečně velkého počtu částic průměrně (x(t))2 = (x(t − 1))2 + ξ2 , x(0) = 0. Střední kvadratrická vzdálenost tedy bude růst úměrně času: (x(t))2 = ξ2 t, tj. (x(t))2 t = 2D ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Konkrétně si to představme na případě bakterií ve vodě v tenké trubici. Bakterie je velká asi 10−4 cm. Difúzní koeficient D je asi D = 10−5 cm2 /s O vzdálenost x přibližně rovné své velikosti se posune asi za t = x2 2D = 5 · 10−4 s, což je půl milisekundy. Na vzdálenost jednoho centimetru bakterie difundují ale až za čas t = 12 2·10−5 = 5 · 104 s. To je skoro 14 hodin. A do dvakrát tak velké vzdálenosti jim to zabere 4× tak dlouhý čas. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Pravděpodobnost, že se částečka bude vyskytovat v čase t na místě x bude p(x, t) = 1 2 p(x + ∆x, t − ∆t) + 1 2 p(x − ∆x, t − ∆t), odečtením p(x, t − ∆t) a podělením ∆t dostaneme p(x,t)−p(x,t−∆t) ∆t = (∆x)2 2∆t p(x+∆x,t−∆t)−2p(x,t−∆t)+p(x−∆x,t−∆t) (∆x)2 . Z definice derivace limitním přechodem ∆t → 0 a ∆x → 0 dostaneme rovnici difúze: ∂p ∂t = D ∂2p ∂x2 , kde D = (∆x)2 2∆t = const. je difúzní koeficient. Z matematického hlediska je rovnice totožná s rovnicí vedení tepla. Zvlášť pro více prostorových dimenzí se používá pro operátor na pravé straně rovnice difúze (nebo vedení tepla) označení ∂2 ∂x2 1 + · · · + ∂2 ∂x2 n = ∇2 = ∆ a nazývá se Laplaceův operátor. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Poznámka 13. Operátor nabla je často používán v zápisech vícerozměrných dynamických systémů. Jde vlastně o zkrácený zápis vektoru parciálních derivací, tj. ∇f = ( ∂ f ∂x1 , ∂ f ∂x2 , . . . , ∂ f ∂xn ). Zápis s tečkou pak značí součet složek jako je tomu u skalárního součinu, tj. ∇ · f = ∂ f ∂x1 + ∂ f ∂x2 + · · · + ∂ f ∂xn , Laplaceův operátor ∆ je tedy formálně skalární součin operátorů ∇ · ∇ = ∆. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Spojitý přístup: Uvažujme proudění tenkou trubicí: x x + ∆x S J(x, t) V = S · ∆x J(x, t) je vektor ve směru toku o velikosti hustoty částic (počet částic v čase t na jednotku plochy), S je plocha řezu. Je-li u(x, t) koncentrace částic v (x, t), pak bude v objemu V změna množství částic za ∆t: S ∆x u(x, t + ∆t) − S ∆x u(x, t) ∆t = S(J(x, t) − J(x + ∆x, t)). Dělením S ∆x a limitním přechodem ∆t → 0 dostaneme zákon zachování hmoty: ∂u ∂t = − ∂J ∂x ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Obecněji, pokud by částice v trubici navíc vznikaly nebo zanikaly s hustotou f (x, t) (počet vzniklých nebo zaniklých částic na jednotku času a objemu), byla by rovnice zákona zachování ve tvaru ∂u ∂t = − ∂J ∂x + f. Hustota proudění částic J(x, t) je nejčastěji ovlivňována dvěma jevy advekcí, tj. přenosem částic v médiu proudícím rychlostí v, pak Jadv = vu a difúzí. Difúzní proudění podléhá empirickému Fickovu zákonu Jdi f f = −D ∂u ∂x , tedy tok částic závisí přímo úměrně na změně koncentrace částic a směřuje k vyrovnání koncentrace. V případě čisté difúze tedy platí ∂u ∂t = D ∂2u ∂x2 = D∇2 u. (19) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Prozatím se spokojíme s jednou prostorovou dimenzí a podíváme se blíže na řešení rovnice difúze (19). To závisí jistě na počátečních podmínkách, ale také na dalších podmínkách. Nejčastěji jsou to podmínky okrajové, tj. koncentrace částic je určena na okraji trubice x0 ∈ {0, l}: • u(x0, t) = φ(t) Dirichletovy, • ∂u ∂x (x0, t) = ψ(t) Neumannovy, • ∂u ∂x (x0, t) = −cu(x0, t) Robinovy Rovnice difúze je separovatelná. Hledáme tedy řešení (19) ve tvaru u(x, t) = F(x)G(t). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Dosazením do (19) dostáváme F(x)G′ (t) = D∇2 F(x)G(t), G′(t) DG(t) = ∇2F(x) F(x) = −λ = const. tj. −∇2 F(x) = λF(x), a G′ (t) = −λDG(t). 59. příklad: Stacionární difúze membránou při D = konst.. Uvažujme membránu o tloušťce l, jejíž povrchy jsou udržovány na konstantních koncentracích. Určete stacionární stav difúze na membráně. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Řešení rovnice difúze pomocí Fourierovy transformace. Fourierova transformace převádí funkci F(x) na funkci F(α). Zjednodušeně řečeno můžeme každou funkci zapsat jako ”součet”funkcí sin a cos s nějakými frekvencemi a amlitudami. Periodické funkce můžeme takto rozložit do tzv. Fourierovy řady (lineární kombinace funkcí sin(αx) a cos(αx) s různými frekvencemi α a dostaneme, tzv. spektrum funkce, které určuje kterým frekvencím přísluší která amplituda. Pro neperiodické funkce můžeme přejít od součtu k integrálu. Podstatnou vlastností této transformace je to, že můžeme zpětně z tohoto spektra dostat funkci původní, tj. existuje inverzní transformace. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Řešení rovnice difúze pomocí Fourierovy transformace. Fourierova transformace převádí funkci F(x) na funkci F(α). Zjednodušeně řečeno můžeme každou funkci zapsat jako ”součet”funkcí sin a cos s nějakými frekvencemi a amlitudami. Periodické funkce můžeme takto rozložit do tzv. Fourierovy řady (lineární kombinace funkcí sin(αx) a cos(αx) s různými frekvencemi α a dostaneme, tzv. spektrum funkce, které určuje kterým frekvencím přísluší která amplituda. Pro neperiodické funkce můžeme přejít od součtu k integrálu. Podstatnou vlastností této transformace je to, že můžeme zpětně z tohoto spektra dostat funkci původní, tj. existuje inverzní transformace. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Řešení rovnice difúze pomocí Fourierovy transformace. Fourierova transformace převádí funkci F(x) na funkci F(α). Zjednodušeně řečeno můžeme každou funkci zapsat jako ”součet”funkcí sin a cos s nějakými frekvencemi a amlitudami. x 6 0 5 1 -1 -3 2 -2 0 2 431 Periodické funkce můžeme takto rozložit do tzv. Fourierovy řady (lineární kombinace funkcí sin(αx) a cos(αx) s různými frekvencemi α a dostaneme, tzv. spektrum funkce, které určuje kterým frekvencím přísluší která amplituda. Pro neperiodické funkce můžeme přejít od⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Řešení rovnice difúze pomocí Fourierovy transformace. Fourierova transformace převádí funkci F(x) na funkci F(α). Zjednodušeně řečeno můžeme každou funkci zapsat jako ”součet”funkcí sin a cos s nějakými frekvencemi a amlitudami. 0,6 -400 0,4 400 -200 200 0,2 -600 0 x 0,8 0 -800 1 600 Periodické funkce můžeme takto rozložit do tzv. Fourierovy řady (lineární kombinace funkcí sin(αx) a cos(αx) s různými frekvencemi α⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Řešení rovnice difúze pomocí Fourierovy transformace. Fourierova transformace převádí funkci F(x) na funkci F(α). Zjednodušeně řečeno můžeme každou funkci zapsat jako ”součet”funkcí sin a cos s nějakými frekvencemi a amlitudami. Periodické funkce můžeme takto rozložit do tzv. Fourierovy řady (lineární kombinace funkcí sin(αx) a cos(αx) s různými frekvencemi α a dostaneme, tzv. spektrum funkce, které určuje kterým frekvencím přísluší která amplituda. Pro neperiodické funkce můžeme přejít od součtu k integrálu. Podstatnou vlastností této transformace je to, že můžeme zpětně z tohoto spektra dostat funkci původní, tj. existuje inverzní transformace. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × F(α) = ∞ −∞ F(x)e−iαx dx F(x) = 1 2π ∞ −∞ F(α)eiαx dα Komplexní zápis eiαx skrývá ony periodické funkce podle Eulerova vztahu eix = cos x + i sin x. Fourierova transformace se používá např. při analýze signálů (zjištění vad strojů), mp3, jpeg, MPEG kompresi obrazu a zvuku nebo při lokalizaci systémem GPS. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Hledejme nyní řešení rovnice difúze (19) s okrajovými podmínkami u(±∞, t) = 0, ∂u ∂x (±∞, t) = 0 a počáteční podmínkou u(x, 0) = M S δ(x), kde M je množství látky vstříknuté do trubice o ploše řezu S v počátku x = 0 a v čase t = 0. ”Funkce”δ je Diracova delta funkce, jde ve skutečnosti o tzv. distribuci, což je rozšíření pojmu funkce pomocí limity. Animace Z tohoto limitního pohledu můžeme napsat ”definici”Diracovy delta funkce takto: δ(x) = ∞ : x = 0 0 : x = 0 přičemž ∞ −∞ δ(x)dx = 1. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Formálně korektně bychom museli definovat δ jako míru množiny A ⊂ R, přitom δ(A) = 1 : 0 ∈ A 0 : 0 /∈ A. Odtud také vyplývá její základní vlastnost ∞ −∞ f (x)(dδ(x)) = ∞ −∞ f (x)δ(x)dx = f (0). Model tedy předpokládá tenkou a nekonečně dlouhou trubici, difúze probíhá pouze ve směru osy x, přitom v čase t = 0 abstrahujeme od skutečnosti, že vstřik látky trvá nějakou dobu a probíhá do objemu, nikoliv plochy. Fakticky je ale tato abstrakce možná, protože předpokládáme velmi dlouhý časový interval pozorování difúze a dostatečně tenkou trubici oproti její délce. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Budeme tedy nyní hledat řešení u(x, t) rovnice difúze (19) ∂u ∂t = D ∂2u ∂x2 = D∇2 u Fourierovou transformací řešení u(x, t) je U(α, t) = ∞ −∞ u(x, t)e−iαx dx 60. příklad: Ukažte, že pro Fourierovu transformaci U(α, t) koncentrace u(x, t) platí diferenciální rovnice ∂U ∂t + Dα2 U = 0. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Ukažte, že pro Fourierovu transformaci U(α, t) koncentrace u(x, t) platí diferenciální rovnice ∂U ∂t + Dα2 U = 0. ∂U ∂t = ∂ ∂t ∞ −∞ u(x, t)e−iαx dx = = ∞ −∞ ∂u(x,t) ∂t e−iαx dx = ∞ −∞ D ∂2u(x,t) ∂x2 e−iαx dx = = D ∂u(x,t) ∂x e−iαx ∞ −∞ − ∞ −∞ (−iα)e−iαx ∂u(x,t) ∂x dx = = D u(x, t)(−iαe−iαx ) ∞ −∞ − ∞ −∞ α2 e−iαx u(x, t)dx , tj. ∂U ∂t + Dα2 U = 0 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Ukažte, že pro Fourierovu transformaci U(α, t) koncentrace u(x, t) platí diferenciální rovnice ∂U ∂t + Dα2 U = 0. ∂U ∂t = ∂ ∂t ∞ −∞ u(x, t)e−iαx dx = = ∞ −∞ ∂u(x,t) ∂t e−iαx dx = ∞ −∞ D ∂2u(x,t) ∂x2 e−iαx dx = = D ∂u(x,t) ∂x e−iαx ∞ −∞ − ∞ −∞ (−iα)e−iαx ∂u(x,t) ∂x dx = = D u(x, t)(−iαe−iαx ) ∞ −∞ − ∞ −∞ α2 e−iαx u(x, t)dx , tj. ∂U ∂t + Dα2 U = 0 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Ukažte, že pro Fourierovu transformaci U(α, t) koncentrace u(x, t) platí diferenciální rovnice ∂U ∂t + Dα2 U = 0. ∂U ∂t = ∂ ∂t ∞ −∞ u(x, t)e−iαx dx = = ∞ −∞ ∂u(x,t) ∂t e−iαx dx = ∞ −∞ D ∂2u(x,t) ∂x2 e−iαx dx = = D ∂u(x,t) ∂x e−iαx ∞ −∞ − ∞ −∞ (−iα)e−iαx ∂u(x,t) ∂x dx = = D u(x, t)(−iαe−iαx ) ∞ −∞ − ∞ −∞ α2 e−iαx u(x, t)dx , tj. ∂U ∂t + Dα2 U = 0 Předpokládáme-li absolutní konvergenci integrál, můžeme zaměnit pořadí integrace. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Ukažte, že pro Fourierovu transformaci U(α, t) koncentrace u(x, t) platí diferenciální rovnice ∂U ∂t + Dα2 U = 0. ∂U ∂t = ∂ ∂t ∞ −∞ u(x, t)e−iαx dx = = ∞ −∞ ∂u(x,t) ∂t e−iαx dx = ∞ −∞ D ∂2u(x,t) ∂x2 e−iαx dx = = D ∂u(x,t) ∂x e−iαx ∞ −∞ − ∞ −∞ (−iα)e−iαx ∂u(x,t) ∂x dx = = D u(x, t)(−iαe−iαx ) ∞ −∞ − ∞ −∞ α2 e−iαx u(x, t)dx , tj. ∂U ∂t + Dα2 U = 0 Dosadíme z rovnice difúze. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Ukažte, že pro Fourierovu transformaci U(α, t) koncentrace u(x, t) platí diferenciální rovnice ∂U ∂t + Dα2 U = 0. ∂U ∂t = ∂ ∂t ∞ −∞ u(x, t)e−iαx dx = = ∞ −∞ ∂u(x,t) ∂t e−iαx dx = ∞ −∞ D ∂2u(x,t) ∂x2 e−iαx dx = = D ∂u(x,t) ∂x e−iαx ∞ −∞ − ∞ −∞ (−iα)e−iαx ∂u(x,t) ∂x dx = = D u(x, t)(−iαe−iαx ) ∞ −∞ − ∞ −∞ α2 e−iαx u(x, t)dx , tj. ∂U ∂t + Dα2 U = 0 Použijeme per partes ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Ukažte, že pro Fourierovu transformaci U(α, t) koncentrace u(x, t) platí diferenciální rovnice ∂U ∂t + Dα2 U = 0. ∂U ∂t = ∂ ∂t ∞ −∞ u(x, t)e−iαx dx = = ∞ −∞ ∂u(x,t) ∂t e−iαx dx = ∞ −∞ D ∂2u(x,t) ∂x2 e−iαx dx = = D ∂u(x,t) ∂x e−iαx ∞ −∞ − ∞ −∞ (−iα)e−iαx ∂u(x,t) ∂x dx = = D u(x, t)(−iαe−iαx ) ∞ −∞ − ∞ −∞ α2 e−iαx u(x, t)dx , tj. ∂U ∂t + Dα2 U = 0 a znovu per partes spolu s okrajovými podmínkami. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Ukažte, že pro Fourierovu transformaci U(α, t) koncentrace u(x, t) platí diferenciální rovnice ∂U ∂t + Dα2 U = 0. ∂U ∂t = ∂ ∂t ∞ −∞ u(x, t)e−iαx dx = = ∞ −∞ ∂u(x,t) ∂t e−iαx dx = ∞ −∞ D ∂2u(x,t) ∂x2 e−iαx dx = = D ∂u(x,t) ∂x e−iαx ∞ −∞ − ∞ −∞ (−iα)e−iαx ∂u(x,t) ∂x dx = = D u(x, t)(−iαe−iαx ) ∞ −∞ − ∞ −∞ α2 e−iαx u(x, t)dx , tj. ∂U ∂t + Dα2 U = 0 ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Obecným řešením rovnice ∂U ∂t + Dα2 U = 0 je zřejmě funkce U(α, t) = e−Dα2t · konst. Přitom konstanta je určena počáteční podmínkou u(x, 0) = M S δ(x) a obecně může záviset na α. V našem případě však platí U(α, 0) = ∞ −∞ u(x, 0)e−iαx dx = ∞ −∞ M S δ(x)e−iαx dx = M S . Partikulárním řešením transformované rovnice je tedy U(α, t) = M S e−Dα2t . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Zpětnou transformací dostáváme u(x, t) = 1 2π ∞ −∞ U(α, t)eiαx dα = M 2πS ∞ −∞ e−Dα2t eiαx dα. Protože eiαx = cos αx + i sin αx, funkce sin je lichá a funkce cos je sudá, platí u(x, t) = M πS ∞ 0 e−Dα2t cos(αx)dα. Poznámka 14. Substituce α = ξ√ Dt , x = η √ Dt převádí výše uvedený integrál na tvar M πS √ Dt ∞ 0 e−ξ2 cos(ηξ)dξ. Označíme a vypočteme následující integrál: I(η) = ∞ 0 e−ξ2 cos(ηξ)dξ. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Nejprve si všimněme, že platí dI dη = − ∞ 0 ξ sin(ηξ)e−ξ2 dξ. Dále pak de−ξ2 dξ = −2ξe−ξ2 , tj. dI dη = 1 2 ∞ 0 sin(ηξ)d(e−ξ2 ) dξ dξ. Metodou per partes nyní můžeme přepsat tento integrál do tvaru dI dη = 1 2 sin(ηξ)e−ξ2 ∞ 0 − ∞ 0 d sin(ηξ) dξ e−ξ2 dξ . Funkce sin(ηξ)e−ξ2 je pro ξ = 0 nulová, pro ξ → ∞ také, protože sin je ohraničený a e−ξ2 → 0. Dostali jsme tedy pro I(η) diferenciální rovnici dI(η) dη = −1 2 ∞ 0 η cos(ηξ)e−ξ2 dξ = −η 2 I(η). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × dI(η) dη = −η 2 I(η) Protože I(η) = ∞ 0 e−ξ2 cos(ηξ)dξ, je počáteční podmínkou zřejmě I(0) = ∞ 0 e−ξ2 dξ. 61. příklad: Ukažte, že platí ∞ 0 e−x2 dx = √ π 2 . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Ukažte, že platí ∞ 0 e−x2 dx = √ π 2 . Místo integrálu I = ∞ 0 e−x2 dx budeme hledat jeho druhou mocninu I2 = ∞ 0 e−x2 dx · ∞ 0 e−y2 dy = Ω e−(x2+y2) dxdy Označme Ω 1. kvadrant dvojrozměrného prostoru s kartézskými souřadnicemi x, y. V přepisu do polárních souřadnic platí: Ω e−(x2+y2) dxdy = π 2 0 ∞ 0 e−r2 rdr dϕ Substitucí r2 = t dostáváme I2 = π 2 0 ∞ 0 1 2 e−t dt dϕ = π 2 0 dϕ −1 2 e−t ∞ 0 = π 4 . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Poznámka 15. Funkce erf (x) = 2√ π x 0 e−ξ2 dξ se nazývá Gaussova chybová funkce a má tvar 1 0 0,5 0 -2 -0,5 -1 -4 x 42 62. příklad: Ukažte, že řešením úlohy (19) s dříve uvedenými okrajovými a počátečními podmínkami je koncentrace u(x, t) = M S √ 4πDt e− x2 4Dt . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Analogicky lze odvodit řešení v trojrozměrném případě jako u(x, t) = M 4πt 4πDxDyDzt e − x2 4Dxt − y2 4Dyt − z2 4Dzt . Poznámka 16. V jednorozměrném případě tedy maximální koncentrace umax látky v čase klesá s 1√ t . V dvojrozměrném případě klesá s 1 t a v trojrozměrném s 1 t √ t . 63. příklad: Vytvořte animaci koncentrace u(x, t) v čase (program Maple, funkce animate). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model difúze s advekcí V této a následující kapitole model difúze rozšíříme. Uvažujme znovu rovnici zákona zachování hmoty ∂u ∂t = − ∂J ∂x + f = − ∂(Jdi f f +Jadv) ∂x + f. V předchozí kapitole studovaný model difúze předpokládal, že • Jadv = vu = 0, tj. v = 0, nedochází k advekci, tj. přenosu látky rychlostí v • f = 0, tj. v systému nevznikají nebo nezanikají žádné částice, nedochází k reakci V této kapitole porušíme první podmínku, v následující kapitole pak druhou. Ještě si povšimněme součtu Jdi f f + Jadv. Tuto superpozici můžeme provést pouze za předpokladu, že difúze a přenos média jsou navzájem nezávislé procesy. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × ∂u ∂t = − ∂(Jdi f f +Jadv) ∂x = − ∂(−D ∂u ∂x +vu) ∂x , tj. v případě, že rychlost média nezávisí na čase a místě (např. nestlačitelná kapalina proudící konstantní rychlostí) ∂u ∂t = D ∂2u ∂x2 − v ∂u ∂x . (20) Ve vícerozměrném případě pak ∂u ∂t = D n ∑ i=1 ∂2u ∂x2 i − n ∑ i=1 vi ∂u ∂xi , což můžeme zkráceně zapsat takto: ∂u ∂t = D∇2 u − v · ∇u. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Pokusme se odhadnout řešení rovnice (20). V případě, že rychlost v byla nulová, vyřešili jsme rovnici difúze a našli řešení úlohy s bodovým zdrojem jako Gaussovu křivku u(x, t) = M S √ 4πDt e− x2 4Dt . Pokud bude difundující látka unášena prostředím konstantní rychlostí v = 0, bude místo x za čas t v místě vt. Substitucí ξ = x − vt pak posuneme toto pohybující se místo do počátku. Řešením by měla být tedy koncentrace u(x, t) = M S √ 4πDt e− (x−vt)2 4Dt . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 64. příklad: Dokažte, že substituce ξ = x − vt převádí rovnici difúze s advekcí na rovnici difúze bez advekce. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model intravenózní injekce Lékař zavede injekci antihistaminika pacientovi s alergickou reakcí. Jaká bude koncentrace chemikálie v krvi za minutu? Koncepce: Stavovou proměnnou bude koncentrace antihistaminika u(x, t), parametry budou rychlost krve v, difúzní koeficient D, počáteční koncentrace antihistaminika u0 a celkový čas zavádění injekce T. Diagram: Předpokládejme, že lékař vstřikuje injekci rychlostí krve, tj. počáteční koncentrace u0 je v žíle v celé délce L = vT. Bez újmy na obecnosti můžeme ”sledovat”trasu krve od místa vpichu x = 0, tedy proměnná x bude svázána s rovnoměrně se pohybujícím tokem krve. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Intoxikovanou část žíly délky L můžeme považovat za infinitesimální součet jejích elementů, kde množství antihistaminika v elementu oblemu žíly je dM = u0Sdξ. Zřejmě platí du(x, t) = dM S √ 4πDt e− (x−ξ)2 4Dt , tj. du(x, t) = u0 √ 4πDt e− (x−ξ)2 4Dt dξ. Superpozicí na kontaminované délce L pak dostáváme následující rovnici: Rovnice: u(x, t) = L 0 u0 √ 4πDt e− (x−ξ)2 4Dt dξ. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 65. příklad: Substitucí η = x−ξ√ 4Dt převeďte integrál do tvaru u(x, t) = u0 2 erf ( x√ 4Dt ) − erf ( x−L√ 4Dt ) . 66. příklad: Nakreslete funkci koncentrace u(x, t) pro v = 0.1m/s, D = 2 · 10−5 m2 /s, T = 5s a u0 = 0.1. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Reakčně-difúzní model Jak už bylo řečeno, v této kapitole model difúze rozšíříme o reakční složku. Uvažujme rovnici zákona zachování hmoty bez advekce ∂u ∂t = − ∂J ∂x + f = − ∂Jdi f f ∂x + f. Funkce f = 0 značí, že v systému vznikají nebo zanikají částice. Typickým příkladem takového systému jsou chemické reakce v tekutinách, kde se reakcí vznikající látka difúzí šíří tekutinou. Tyto modely bývají většinou vícerozměrné a jejich řešení značně komplexní. Pro jednoduchost si ukážeme model jednorozměrný, ve kterém vzniká typický jev - postupující vlna. Už tento jednoduchý model má samozřejmě v závislosti na tvaru funkce f, okrajových a počátečních podmínkách velice komplexní chování. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Postupující vlna Uvažujme reakčně difúzní rovnici tvaru ∂u ∂t = D ∂2u ∂x2 + f (u), (21) kde se změna koncentrace u ∈ 0, 1 v čase závisí na difúzi a na reakci, která není funkcí času, pouze koncentrace, přitom f (0) = 0 a f (1) = 0. Předpokládejme nyní, že má rovnice řešení v konkrétním tvaru u(x, t) = U(x − vt) = U(ξ), kde ξ = x − vt je podobně jako u rovnice s advekcí transformace, která bude x posouvat rychlostí v ve směru osy x. Oproti rovnici s advekcí, ale nyní v není rychlostí média (ta je nulová), ale libovolně zvolenou hodnotou. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Pokud bychom tedy našli nějaké řešení U(ξ), ve skutečnosti by představovalo řešení posouvající se rychlostí v po ose x. Řešení u(x, t) odpovídající U(ξ), které ”spojuje”rovnovážné body se nazývá postupující vlnou (travelling wave). 67. příklad: Ukažte, že transformace ξ = x − vt převádí rovnici (21) do tvaru DU′′ + vU′ + f (U) = 0. Zavedením U′ = V dostáváme dvojrozměrný systém U′ = V V′ = − v D V − f (U) D . (22) Stacionárními body pak budou zřejmě nulové body funkce f, které jsou minimálně dva: U = 0 a U = 1, jeden z nich bývá sedlem a separatrix sedla tyto rovnováhy ”spojuje”. Můžeme tak nalézt posupující vlnu u(x, t). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Typickým příkladem takovéto rovnice je Fisherova-Kolmogorovova rovnice, kde reakční funkce je tvaru f (u) = ru(1 − u). 68. příklad: Tuhle funkci už jsme někde měli, že? Co by asi mohla modelovat Fisherova-Kolmogorovova rovnice? Co by představovala proměnná u? 69. příklad: Pro Fisherovu-Kolmogorovovu rovnici proveďte předchozí transformace na dvojrozměrný systém (22), najděte jeho stacionární body a určete jejich typ. Nakreslete fázový portrét a simulujte v Xppautu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Jiným příkladem reakčně-difúzní rovnice je FitzHughova-Nagumova rovnice, kde reakční funkce je tvaru f (u) = ku(u − a)(1 − u), která popisuje model šíření vzruchu v neuronu. Proměnná u představuje normalizované membránové napětí neuronu (membránové napětí představuje rozdíl potenciálů uvnitř a vně neuronu, je určeno koncentrací iontů). 70. příklad: Ukažte, že FitzHughova-Nagumova rovnice má stabilní i nestabilní řešení, u = a, u = 0 a u = 1. Nakreslete řešení s různými počátečními podmínkami a pokuste se je interpretovat vzhledem k času a prostoru. V případě u = 0 mluvíme o nervu v depolarizovaném stavu, v případě u = 1 mluvíme o nervu v polarizovaném stavu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model šíření epidemie typu SIR Koncepce: Uvažujme epidemiologický model SIR se smrtelnou chorobou: S′ = −βSI I′ = βSI − νI R′ = νI, (23) kde β, ν > 0 jsou parametry a stavové proměnné reprezentující okamžitý počet S náchylných, I infekčních a R uhynulých jedinců v čase. Uvažujme ovšem situaci, kdy se nakažení jedinci přemisťují v prostoru. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Rovnice: ∂S ∂τ = −βSI ∂I ∂τ = βSI − νI + D∇2 I ∂R ∂τ = νI, (24) kde člen s D > 0 odpovídá difúzi (šíření) nemocných do prostoru. Pro jednoduchost budeme předpokládat, že prostorová proměnná je pouze jednorozměrná (ξ) a že se pohybují pouze nemocní jedinci. Může jít např. o model vztekliny u lišek v kaňonu :-), které se přemisťují v důsledku konfliktů ve svém původním teritoriu. Předpokládáme, že po skončení epidemie bude populace lišek prostorově homogenní, ale zmenšená o uhynulé jedince, tj. bude platit lim τ→∞ S(ξ, τ) = S∞ < N, lim τ→∞ I(ξ, τ) = 0, lim τ→∞ R(ξ, τ) = N − S∞ < N. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 71. příklad: Ukažte, že substituce s = S N , i = I N , r = R N , t = ντ, x = ν D ξ, R0 = βN ν převádí předchozí systém na ∂s ∂t = −R0si, ∂i ∂t = R0si − i + ∂2i ∂x2 , ∂r ∂t = i. (25) 72. příklad: Hledejte postupující vlnu, tedy zaveďte novou proměnnou z = x − vt a přejděte k diferenciálnímu systému s proměnnými U(z) = s(x, t), V(z) = i(x, t), W(z) = r(x, t). Vysvětlete proč hledáme řešení splňující lim z→∞ U(z) = 1, lim z→−∞ U(z) = S∞ N , lim z→±∞ V(z) = 0, lim z→±∞ V′ (z) = 0, lim z→∞ W(z) = 0, lim z→−∞ W(z) = N−S∞ N . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Systém je tvaru U′ = R0 v UV, V′ = −R0 v UV + 1 v V − 1 v V′′ W′ = −1 v V. (26) Přitom V′′ = dV′ dU dU dz = dV′ dU R0 v UV, tedy dV dU = −1 + 1 R0U − 1 v dV′ dU . Integrací podle U pak V = −U + 1 R0 ln U − 1 v V′ + c, kde pro z → ∞ platí 0 = −1 + 0 − 0 + c. Dostáváme tak systém U′ = R0 v UV, V′ = v(−U − V + 1 R0 ln U + 1). (27) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Z podmínek pro z → −∞ pak 0 = v(1 − S∞ N + 1 R0 ln S∞ N ). Odtud R0 = ln S∞ N S∞ N − 1 > 1 73. příklad: Nakreslete fázový portrét. Ukažte, že pro rychlost šíření epidemie musí platit v ≥ 2 R0 − 1. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Model šíření genu v populaci Koncepce: Uvažujme populaci, ve které se šíří mutace genu — alela a namísto původní alely A. Frekvence alely a je p, frekvence alely A je q = 1 − p. Genová mutace se navíc šíří do prostředí náhodně migrací (náhodná procházka, difúze). Budeme uvažovat pro jednoduchost pouze jednu prostorovou proměnnou a aditivní selekční koeficient genu, tedy půjde o případ, kdy gen A ani a není dominantní, ale míra selekce genotypů AA, Aa a aa aditivně roste, tj. relativní fitness je 1, 1 + s a 1 + 2s. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Fitness genu v populační genetice Pravděpodobnost, že nějaký fenotyp přežije a zanechá potomky je mírou jeho fitness. Budeme předpokládat, že fitness genotypů je konstantní a je rovna pravděpodobnosti jeho přežití. Mluvíme o tzv. absolutním fitness, protože jeho hodnota je závislá na fitness ostatních genotypů. Obvykle však známe hodnotu životaschopnosti každého genotypu vztaženého relativně k ostatním vybraným genotypům jako standard k porovnání. Relativní fitness w tedy vyjadřuje podíl potomků produkovaných jedním genotypem v porovnání s genotypem jiným, jakousi reprodukční způsobilost genotypu. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: A a A a AA Aa aA aa Rovnice: genotyp aa aA AA původní frekvence p2 2pq q2 podíl po selekci p2 (1 + 2s) 2pq(1 + s) q2 frekvence genotypu po selekci p2(1+2s) w 2pq(1+s) w q2 w kde w = p2 (1 + 2s) + 2pq(1 + s) + q2 = 1 + 2sp. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × V následující generaci bude tedy frekvence mutace genu a rovna p2(1 + 2s) + pq(1 + s) w . Změna (derivace) bude tedy odpovídat rozdílu frekvencí alely genu a ˙p = p2(1 + 2s) + pq(1 + s) w − p = p2(1 + 2s) + pq(1 + s) − p(1 + 2sp) 1 + 2sp = sp(1 − p) 1 + 2sp . Bude-li p(x, t) představovat frekvenci alely a na daném místě v daném čase a mutace genu se bude difúzně šířit do okolí (náhodná procházka), pak musí platit rovnice ˙p = D ∂2p ∂x2 + sp(1 − p) 1 + 2sp , (28) ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 74. příklad: Hledejte řešení jako postupující vlnu a převeďte rovnici na systém ODR. Najděte rovnovážné body rovnice a určete jejich stabilitu. Nakreslete nulkliny a fázový diagram systému. Vyhodnocení: Z fázového diagramu je zřejmé, že frekvence alely P(ξ) = P(x − vt) = p(x, t) klesá k nulové stabilní rovnováze pro ξ → ∞. Vzhledem k času jde ale o postupující vlnu šíření genu, x v · t P(x − vt) což je zcela v souladu s očekáváním. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Šíření kolonií mikroorganismů Koncepce: Kvasinky jsou jednobuněčné houbové mikroorganismy, množí se zejména nepohlavně a je pro ně charakteristický způsob dělení buněk, takzvané pučení. Buňky kvasinek potřebují ke svému dělení energii, kterou získávají z cukru - glukózy. Pokusíme se vytvořit model šíření kolonie kvasinek. Pro jednoduchost budeme uvažovat jen jednu prostorovou proměnnou x. Počet buněk v jednotce objemu (hustotu buněk) označíme n(x, t), koncentraci glukózy g(x, t). Glukóza se ve vodě šíří difúzí. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Diagram: Rovnice: ˙n = kn(g − g∗ ) ˙g = D ∂2g ∂x2 − ckn(g − g∗ ). ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 75. příklad: Zdůvodněte uvedený tvar rovnic. Zavedením nových proměnných z = x − vt, N(z) = n(x, t), G(z) = g(x, t) − g∗ dostáváme −v dN dz = kNG −v dG dz = D d2G dz2 − ckNG. Přičtením c-násobku 1. rovnice k druhé pak dostáváme −vc dN dz − v dG dz = D d2G dz2 . Integrací pak dostáváme konst − vcN − vG = D dG dz . ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Konstantu můžeme dopočítat např. touto úvahou: v případě, že G ≡ 0, je také dG dz ≡ 0, proto konst = vcN0 pro G(z) ≡ 0. N0 tedy představuje jakési hraniční maximální množství (hustota) kvasinek, když je glukóza již vypotřebována. Dostáváme tedy systém rovnic dN dz = − k v NG dG dz = vcN0 D − vc D N − v D G. 76. příklad: Najděte rovnovážné body rovnice a určete jejich stabilitu. Nakreslete nulkliny a fázový diagram systému. Nakreslete graf řešení N versus z a G versus z ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Vyhodnocení: Z fázového diagramu je zřejmé, stabilní rovnováhou je bod [0, cN0], přitom hustota buněk řešení klesá se z → ∞ k nule a koncentrace glukózy roste k cN0. Vzhledem k času jde o postupující vlnu šířící se kolonie kvasinek, která vypotřebovává glukózu, což je zcela v souladu s očekáváním. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × 77. příklad: Podívejte se na článek Gray-Scottův model a jeho simulaci Postupující vlny a vznik vzorů k nakouknutí: Nerovnovážná termodynamika a její aplikace, ZČU v Plzni ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × Reference [1] Nicholas Britton. Essential mathematical biology. Springer Science & Business Media, 2012. [2] AC Chiang and Kevin Wainwright. Fundamental methods of mathematical economics. 2005. [3] Leah Edelstein-Keshet. Mathematical models in biology. Siam, 1988. [4] Stephen P Ellner and John Guckenheimer. Dynamic models in biology. Princeton University Press, 2011. [5] G Bard Ermentrout and David H Terman. Mathematical foundations of neuroscience, volume 35. Springer Science & Business Media, 2010. [6] Y. A. Kuznetsov. Elements of applied bifucation theory, Second edition, Applied Mathematical Sciences 112. Springer-Verlag, Berlin, Heidelberg, New York, 1998. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × [7] Stephen Lynch. Dynamical systems with applications using MapleTM. Springer Science & Business Media, 2009. [8] James D Murray. Mathematical Biology II: Spatial Models and Biomedical Applications, vol. 18 of Interdisciplinary Applied Mathematics. Springer-Verlag New York Incorporated, 2001. [9] James D Murray. Mathematical Biology I: An Introduction, vol. 17 of Interdisciplinary Applied Mathematics. Springer, New York, NY, USA, 2002. [10] Paul E Phillipson and Peter Schuster. Modeling by nonlinear differential equations. World Scientific, 2009. [11] Pierre NV Tu. Dynamical systems: An introduction with applications in economics and biology. Springer Science & Business Media, 2012. [12] Edward K Yeargers, James V Herod, and Ronald W Shonkweiler. An introduction to the mathematics of biology: with computer algebra models. Springer Science & Business Media, 2013. ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita × KONEC ⊳⊳ ⊳ ⊲ ⊲⊲ c 2019 Masarykova univerzita ×