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).

Matematické modely reálných situací

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.

Trocha historie

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.

• Někdy se sestavovaly modely stejného procesu v menším měřítku. Podívejte se na příběh Vteřiny před katastrofou - Tsunami v horách. Seznamte se s tím, co se tehdy odehrálo a poté dokoukejte (nebo přetočte) na čas 39 minut. Na zmenšeném modelu údolí Vajont tam shazují kameny jako simulaci pádu hory do nádrže. Klíčová slova jsou "Tým inženýrů zasvěcenně odhadne, že jedna minuta je nejrychlejší možnou variantou skutečného sesuvu. V přepočtu pro simulaci se jedná o zhruba 4 sekundy." Skutečně, ve zmenšeném modelu plyne čas jinou rychlostí a k odhalení tohoto efektu nám stačí znalost diferenciální rovnice a toho, jak se tato rovnice mění při změnách jednotek. Nemusíme rovnici vůbec řešit!
• Nekdy se podařilo vybudovat systémy, ve kterých se tok tekutiny řídí stejnými rovnicemi, jako jsou rovnice studovaného modelu. To je podstatou hydromechanických počítačů. Například počítač pro model ekonomiky.
• Totéž čistě mechanicky vedlo ke konstrukci mechanických počítačů. U jejich sestavení stáli Vanevar Bush (duchovní otec Internetu) a Antonín Svoboda (významný český matematik, o jeho životě viz krátký seriál na root.cz: díl 1 o počátku kariéry, díl 2 o válečných úspěších s mechanickými počítači, díl 3, díl 4 atd). První mechanické počítače byly sestavovány ze stavebnice, viz obrázek.

MENDELU učí distančně, tak pojďme na Harvard

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.

• První video Spruce Budworm Outbreaks nás seznámí s problémem lesního škůdce, který dokáže totálně zničit lesy.
• Další dvě videa, 5.1 - Our Dynamic Earth and Malthus' Population Model a 5.2 - Limits on Growth - Verhulst's Model, ukazují, jak se dá obecně růst populace modelovat pomocí diferenciálních rovnic a jak je možné modelovat růst populace v prostředí s omezenou nosnou kapacitou. Je zde zdůvodněno, proč je $$\frac{\mathrm dP}{\mathrm dt}=rP\left(1-\frac PK\right)$$ rozumným matematickým modelem pro růst populace.
• Od videa Change and Bifurcation se soustředíme na to, že změny ve vstupních datech způsobují změny v chování modelu a někdy tato změna může být dramatická. Někdy i mírné zvýšení lovu způsobí překročení meze, za kterou již nastává nevyhnutelně vymření populace. Takovým stavů se říká bifurkace. Ukázky jsou ve videu More Examples of Bifurcations - Fishing and Weather.
• Ve videích od 7.2 Modelling Fish Populations přes 7.3 - How Does Fishing Affect the Population až po 7.4 - Mathematics and Biology - A Symbiotic Relationship se zabýváme lovem populace a jeho udržitelností. Je zde na modelu obsahujícím derivaci (ale jinak středoškolskými metodami) rozebráno, jaký má intenzita lovu $\alpha$ v rovnici $$\frac{\mathrm dP}{\mathrm dt}=rP\left(1-\frac PK\right)-\alpha$$ vliv na to, jestli populace dlouhodobě přežívá, či nikoliv.
• Ve videích dále se věnujeme růstu populace obaleče. Jedná se vlastně o nekonstantní lov podle modelu $$\frac{\mathrm dP}{\mathrm dt}=rP\left(1-\frac PK\right)-\frac{P^2}{P^2+1}.$$ Tento lov zprostředkují ptáci. Ti při malých populacích obaleče loví jinou, lépe dostupnou potravu. Situace se změní, pokud je obaleče dostatek. Poté jej ptáci ochotně využijí jako hlavní potravu svého jídelníčku. Tím pomáhají při stabilizaci systému. Situace se dále mění při přemnožení obaleče nad hladinu, kdy jsou predátoři saturováni. Ptáci nezvládnou sezobat celý švédský stůl a populace obaleče, která byla dosud držena na uzdě, se rozmnoží až na nosnou kapacitu prostředí. To ovšem vede k tomu, že jsou napadeny všechny stromy a dojde k destrukci lesa. Ve videu 8.1 - A Population Model for Budworms je zopakován základní logistický model. Začlenění vlivu predátorů je diskutováno ve videu 8.2 - A Budworm Population Model with Predation. Zejména si všimněte, že predátoři nejsou nenasytní a jejich vliv je shora ohraničený. Proto musíme pro modelování vlivu predátorů hledat mezi funkcemi, které nerostou neohraničeně. Prolematika stabilních bodů systému při rostoucí nosné kapacitě prostředí (tj. zhruba v rostoucím lese) je obsahem videa 8.3 - Equilibrium Points and Carrying Capacity. A na závěr video 8.4 - Increasing the Carrying Capacity Further, které popisuje mechanismus, jak relativně malé rozkolísání může překopit populaci obaleče ze stabilního stavu s nízkou velikostí populace do stavu s mnohonásobně vyšší velikostí.

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í.

