Na konci textu je zadání práce, za kterou je možno získat 4 body do celkového hodnocení. Bodové hodnocení je poměrně luxusní, nevyjadřuje však obtížnost úkolu. Jeho cílem je motivovat. Úkol sám o sobě težký není, spočívá ve změnách konstant v předpřipraveném skriptu (první část) a ve změně modelu a parametrů (druhá část).
Co trápí svět dneška? Určitě koronavirus, který jsme modelovali zde a také devastace životního prostředí při přemnožení škůdců. V současnosti Českou republiku trápí mimo jiné třeba kůrovec. Ale i jiné země mají podobné problémy.
Minulý týden jsme se seznámili s tím, jak je možno pomocí derivací a rovnic obsahujících derivace modelovat reálné děje. Nyní na to navážeme. Ukážeme si, jak modelovat růst populací, vyzkoušíme si to na modelu trvale udržitelného rybolovu a poté se podíváme na model obaleče v kanadských lesích a ve stručnosti i na model kůrovce v českých lesích.
Mnoho rovnic z minulého týdne je podobných nebo stejných a to umožňuje přenášet znalosti z jednoho děje na jiný. Pokud totiž děje probíhají stejně, nemusí nás zajímat, že zákonitosti jsou jiné. Že jednou jde o fyzikální proces, jednou o chemickou reakci a jindy o tok peněz nebo mechanický pohyb. To umožnilo některé zajímavé počiny.
Kvalitní série videí z Harvardovy univerzity. Pokud máte problémy s mluvenou angličtinou, zapněte si titulky. Pokud máte problém i s titulky, stáhněte si je v txt souboru a protáhněte je přes Google překladač.
Úkol. Shlédněte videa z tohoto seznamu videí. Snažte se jim maximálně porozumět. Matematicky tam není nic náročného, jde hlavně o vysvětlování, povídání, principy. V seznamu je 12 videí, jen jedno přesáhne délkou 5 minut.
V následujícím textu se k videím budeme vracet a budeme používat myšlenky z nich. Nejprve však obecný přístup k numerickému modelování chování rovnice s derivací.
Vývoj systému popsaného matematickým modelem $$\frac{\mathrm dy}{\mathrm dt}=f(y)\quad y(0)=y_0$$ můžeme modelovat tak, že okamžitou rychlost nahradíme průměrnou rychlostí $$\frac{\Delta y}{\Delta t}$$ za čas $\Delta t$. Model má potom tvar $$\frac{\Delta y}{\Delta t}=f(y)$$ a změnu veličiny $y$ můžeme vyjádřit vztahem $$\Delta y=f(y)\Delta t.$$ Veličina $\Delta y$ je přírůstek veličiny $y$, veličina $\Delta t$ je přírůstek času. Jedná se vlastně o krok, po kterém se posunujeme v čase. Pseudokód pro aproximaci modelu by potom vypadal následovně.
t=0, y=y0 # inicializace krok=0.1 # nastavení kroku opakuj dy=f(y)*krok # vypocet zmeny veliciny y y=y+dy # vypocet nove hodnoty veliciny t=t+krok # posun casuTato metoda se nazývá Eulerova metoda. Na rovinu si řekněme, že není moc dobrá, je primitivní a používá se jenom z pedagogických důvodů. Pro řešení diferenciálních rovnic existuje řada robustních metod předprogramovaných specialisty na numerickou matematiku. Tyto metody jsou okamžitě k použití v programech, ve kterých se očekává modelování pomocí diferenciálních rovnic a jsou rychlé a přesné. Pro naše primitivní pokusy však postačí i Eulerova metoda.
V modelech s konstantním lovem (7.3 - How Does Fishing Affect the Population) byl studován vliv parametru $\alpha$ na to, do jakého stavu bude model $$\frac{\mathrm dP}{\mathrm dt}=rP\left(1-\frac PK\right)-\alpha$$ konvergovat. V programu Sage to může vypadat například takto. Vyzkoušejte si následující modifikace modelu. (Model otevřte v novém panelu, poeditujte dle potřeby a klikněte na tlačítko pro odeslání výpočtu.)
Program Sage je doporučený, protože používá vyspělý a snadno použitelný skriptovací jazyk Python. Pro mnoho úkolů je skriptování lepší než klikání v grafickém prostředí a matematické simulace jsou typickým představitelem takového úkolu. Můžete si ale výpočty nasimulovat v libovolném tabulkovém procesoru (Excel, Calc) nebo v libovolném jiném programu, například R, XPPAUT nebo Maxima.
Ve videu od "8.1. A population model for budworms" výzkumníci z Harvardu vysvětlují, proč je rovnice $$\frac{\mathrm dP}{\mathrm dt}=rP\left(1-\frac Pq\right)-\frac{P^2}{P^2+1}$$ dobrým modelem pro populaci obaleče. Tento model je již zjednodušený, ve skutečnosti jsou i ve druhém členu parametry udávající míru působení predátorů na obaleče a člen má tvar například $\frac{\alpha P^2}{P^2+\beta}$. Analýza je ve videu 8.2. vedena tak, že je rovnice upravena do tvaru $$\frac{\mathrm dP}{\mathrm dt}=P\left[r\left(1-\frac Pq\right)-\frac{P}{P^2+1}\right].$$ Znaménko výrazu na pravé straně sledujeme tak, že porovnáváme graf funkce $y(P)=\frac{P}{P^2+1}$ s grafem klesající přímky $y(P)=r\left(1-\frac Pq\right)$ pro $r=\frac 12$ a proměnlivý parametr $q$. Viz obrázek. Pokud je modrá přímka nad grafem červené funkce, závorka je kladná, derivace populace obaleče podle času je kladná a velikost populace roste. V opačném případě velikost populace klesá.
Ve videích 8.3 a 8.4 je ukázáno, že křivky mohou mít jeden průsečík nebo tři průsečíky (hraniční případ se dvěma průsečíky je výjimečný a pravděpodobnost jeho realizace je nulová). V prvním případě je průsečík pro relativně malé hodnoty a tomu odpovídá relativně malá velikost populace obaleče. Ve druhém případě existují dva stabilní stavy, druhý z nich je násobně větší než první a odpovídá dramatickému zvýšení populace obaleče. Oba stavy se liší nosnou kapacitou prostředí $q$ a ta souvisí s velikostí lesa.
Proces tedy funguje tak, že v mladém lese je obaleč v malé míře. Společné působení nosné kapacity prostředí a ptáků, působících jako predátoři, mu nedovolí se přemnožit. Jak les roste, roste i nosná kapacita prostředí a objeví se možnost dalšího stabilního stavu (v obrázku pravý průsečík křivek). Obaleč však je stále na nižší hladině (v obrázku levý průsečík křivek). Na uzdě jej drží predátoři. Stačí však malý výkyv, aby populace přeskočila nad nestabilní stav (v obrázku prostřední průsečík křivek) a škůdce se přemnoží. Ptáci jsou saturovaní a nestihnou populaci obaleče zredukovat.
Ve videu "8.4. Increasing the Carrying Capacity Further" je pevně zvolena konstanta $r=\frac 12$. To s sebou nese, že analýza nedokáže zachytit ještě jeden jev: pokud přímku posuneme poněkud nahoru, dolní stabilní stav zmizí a model má pouze jeden stabilní stav s obrovskou velikostí populace. To bude úkol pro nás.
Úkol. Nakreslete funkce $y(P)=\frac{P}{P^2+1}$ a $y(P)=r\left(1-\frac Pq\right)$ do jednoho obrázku a pokuste se najít nějakou hodnoty $r$ a $q$, pro které jsou tři průsečíky (tak jako na obrázku) a poté hodnoty, pro které je jeden průsečík, ale v oblasti klesání červené křivky. Můžete vyjít z předpřipraveného příkladu. Tento příklad obsahuje křivky s jedním průsečíkem v rostoucí části červené křivky pro $r=\frac 12$. Měňte postupně $q$, dokud nedostanete průsečíky tři. Poté měňte postupně $r$, dokud dva levé průsečíky nezmizí a nezůstane pouze ten v klesající části červené křivky.
Zvýšení parametru $r$ odpovídá rychlejší reprodukci obaleče. Například pokud se za stejnou dobu vylíhnou tři populace namísto dvou (jako jsme viděli u kůrovce v létě 2018), je parametr $r$ o polovinu větší. Dále se parametr může navýšit oproti minulým letům vlivem velmi mírného průběhu zimy, kdy nedojde k takové redukci populace, na jakou jsme byli v minulosti zvyklí. Původní článek s popsaným modelem D. Ludwig, D. D. Jones and C. S. Holling, Qualitative Analysis of Insect Outbreak Systems: The Spruce Budworm and Forest Journal of Animal Ecology 1978, pp. 315-332 vysvětluje nárůst tohoto parametru tak, že les prospívá a umožňuje rychlejší rozmnožování obaleče.
Popsaný model sestavili matematikové Ludwig a Jones s ekologem Hollingem, aby nahradili starý počítačový model Canadian Forestry Service, který používal 30 654 proměnných spojených diferenčními rovnicemi. Nový model pracoval jenom se třemi rovnicemi. (Zde jsme představili zjednodušenou verzi o jediné rovnici. Model obsahuje ještě další dvě rovnice, které ovliňují nosnou kapacitu prostředí.) Model kromě dobrých výsledků podal i objasnění příčiny náhlého masového přemnožení obaleče, ke kterému docházelo přibližně vždy po 35 letech. Původní počítačový model dokázal pouze počítat simulace.
Dynamika kůrovce je jiná než u obaleče. U obaleče hraje roli rostoucí nosná kapacita prostředí a predace ptáky. Kůrovec nemá významné predátory a jeho přemnožení s nosnou kapacitou tolik nesouvisí. Matematický model, jehož autorem je přední český odborník na matematickou biologii, prof. Vlastimil Křivan, je publikován v Krivan, V., Lewis, M., Bentz, B., Bewick, S., Lenhart, S., Liebhold, A. 2016. A dynamical model for bark beetle outbreaks. Journal of Theoretical Biology 407:25-37. 10.1016/j.jtbi.2016.07.009. Nahlédnutím do modelu zjistíte, že se jedná o něco jako SIR model, který jsme studovali v minulém týdnu. Opět jsou populace, které jsou ve středu našeho zájmu, rozděleny do několika skupin. Pro každou z těchto skupin jsou definovány procesy, které tuto skupinu obohacují a procesy, které tuto skupinu ochuzují. Výsledkem je soustava rovnic, kterou je možno studovat po matematické stránce tak, jak je uvedeno v článku, nebo je možno ji prozkoumávat numericky tak, jak jsme dělali s SIR modelem. Modely založené na myšlence rozdělení populace na skupiny a sledování rychlosti změn těchto skupin se nazývají kompartmentové modely a jsou velice silným nástrojem pro simulace. Bohužel rostoucí počet skupin dramaticky zvětšuje složitost a odkazuje nás v podstatě pouze právě na numerické simulace. Proto jsou v modelu kůrovce použity triky, které umožní soustavu čtyř rovnic (1) zredukovat na soustavu dvou rovnic (10), kterou je již možné prozkoumat metodami založenými na vlastnostech rovinných křivek (v publikaci o kůrovcovi celostránkový obrázek na straně 30).
Pro předmět Matematika a Aplikovaná matematika. Dotace za splnění jsou 4 body do celkového hodnocení ke zkoušce při odevzdání do 17.3.2020, 9:00 CET.
Text musí obsahovat požadované položky, stačí však velmi stručně (odpovědi několika slovy, jedna krátká věta může zahrnovat odpověď na více úkolů). Musí jít o jediné PDF, které je buď výstupem textového procesoru nebo jde o kvalitně oskenovaný text napsaný úhledně rukou a s grafy načrtnutými od ruky. V případě skenu prosím o vyváženost mezi dobrou čitelností a velikostí výsledného souboru. Nejsou přípustné seskládané fotky z mobilu ani jiné obtížně čitelné podklady. Předpokládaný rozsah při zpracování v textovém editoru: dvě strany A4, každý model přibližně na jednu stranu. Záleži však na velikosti obrázků.
UPDATE 7.3.2020: Doporučená strukutra odevzdávaného textu je zde
UPDATE 14.3.2020: Kdo se trápí s programem Sage, pro toho jsem nahrál video s ukázkami práce přesně na tom, co máme plnit v tomto úkolu. V prvním modelu jsem zvyšoval lov, aby populace vymřela, ale mělo se manipulovat s nosnou kapacitou K. Omlouvám se, alespoň máte vzor jak měnit parametr a místo lovu si měňte nosnou kapacitu. Měli byste odhadnout z logiky věci, jestli zvyšovat nebo snižovat a potvrdit si to spuštěním výpočtu.
UPDATE 17.3.: Vystavil jsem moje řešení opatřené ještě krátkou shrnující poznámkou. Zrevidujte si prosím svá řešení, upravte případné chyby, vyřešte případné nesrovnalosti a přineste ke zkoušce.