Matematický model obecně

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 casu

Tato 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.

Matematický model konstantního lovu

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.)

• Pro populaci bez lovu ($\alpha=0$) populace doroste do nosné kapacity prostředí i pro minimální počáteční hodnotu populace.
• S přednastaveným lovem $\alpha=500$ populace vymře při malé počáteční hodnotě velikosti populace a je udržitelně lovena při velké počáteční hodnotě populace. Chování modelu se mění někde okolo hodnoty počáteční velikosti populace 5800. Abyste viděli nejenom pokles, ale stav populace až do okamžiku vyhynutí, nastavte proměnnou konec na hodnotu 50.
• Pro zvýšenou intenzitu lovu existuje hranice, kdy populace nemůže za žádných okolností dlouhodobě přežívat. Křivka vlevo se celá dostane pod osu $x$. Pro dlouhodobé sledování vývoje pomalu se měnící populace možná budete muset prodloužit simulaci. Například pro $\alpha=1020$ můžete nastavit $konec=500$.

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.

Matematický model přemnožení obaleče

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.

Matematický model přemnožení kůrovce

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).

Zadání domácí samostatné práce

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.

1. Prozkoumejte model $$\frac{\mathrm dP}{\mathrm dt}=rP\left(1-\frac PK\right)-\alpha$$ z hlediska změn nosné kapacity prostředí. Vyjděte z hodnot v našem základním nastavení, zvyšte lov na $\alpha = 950$ a zjistěte, jestli populace vymře či nikoliv. Poté ponechte všechny parametry stejné, i parametr $\alpha$ ponechejte na hodnotě $950$ a měňte pouze nosnou kapacitu prostředí $K$. Sledujte, jak se mění chování modelu a zda se může přežívání změnit ve vyhynutí nebo naopak. Pravděpodobně budete muset pro sledování vývoje během delšího časového úseku navýšit hodnotu konec z předvolené hodnoty 20 nejméně na desetinásobek. Sepište své pozorování do relativně uceleného celku: jakou rovnici používáme, co tato rovnice modeluje, jaké máme přednastavené hodnoty, jak se při těchto hodnotách chová model z hlediska trvalé udržitelnosti. Dále napište výsledek svého experimentování, tj. jak se na této trvalé udržitelnosti projeví změna nosné kapacita prostředí, zda je matematicky potřeba pro změnu chování nosnou kapacitu zvětšovat či zmenšovat, kde je přibližně kritická hodnota, která mění chování systému a zda pozorování souhlasí s očekáváními, která máme z biologického úhlu pohledu. Obrázky vložte do textu (pravé tlačítko na obrázek a buď Copy image nebo Save image as), pokud píšete text na počítači. Obrázek charakterizující vyhynutí by mohl vypadat například tak, jak je na ukázáno v tomto odstavci.
2. Po splnění úkolu v modelu obaleče v Kapitole 5 jste získali dvě sady koeficientů $r$ a $q$. Koeficient $q$ je stejný a koeficient $r$ odlišuje stavy s jedním průsečíkem (jediný stabilní stacionární stav) a třemi průsečíky (jediný nestabilní a dva stabilní stacionární stavy) křivek. Použijte tyto parametry pro simulaci vývoje populace z malé hodnoty, například $0.01$. Tak je zajištěno, že se populace ustálí na nejmenším stabilním stavu. Simulace proveďte Eulerovou metodou. Snadný postup je modifikovat model lovu s logistickým růstem. Nastavte dostatečně dlouhou dobu simulace tak, aby šlo zřetelně vidět, že v jednom případě dojde k dramatickému nárůstu a ve druhém zůstane populace stabilizovaná na relativně malé hodnotě. Výstup bude obsahovat: formulaci modelu, slovní popis co tento model popisuje, specifikaci počáteční podmínky, nastavení parametrů $r$ a $q$ pro stav, kdy nedojde k přemnožení, příslušný obrázek, nastavení kdy dojde k přemnožení, příslušný obrázek. Srovnání obou křivek by mohlo vypadat tak jak je na obrázku připojeném k této části. Červená křivka odpovídá rozšíření populace téměř k nosné kapacitě prostředí, která je v modelu nastavena na hodnotu 18. Černá křivka odpovídá stejné nosné kapacitě, ale hodnota $r$ je dostatečně malá na to, aby predátoři populaci udrželi na uzdě. (Výstup je z programu XPPAUT, který je přímo stvořen pro simulace a práce s ním je extrémně rychlá, ovšem až po zapracování a získání zkušeností). Nemusíte však kreslit obě křivky do jednoho obrázku.

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.