Obsah Prolog 3 Fibonacciovi králíci a jejich modifikace . . . . . . . . . . . . . . . . . . . . . . . . . 3 Leslieho model růstu populace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Vektory, matice a operace s nimi . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1 Konstrukce modelů 21 1.1 Stavové proměnné . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.2 Modely disperse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2 Modely s konstantní projekční maticí 33 2.1 Příklad — populace strukturovaná podle plodnosti . . . . . . . . . . . . . . . 33 2.2 Perronova-Frobeniova teorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.3 Řešení projekční rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.4 Transientní dynamika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.5 Analýza citlivosti a pružnosti . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.6 Analýza věkově strukturované populace . . . . . . . . . . . . . . . . . . . . . 61 2.7 Události v životním cyklu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3 Identifikace parametrů modelu 79 3.1 Inversní metody časových řad . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 3.2 Parametry populace se stabilizovanou věkovou strukturou . . . . . . . . . . . 85 4 Modely s externí variabilitou 89 4.1 Sezónní variabilita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.2 Periodická variabilita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.3 Aperiodická variabilita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5 Modely s interní variabilitou 99 5.1 Příklad — populace strukturovaná podle plodnosti . . . . . . . . . . . . . . . 99 5.2 Konstrukce modelů . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.3 Asymptotické vlastnosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6 Modely dvojpohlavní populace 115 6.1 Populace strukturovaná podle plodnosti . . . . . . . . . . . . . . . . . . . . . 115 6.2 Věkově strukturovaná dvojpohlavní populace . . . . . . . . . . . . . . . . . . 118 1 2 OBSAH Prolog Fibonacciovi králíci a jejich modifikace Leonardo Pisánský, známější jako Fibonacci, se narodil kolem roku 1170 v italské Pise a zemřel roku 1250. Vzdělání získal v severní Africe, kde jeho otec Guilielmo Bonacci působil jako diplomat. Svoje vědomosti sepsal do knihy Liber abaci. Toto dílo publikované roku 1202 má hlavní zásluhu na tom, že v Evropě byl přijat poziční systém zápisu čísel (pomocí indických symbolů, kterým dnes říkáme arabské číslice). Ve třetí části knihy Fibonacci zformuloval a řešil úlohu: Kdosi umístil pár králíků na určitém místě, se všech stran ohrazeném zdí, aby poznal, kolik párů králíků se při tom zrodí průběhem roku, jestliže u králíků je tomu tak, že pár králíků přivede na svět měsíčně jeden pár a že králíci počínají rodit ve dvou měsících svého věku.1 Tuto úlohu a její řešení lze považovat za jeden z prvních matematických modelů růstu populace. Budeme ji řešit s použitím současné symboliky. Ze zadání úlohy plyne, že králíky můžeme rozdělit do dvou kategorií (tříd) — na ty, kteří jsou mladší než dva měsíce a tedy dosud „nerodí potomky, a na ty staré aspoň dva měsíce a tedy plodné. Označme x(t), resp. y(t), počet párů juvenilních (mladých, dosud neplodných), resp. dospělých (plodných), králíků v t-tém měsíci. Z poněkud vágního Fibonacciova popisu však není jasné, co přesně má vyjadřovat „počet párů králíků v t-tém měsíci . Budeme si tedy představovat, že každý měsíc v určený den proběhne sčítání králíků, kterým získáme hodnoty x(t) a y(t). Nyní je potřeba vyjasnit, kdy se nové páry rodí. Jedna z možností je, že také k porodům dochází určitý den v měsíci. Abychom úvahy dále zjednodušili (a zreprodukovali Fibonacciův výsledek) budeme předpokládat, že králíci se rodí první den a jejich sčítání provádíme poslední den měsíce. Při sčítání mají tedy novorození králíci věk již jeden měsíc. Při sčítání následujícího měsíce mají tito králíci již věk dva měsíce a patří tedy mezi plodné. Poněvadž pár plodných králíků „zrodí (tj. zplodí a porodí) jeden pár mladých, bude počet párů mladých v t-tém měsíci stejný jako počet párů plodných v měsíci předchozím, x(t) = y(t − 1). (1) Králíci jsou na místě ohrazeném zdí. Tomu můžeme rozumět tak, že jsou chráněni před predátory a tedy neumírají, a také, že nemohou nikam utéci. Proto bude počet plodných v t-tém měsíci roven jejich počtu v předchozím měsíci zvětšenému o počet mladých, kteří se v předchozím měsíci narodili a během měsíce dospěli, y(t) = y(t − 1) + x(t − 1). (2) 1 Překlad E. Čecha. Citováno dle J. Bečvář a kol., Matematika ve středověké Evropě. Praha: Prometheus 2001, str. 277. 3 4 PROLOG měsíc t 0 1 2 3 4 5 6 7 8 9 10 11 12 počet juvenilních párů x(t) 0 1 1 2 3 5 8 13 21 34 55 89 144 počet plodných párů y(t) 1 1 2 3 5 8 13 21 34 55 89 144 233 celkový počet párů z(t) 1 2 3 5 8 13 21 34 55 89 144 233 377 Tabulka 1: Řešení Fibonacciovy úlohy o králících za předpokladu, že k rození dochází na začátku měsíce, počty zjišťujeme na konci měsíce, tj. používáme model (3). Rovnice (1) a (2) můžeme považovat za model růstu populace králíků; její aktuální velikost počítáme z velikosti v minulosti. Při matematickém modelování nějakých procesů je ovšem obvyklé usuzovat na budoucnost z přítomnosti. V rovnicích (1) a (2) budeme psát t+1 místo t, rovnice tedy přepíšeme do tvaru x(t + 1) = y(t), y(t + 1) = x(t) + y(t), nebo ve vektorovém zápisu x y (t + 1) = 0 1 1 1 x y (t). (3) Měsíc, ve kterém „kdosi umístil pár králíků na určitém místě , budeme považovat za nultý, onen „umístěný pár za dospělé. Máme tedy počáteční podmínku x(0) = 0, y(0) = 1. Odtud již můžeme postupně počítat počty x(t) a y(t) pro libovolné t = 1, 2, 3, . . . a z nich celkový počet párů z(t) = x(t) + y(t). Výpočet je shrnut v tabulce 1. Výsledek 377 párů odpovídá výsledku v Liber abaci.2 Jiná z možností, jak zadání porozumět, je mírně realističtější představa, že králíci se rodí kdykoliv, ale opět je sčítáme v určitý den měsíce. Při sčítání tedy mohou mít novorozenci, tj. králíci narození od předchozího sčítání, věk z intervalu [0, 1) a starší, ale dosud neplodní králíci věk z intervalu [1, 2). Při této interpretaci rozdělíme třídu juvenilních párů na dvě a označíme x0(t) počet novorozených párů a x1(t) počet neplodných párů věku alespoň jeden měsíc, ale méně než dva měsíce. Poněvadž novorozenci jsou bezprostředními potomky plodných párů, mladí jsou ti, kteří se v předchozím měsíci narodili, a počet plodných je počtem plodných z předchozího měsíce zvětšeným o počet mladých, kteří dosáhli věku aspoň dva měsíce, dostaneme model x0(t + 1) = y(t) x1(t + 1) = x0(t) y(t + 1) = x1(t)+y(t), který opět můžeme přepsat do vektorového tvaru   x0 x1 y   (t + 1) =   0 0 1 1 0 0 0 1 1     x0 x1 y   (t). (4) Při počátečních podmínkách x0(0) = 0, x1(0) = 0, y(0) = 1 a označení celkového počtu párů jako z(t) = x0(t) + x1(t) + y(t), dostaneme počty králíků, jak je uvedeno v tabulce 2. Výsledný počet párů králíků za rok je při této interpretaci téměř třikrát menší, než původní 2 To nemusí znamenat, že by si Fibonacci představoval rození na začátku měsíce a sčítání na jeho konci. Pravděpodobnější je, že si neuměl představit nulový věk a proto jeho novorozenci měli hned věk 1 a v následujícím měsíci tak byli dvouměsíční a tedy již plodní. FIBONACCIOVI KRÁLÍCI 5 měsíc t 0 1 2 3 4 5 6 7 8 9 10 11 12 počet novorozených párů x0(t) 0 1 1 1 2 3 4 6 9 13 19 28 41 počet neplodných párů x1(t) 0 0 1 1 1 2 3 4 6 9 13 19 28 počet plodných párů y(t) 1 1 1 2 3 4 6 9 13 19 28 41 60 celkový počet párů z(t) 1 2 3 4 6 9 13 19 28 41 60 88 129 Tabulka 2: Řešení Fibonacciovy úlohy o králících za předpokladu, že k rození dochází kdykoliv v průběhu měsíce a králíky sčítáme v pevně určený den měsíce, tj. používáme model (4). Fibonacciův výsledek. Prvním obecným poučením tedy může být to, že sestavení modelu růstu populace je potřebné věnovat pozornost, přesně formulovat a zdůvodnit předpoklady, za kterých je model sestaven. Různé modely téhož procesu mohou totiž dávat různé výsledky. Druhým závěrem je pozorování, že vývoj populace králíků můžeme popsat matematickým modelem (3) nebo (4), které jsou zvláštními případy obecného modelu ve tvaru n(t + 1) = An(t); (5) přitom složky vektoru n(t) vyjadřují časově závislé velikosti jednotlivých tříd, do kterých je populace strukturována (rozdělena) a A je nějaká matice. Fibonacciův model je krásný matematicky, není ovšem příliš realistický biologicky. Králíci neumírají, dospívají v přesně určených časech, plodí přesně určený počet potomků v pravidelných intervalech. Fibonacci samozřejmě nepředstíral, že popisuje vývoj populace králíků, vytvořil jakousi umělou skutečnost — jeho králíci žijí a množí se na „místě ohraženém zdí . Myšlenka modelovat pomocí rovnice (5) vývoj populace rozdělené na několik disjunktních tříd, přičemž čas plyne v diskrétních krocích, je však velmi plodná. Pokusíme se modelovat vývoj populace za realističtějších předpokladů. Ponecháme původní představu času plynoucího v diskrétních krocích (nejedná se tedy o čas fyzikální) a zvolíme nějakou časovou jednotku (ve Fibonacciově úloze jí byl jeden měsíc). Populaci si budeme představovat jako tvořenou velkým počtem jedinců (v případě organismů rozmnožujících se pohlavně budeme za „jedince považovat páry nebo samice). Každý z jedinců může být jednoho z typů — juvenilní (mladý, neplodný) nebo dospělý (plodný). Jinak jsou jedinci nerozlišitelní. V populaci probíhají tři procesy — rození (vznik nových jedinců), dospívání (maturace, přeměna juvenilního jedince na plodného) a umírání (nebo z jiného pohledu přežívání). Narození jedince, jeho přeměnu na plodného a jeho úmrtí považujeme za náhodné jevy. O umírání (přežívání) a dospívání budeme předpokládat, že se jedná o jevy stochasticky nezávislé. Označme σ1 . . . pravděpodobnost, že juvenilní jedinec přežije jedno období, σ2 . . . pravděpodobnost, že plodný jedinec přežije jedno období, γ . . . pravděpodobnost, že juvenilní jedinec během období dospěje, ϕ . . . střední počet potomků plodného jedince za jedno období. O pravděpodobnostech přežití σ1 a σ2, pravděpodobnosti maturace γ a fertilitě ϕ budeme předpokládat 0 < σ1 < 1, 0 ≤ σ2 < 1, 0 < γ ≤ 1, 0 < ϕ; (6) v reálně existující populaci totiž musí být možné, že se juvenilní jedinec dožije plodnosti (σ1 > 0, γ > 0) a že se nějací noví jedinci rodí (ϕ > 0), přežití nikdy není jisté (σ1 < 1, 6 PROLOG σ2 < 1). Nevylučujeme možnost σ2 = 0, tj. že jedinci po „produkci potomků (porodu, nakladení vajíček a podobně) hynou; taková populace se nazývá semelparní. Nevylučujeme však ani možnost σ2 > 0, tj. že dospělí jedinci plodí po delší úsek života; taková populace se nazývá iteroparní. Jedinci mohou dospívat bezprostředně po narození, tj. v čase kratším, než je zvolené období. V období po narození tedy takový jedinec, pokud nezemře, jistě dospěje, γ = 1. Jedinci z populace mohou dospívat i s jistým zpožděním, γ < 1. Zhruba řečeno, při délce časového kroku jeden rok jsou jednoleté organismy semelparní s bezprostředním dospíváním, drobní ptáci a savci jsou iteroparní s bezprostředním dospíváním, lososi nebo cikády jsou semelparní se zpožděným dospíváním, velcí ptáci a savci (včetně člověka) jsou iteroparní se zpožděným dospíváním. Snažíme se tedy modelovat dosti obecnou populaci. Označme dále n1(t), resp. n2(t), velikost (počet jedinců, populační hustotu, celkovou biomasu a podobně) části populace tvořené juvenilními, resp. plodnými, jedinci v t-tém časovém kroku. Juvenilní část populace je tvořena jedinci, kteří se za poslední období narodili, a jedinci, kteří již tuto třídu populace tvořili, přežili období a nedospěli v něm. Očekávaná velikost juvenilní části populace v následujícím období tedy bude n1(t + 1) = σ1(1 − γ)n1(t) + ϕn2(t). (7) Plodná část populace bude tvořena jedinci, kteří byli juvenilní, nezemřeli a dospěli, a jedinci, kteří již dospělí byli a přežili. Očekávaná velikost plodné části populace v následujícím období tedy bude n2(t + 1) = σ1γn1(t) + σ2n2(t). (8) Tyto rovnice opět můžeme přepsat ve vektorovém (maticovém) tvaru n(t + 1) = n1 n2 (t + 1) = σ1(1 − γ) ϕ σ1γ σ2 n1 n2 (t) = An(t), (9) kde jsme označili A = σ1(1 − γ) ϕ σ1γ σ2 . Poznamenejme ještě, že kdybychom připustili σ1 = σ2 = 1 a položili γ = ϕ = 1 (jedinci jistě přežívají, tj. neumírají, jistě během období dospějí a dospělí vždy vyprodukují právě jednoho potomka), dostaneme původní Fibonacciův model (3). Leslieho model růstu populace Patrick Holt Leslie (1900–1972) absolvoval bakalářské studium fyziologie na Christ Church College oxfordské university. Zdravotní problémy mu však neumožnily dokončit studium medicíny. Několik let pracoval v bakteriologické laboratoři na katedře patologie, pak ale svůj zájem obrátil k matematické statistice a přijal místo v nově založeném Odboru populací živočichů (Bureau of Animal Population). Za druhé světové války se toto centrum zabývalo metodami regulace hlodavců v zásobárnách obilí. Roku 1945 publikoval Leslie v časopisu Biometrika svůj nejslavnější článek On the use of matrices in certain population mathematics. V něm sestavil a analyzoval model růstu počtu samic v populaci potkanů (Rattus norvegicus); jeho model ovšem může být stejně dobře použit pro lidskou populaci. Za jedinou charakteristiku samice budeme považovat její věk. Ten je při narození nulový a postupně narůstá s plynutím času. Fyzikální čas považujeme za kontinuum, lze ho vyjádřit LESLIEHO MODEL RŮSTU POPULACE 7 reálným číslem. Věk ovšem obvykle vyjadřujeme číslem přirozeným; u lidí v novorozeneckém věku počtem dní, poté počtem týdnů, ještě později počtem měsíců, a nakonec počtem let od narození. U praxe vyjadřovat věk přirozeným číslem zůstaneme i při sestavování modelu růstu populace. Nejprve je tedy potřeba stanovit časovou jednotku. Může to být den, týden, měsíc, rok, pět let nebo desetiletí; záleží na tom, jak přesně (jemně) chceme věkovou strukturu populace popisovat. Spojitou veličinu věk rozdělíme na diskrétní věkové třídy — v první třídě budou samice ve věku od narození do jedné (časové jednotky), ve druhé samice od věku jedna do věku dvě, atd. Aby bylo přiřazení samic do věkové třídy jednoznačné, dohodneme se, že obecně v i-té věkové třídě jsou zařazeny samice, jejichž věk je aspoň i − 1 a méně než i, tj. samice, jejichž věk je z intervalu [i − 1, i).3 Označme nyní ni(t) množství samic v i-té věkové třídě v čase t (v časovém okamžiku t). Čas vyjadřujeme ve stejných jednotkách jako věk a uvažujeme ho pouze s celočíselnými hodnotami, t = 0, 1, 2, . . . . „Množství samic může být vyjádřeno počtem jedinců nebo nějakou jinou jednotkou, např. počtem desetitisíců jedinců, počtem jedinců vztaženým na jednotku plochy a podobně. Ze všech procesů, které probíhají v populaci se omezíme na dva nejdůležitější z hlediska jejího růstu — rození a umírání. Budeme předpokládat, že z hlediska těchto procesů se samice v jedné věkové třídě neliší jedna od druhé. Označme Pi pravděpodobnost jevu, že samice, která v čase t patřila do i-té věkové třídy, přežije interval jednotkové délky, tedy bude žít i v čase t + 1. V dostatečně velké populaci bude tato pravděpodobnost rovna relativní četnosti samic z i-té věkové třídy, které přežily jedno období, tedy Pi = ni+1(t + 1) ni(t) . Pokud bychom proces umírání považovali za deterministický, lze parametr Pi interpretovat jako podíl samic, které během časového intervaly (t, t + 1] nezemřely, mezi všemi samicemi, které v čase t patřily do i-té věkové třídy. Hodnoty Pi se nazývají projekční koeficienty — promítají totiž aktuální množství samic na jejich množství v následujícím časovém období. Budeme předpokládat, že parametr Pi nezávisí na čase, že přežití jednoho období závisí pouze na věku, nikoliv na nějakých vnějších vlivech. Dále budeme předpokládat, že existuje nějaký maximální možný věk, kterého se lze dožít, ale který není možné překročit. Tedy že existuje přirozené číslo k takové, že Pi > 0 pro každé i = 1, 2, 3, . . . , k − 1 a Pk = 0. To znamená, že uvažujeme právě k věkových tříd. Z předchozí rovnosti nyní vyjádříme ni+1(t + 1) = Pini(t), i = 1, 2, 3, . . . , k − 1. (10) Poněvadž samice z i-té věkové třídy může (ale nemusí) uhynout během časového intervalu jednotkové délky, platí ni+1(t + 1) ≤ ni(t). Projekční koeficienty Pi tedy splňují podmínky 0 < Pi ≤ 1, i = 1, 2, 3, . . . , k − 1. (11) 3 Uvedené rozdělení do věkových tříd samozřejmě není jediné možné. Stejně dobře by v první věkové třídě mohly být samice věku z intervalu [0, 1], ve druhé samice věku z intervalu (1, 2]; obecně v i-té věkové třídě by byly samice věku z intervalu (i − 1, i], i = 2, 3, . . . . Jiná možnost třídění podle věku by mohla být [0, 1 2 ], (1 2 , 3 2 ], (3 2 , 5 2 ], atd. Věkové třídy uvažované v t textu mají jistou „estetickou přednost — intervaly věku v každé třídě jsou stejného typu, jsou polouzavřené zleva. 8 PROLOG O procesu rození budeme předpokládat, že počet dcer, které se narodí samicím z i-té věkové třídy během jednotkového časového intervalu je úměrný počtu těchto samic. Formálněji: počet novorozených samiček, tj. samic z první věkové třídy, které se narodí v časovém intervalu (t.t + 1] a jejichž matky patřily v čase t do i-té věkové třídy je roven Fini(t), kde Fi je nějaká nezáporná konstanta. Hodnotu Fi lze považovat za očekávaný počet živých dcer jedné typické samice z i-té věkové třídy za jednotkový časový interval; někdy se nazývá míra reprodukce ve věku i. Celkový počet novorozených samic v čase t + 1 tedy bude n1(t + 1) = F1n1(t) + F2n2(t) + · · · + Fknk(t) = k i=1 Fini(t). (12) V přirozené populaci nebudou mít potomky samice velmi mladé (nedospělé) a samice po menopauze. Existují tedy jisté hodnoty m, M ∈ {1, 2, . . . , k} takové, že m ≤ M a 0 = F1 = F2 = · · · = Fm−1 = FM+1 = FM+2 = · · · = Fk, Fm > 0, Fm+1 > 0, Fm+2 > 0, . . . , FM−1 > 0, FM > 0; (13) m označuje věk dosažení pohlavní dospělosti, M věk menopauzy; jsou-li samice plodné celý život, je M = k. Model vývoje populace je tedy dán rovnostmi (10), (12), tj. n1(t + 1) = F1n1(t) + F2n2(t) + F3n3(t) + · · · + Fk−1nk−1(t) + Fknk(t) n2(t + 1) = P1n1(t) n3(t + 1) = P2n2(t) ... nk(t + 1) = Pk−1nk−1(t). (14) Při označení A =          F1 F2 F3 . . . Fk−1 Fk P1 0 0 . . . 0 0 0 P2 0 . . . 0 0 ... ... ... ... ... ... 0 0 0 . . . 0 0 0 0 0 . . . Pk−1 0          , n(t) =          n1(t) n2(t) n3(t) ... nk−1(t) nk(t)          lze rovnosti (14) zapsat v maticovém tvaru n(t + 1) = An(t), tedy ve tvaru stejném, jako model Fibonacciových králíků (5). Matice A uvedeného tvaru se nazývá Leslieho matice. Na závěr ještě shrneme předpoklady, za kterých byl model vývoje populace (14) sestaven: • je adekvátní klasifikovat jedince podle věku, • informace o věku jedinců uvnitř věkové třídy je irelevantní, • plodnosti a úmrtnosti závisí pouze na věku. LESLIEHO MODEL RŮSTU POPULACE 9 t 0 1 2 3 4 5 1 2 3 4 5 n n1 n2 n3 Obrázek 1: Vývoj populace podle modelu (15). Ještě poznamenejme, že stejným způsobem lze modelovat nejen populaci savců, ale i populace jiných obratlovců, u nichž plodnost a přežívání závisí především na věku. Vždy je důležité zvolit časové měřítko (u velkých savců roky až desetiletí, u drobných savců měsíce až roky) a vhodnou jednotku velikosti populace (může to být počet jedinců, populační hustota, celková biomasa a podobně). Příklad Uvažujme jednoduchý Leslieho model populace strukturované do tří věkových tříd, kde jsou plodní jedinci druhé a třetí věkové třídy. Jedná se tedy o model (14) s k = M = 3, m = 2. Položme konkrétně F2 = 1, F3 = 5, P1 = 0,3, P2 = 0,5. Budeme předpokládat, že na počátku se populace skládala pouze z jednotkového množství jedinců třetí věkové třídy. Jedná se tedy o model   n1(t + 1) n2(t + 1) n3(t + 1)   =   0 1 5 0,3 0 0 0 0,5 0     n1(t) n2(t) n3(t)   ,   n1(0) n2(0) n3(0)   =   0 0 1   . (15) Snadno spočítáme, že   n1(1) n2(1) n3(1)   =   5 0 0   ,   n1(2) n2(2) n3(2)   =   0,0 1,5 0,0   ,   n1(3) n2(3) n3(3)   =   1,50 0,00 0,75   ,   n1(4) n2(4) n3(4)   =   3,75 0,40 0,00   ,   n1(5) n2(5) n3(5)   =   0,450 1,125 0,225   a tak můžeme pokračovat dále. Tento výsledek nám ovšem neposkytuje žádnou souhrnnou informaci, která by umožnila nějaký vhled do vývoje populace. A to ani v případě, že ho zobrazíme graficky, viz obr. 1. Vypočítáme proto vývoj populace podle modelu (15) na delším časovém úseku. Výsledek je prezentován graficky na obr. 2, na svislé ose (velikosti složek populace) je logaritmická 10 PROLOG t10 20 30 40 50 n1 n2 n3 n 0,1 1 10 Obrázek 2: Vývoj populace podle modelu (15). Na svislé ose je logaritmická stupnice. t 10 20 30 40 500 n1/N n2/N n3/N relativní četnost 0,2 0,4 0,6 0,8 1,0 Obrázek 3: Vývoj věkové struktury populace podle modelu (15), N = n1 + n2 + n3 LESLIEHO MODEL RŮSTU POPULACE 11 E t10 20 300 T N 1 2 3 4 T n1/N 1,0 0,5 0 Et 10 20 30 T n2/N 1,0 0,5 0 Et 10 20 30 T n3/N 1,0 0,5 0 Et 10 20 30 Obrázek 4: Vliv počátečních podmínek na vývoj populace. stupnice. Vidíme, že kolísání velikostí jednotlivých věkových tříd se po jisté době (asi po 30 časových jednotkách) ustálí a velikost populace roste exponenciálně, velikosti jednotlivých věkových tříd rostou stejně rychle. To znamená, že po jisté době vývoje se věková struktura populace stabilizuje, relativní zastoupení jednotlivých věkových tříd zůstává konstantní. Tuto úvahu potvrzuje i grafické znázornění vývoje relativních velikostí jednotlivých věkových tříd na obr. 3. Ještě poznamenejme, že velikost populace modelovaná rovnostmi (15) s časem roste, relativní zastoupení jednotlivých věkových tříd ve věkově stabilizované populaci (tj. v populaci po dostatečně dlouhém vývoji) se s rostoucím věkem snižuje. Vliv počátečních podmínek V modelu (15) budeme měnit počáteční podmínky. Postupně budeme předpokládat, že na počátku se populace skládala z jednotkového množství jedinců právě jedné ze tří věkových tříd, nebo že v populaci jednotkové velikosti mohly být zastoupeny všechna věkové třídy. Uvažujme tedy první rovnost z modelu (15) spolu s některou z počátečních podmínek   n1(0) n2(0) n3(0)   =   1 0 0   ,   n1(0) n2(0) n3(0)   =   0 1 0   ,   n1(0) n2(0) n3(0)   =   0 0 1   ,   n1(0) n2(0) n3(0)   =   ξ1 ξ2 ξ3   , (16) kde ξ1 ≥ 0, ξ2 ≥ 0, ξ3 ≥ 0, ξ1 + ξ2 + ξ3 = 0. Výsledky jsou graficky znázorněny ma obr. 4. Vidíme, že ve všech případech po asi 25 časových jednotkách velikost populace roste exponenciálně a že relativní velikosti jednotlivých věkových tříd se ustálí na stejných hodnotách. Počáteční struktura populace tedy neovlivňuje ani rychlost růstu populace, ani její stabilizovanou strukturu. Ovlivňuje pouze celkovou velikost populace. 12 PROLOG původní F2 F3 P2 P1 0 E t10 20 30 40 50 TN 2 4 6 Obrázek 5: Vývoj celkové velikosti populace (N = n1 +n2 +n3) podle modelu (15), ve kterém je právě jeden z parametrů zmenšen o 10%. Vliv změny parametrů V modelu (15) postupně zmenšíme každý z parametrů o 10%. Vývoj celkové velikosti populace v takových případech vidíme na obr. 5. Zmenšení plodnosti ve druhé věkové třídě zpomalí růst populace, ale celková velikost populace stále roste exponenciálně. Zmenšení plodnosti nejstarších jedinců nebo zmenšení některé z pravděpodobností přežití však způsobí, že populace bude vymírat; vymírání je nejrychlejší v případě zvětšení úmrtnosti (omezení přežití) nejmladší věkové třídy. Z těchto výsledků vidíme, že rychlost růstu populace je nejcitlivější na změny novorozenecké úmrtnosti, nejméně je citlivá na změny plodnosti střední věkové třídy. Vliv kolísání velikosti parametrů Z předchozího výpočtu víme, že zmenšení plodnosti třetí věkové třídy vede k vymření populace, zmenšení plodnosti druhé věkové třídy vede ke zpomalení růstu populace. Podívejme se nyní, jak se bude vyvíjet velikost populace, pokud se plodnosti mění v průběhu času. Uvažujme tedy model   n1(t + 1) n2(t + 1) n3(t + 1)   =   0 h(t) 5h(t) 0,3 0 0 0 0,5 0     n1(t) n2(t) n3(t)   , (17) kde h(t) = 1 + 1 10 cos t. Plodnosti tedy kolísají kolem původních hodnot v rozmezí ±10%. Vývoj populace budeme opět vyšetřovat pro různé počáteční hodnoty (16). Výsledný růst populace je znázorněn na obr. 6. Nyní vidíme, že struktura populace se neustálí na nějakých hodnotách relativního zastoupení jednotlivých věkových tříd, ale kolem těchto hodnot kolísají. Populace s parametry závislými na její velikosti Přirozené populace žijí v prostředí s omezenými zdroji. Budeme pro jednoduchost předpokládat, že velká populace spotřebovává více zdrojů a že jejich následný nedostatek způsobí LESLIEHO MODEL RŮSTU POPULACE 13 E t10 20 30 40 50 T1,0 0,5 0 n1/N n2/N n3/N E t10 20 30 40 50 T 5 4 3 2 1 n1 n2 n3 Obrázek 6: Vývoj populace podle modelu (17). Nahoře velikosti jednotlivých věkových tříd, dole relativní zastoupení jednotlivých věkových tříd v populaci. zmenšení plodností jednotlivých věkových tříd. Označme tedy N = n1 + n2 + n3 a uvažujme model   n1(t + 1) n2(t + 1) n3(t + 1)   =   0 g(N) 5g(N) 0,3 0 0 0 0,5 0     n1(t) n2(t) n3(t)   ,   n1(0) n2(0) n3(0)   =   0 0 1   . (18) Přitom funkce g je nezáporná, klesající a taková, že g(0) = 1; pokud by tedy nedocházelo k uvažované vnitrodruhové konkurenci, populace by se vyvíjela podle modelu (15). Zvolme pro určitost g(N) = e−bN , kde b > 0. Vývoj celkové velikosti populace i jejích jednotlivých složek při použití modelu (18) s touto funkcí g a s parametrem b = 0,005 je znázorněn na obr. 7 nahoře. Vidíme, že (absolutní) velikosti jednotlivých věkových tříd se ustálí na jistých hodnotách a také celková velikost populace naroste do jisté limitní hodnoty. Tato mezní velikost populace představuje kapacitu (úživnost) prostředí pro modelovanou populaci. Na obr. 7 dole jsou zobrazeny průběhy velikostí jednotlivých složek takto modelované populace pro různé počáteční podmínky (16). Tento výsledek naznačuje, že limitní velikosti jednotlivých věkových tříd populace nezávisí na počátečních podmínkách. Ty ovlivňují pouze rychlost konvergence. Uvažujme ještě populaci, v níž plodnost jedinců druhé věkové třídy může být větší než jednotková. Vývoj populace budeme tedy modelovat rovnostmi (18), funkce g však bude dána rovností g(N) = Re−bN , kde R ≥ 1. vývoj celkové velikosti populace pro různé hodnoty parametru R vidíme na obr. 8. Pro R = 3 se velikost populace ustálí na hodnotě kapacity prostředí a velikost populace konverguje k této limitní hodnotě s tlumenými oscilacemi. Pro R = 15 a pro R = 45 populace po jisté době vývoje začne kolem mezní hodnoty pravidelně oscilovat (velikost populace je od 14 PROLOG Et 0 100 200 300 400 500 T 0 2 4 6 8 N n1/N n2/N n3/N E 0 100 200 t T n3/N E 0 100 200 t T n2/N E 0 100 200 t T n1/N 5 3 1 Obrázek 7: Nahoře: Vývoj populace s interní variabilitou podle modelu (18) s g(N) = e−0,005N . Dole: Vývoj velikosti jednotlivých věkových tříd populace podle první rovnosti modelu (18) s různými počátečními podmínkami (16). — první, — druhá, — třetí, — čtvrtá podmínka; se čtvrtou z podmínek bylo provedeno 20 simulací. jistého času periodická, v případě R = 15 je perioda rovna 2, v případě R = 45 je rovna 4), pro R = 250 začne kolísat a v tomto kolísaní již není žádná pravidelnost vidět. Vektory, matice a operace s nimi Vektory budeme chápat jako sloupcové. Budeme je označovat tučnými malými písmeny, jejich složky budeme převážně značit stejným písmenem jako je označen vektor, doplněným dolním indexem, nebo znakem vektoru v kulaté závorce s dolním indexem. Tedy v =      v1 v2 ... vk      , (v)i = vi. Symbolem 1 označíme vektor, jehož všechny složky jsou rovny jedné, symbolem o vektor, jehož všechny složky jsou rovny 0, 1 =      1 1 ... 1      , o =      0 0 ... 0      . VEKTORY, MATICE A OPERACE S NIMI 15 R = 250 R = 45 R = 15 R = 3 0 10 20 30 40 50 60 0 10 3020 40 50 60 0 10 20 30 40 50 60 0 10 20 30 40 50 60 0 3 600 7 200 10 800 14 400 0 50 100 150 200 0 500 1 000 1 500 2 000 0 280 560 840 1 120 E t T N E t TN E t TN E t TN Obrázek 8: Vývoj populace s interní variabilitou 16 PROLOG Matice budeme převážně označovat velkými bezpatkovými písmeny, jejich složky stejným písmenem malým s dvojicí indexů, případně znakem matice v kulatých závorkách s dvojicí indexů; první index je řádkový, druhý sloupcový. Matice typu m × n je A =      a11 a12 . . . a1n a21 a22 . . . a2n ... ... ... ... am1 am2 . . . amn      , (A)ij = aij. Vektor s k složkami (k-rozměrný vektor) je tedy matice typu k × 1. Symboly O a I označíme nulovou a jednotkovou matici, O =      0 0 . . . 0 0 0 . . . 0 ... ... ... ... 0 0 . . . 0      , I =      1 0 . . . 0 0 1 . . . 0 ... ... ... ... 0 0 . . . 1      , (O)ij = 0, (I)ij = δij = 1, i = j, 0, i = j; přitom δij je Kroneckerův symbol. Bude-li potřebné zdůraznit, že nulová nebo jednotková matice je současně čtvercovou maticí dimense n, budeme psát On nebo In. Součet matic A, B stejného typu je definován po složkách, tj.      a11 a12 . . . a1n a21 a22 . . . a2n ... ... ... ... am1 am2 . . . amn      +      b11 b12 . . . b1n b21 b22 . . . b2n ... ... ... ... bm1 bm2 . . . bmn      = =      a11 + b11 a12 + b12 . . . a1n + b1n a21 + b21 a22 + b22 . . . a2n + b2n ... ... ... ... am1 + bm1 am2 + bm2 . . . amn + bmn      , (A + B)ij = aij + bij. Násobek matice skalárem γ je matice γA stejného typu, γ      a11 a12 . . . a1n a21 a22 . . . a2n ... ... ... ... am1 am2 . . . amn      =      γa11 γa12 . . . γa1n γa21 γa22 . . . γa2n ... ... ... ... γam1 γam2 . . . γamn      , (γA)ij = γaij. Pro lineární kombinaci matic A, B stejného typu platí (αA + βB)ij = αaij + βbij. Také platí (α + β)A ij = (α + β)aij. Násobení matic: Násobením matice A typu m×n maticí B typu n×p (zprava) dostaneme matici C = AB typu m × p, pro jejíž složky platí cij = n ι=1 aiιbιj. VEKTORY, MATICE A OPERACE S NIMI 17 Pro násobení matice A typu m × n vektorem v o n složkách platí (Av)i = n j=1 aijvj. Transpozice matice: Matice transponovaná k matici A typu m×n je matice AT typu n×m, pro niž platí AT ij = aji. Pro násobení matic platí (AB)T = BTAT. Skalární součin vektorů v, w o k složkách je definován jako maticový součin vT w = k i=1 viwi. Hadamardův součin matic A, B stejného typu je „součin po složkách , tj.      a11 a12 . . . a1n a21 a22 . . . a2n ... ... ... ... am1 am2 . . . amn      ◦      b11 b12 . . . b1n b21 b22 . . . b2n ... ... ... ... bm1 bm2 . . . bmn      =      a11b11 a12b12 . . . a1nb1n a21b21 a22b22 . . . a2nb2n ... ... ... ... am1bm1 am2bm2 . . . amnbmn      , (A ◦ B)ij = aijbij. Čtvercovou diagonální matici, která má v diagonále složky vektoru v značíme diag v. Je tedy diag v =      v1 0 . . . 0 0 v2 . . . 0 ... ... ... ... 0 0 . . . vk      . Pro vektory v, w stejné dimenze platí v ◦ w = v diag w = (diag v)w. Kroneckerův součin matic: Nechť matice X = (xij) je typu µ × ν a matice Y = (yij) je typu κ×λ. Jejich Kroneckerův součin X⊗Y je matice typu µκ×νλ, kterou lze blokově zapsat ve tvaru X ⊗ Y =      x11Y x12Y . . . x1νY x21Y x22Y . . . x2νY ... ... ... ... xµ1Y xµ2Y . . . xµνY      . Například tedy   2 0 0 1 0 −1 0 3 0   ⊗ 3 1 −2 2 =         6 2 0 0 0 0 −4 4 0 0 0 0 3 1 0 0 −3 −1 −2 2 0 0 2 −2 0 0 9 3 0 0 0 0 −6 6 0 0         . 18 PROLOG Operace vec „poskládá sloupce matice nad sebe , přesněji: z matice A typu m × n vytvoří mn-rozměrný vektor vec A, vec A = vec      a11 a12 . . . a1n a21 a22 . . . a2n ... ... ... ... am1 am2 . . . amn      =                           a11 a21 ... am1 a12 a22 ... am2 ... a1n a2n ... amn                           , (vec A)k = aij kde i = 1 − k m , j = k m ; přitom [ξ] označuje celou část z reálného čísla ξ. S využitím operací vec a ⊗ můžeme přepsat součin matice A typu m × n a n-rozměrného vektoru: Av =      a11v1 + a12v2 + · · · + a1nvn a21v1 + a22v2 + · · · + a2nvn ... am1v1 + am2v2 + · · · + amnvn      =      v1a11 + v2a12 + · · · + vna1n v1a21 + v2a22 + · · · + vna2n ... v1am1 + v2am2 + · · · + vnamn      = =      v1 0 . . . 0 0 v1 . . . 0 ... ... ... ... 0 0 . . . v1           a11 a21 ... am1      +      v2 0 . . . 0 0 v2 . . . 0 ... ... ... ... 0 0 . . . v2           a12 a22 ... am2      + · · · · · · +      vn 0 . . . 0 0 vn . . . 0 ... ... ... ... 0 0 . . . vn           a1n a2n ... amn      = = v1Im      a11 a21 ... am1      + v2Im      a12 a22 ... am2      + · · · + vnIm      a1n a2n ... amn      = (vT ⊗ Im) vec A. Symbolem |A| označíme matici, jejíž složky jsou absolutními hodnotami složek matice A, |A| ij = |aij|. Podobně |v| označuje vektor se složkami |v| i = |vi|. Euklidovská norma k-rozměrného vektoru v je definována vztahem v 2 = k i=1 v2 i = √ vTv , VEKTORY, MATICE A OPERACE S NIMI 19 součtová (taxíkařská) norma k-rozměrného vektoru v je definována vztahem v 1 = k i=1 |vi| = 1T |v|. 20 PROLOG Kapitola 1 Konstrukce modelů 1.1 Stavové proměnné „Stav nějakého systému určuje jeho „chování . Například v mechanice je stav systému definován pomocí poloh a hybností všech částic, které ho tvoří; v etologii je stav jedince vyjádřen jeho „motivací ; stav ekosystému je popsán množstvím hmoty a energie, kterou si jeho jednotlivé složky vyměňují; v demografii je stav populace dán velikostí jednotlivých tříd (např. věkových), do nichž můžeme jedince rozdělit a podobně. 1.1.1 Zadehova teorie stavové proměnné Výchozím bodem teorie je pojem abstraktního objektu O, který interaguje s okolím pomocí stimulů (podnětů, buzení; excitation), které na něho působí, a odezev (response), kterými se projevuje navenek. Předpokládejme, že stimuly i odezvy lze nějak kvantifikovat. Přesněji: nechť objekt O pozorujeme v časovém intervalu [t0, t1), kde t0 < t1 ≤ ∞, a stimuly a odezvy objektu v tomto časovém intervalu lze popsat funkcemi e : [t0, t1) → E a r : [t0, t1) → R, kde E a R jsou nějaké podmnožiny Banachova prostoru. Pak lze objekt O ztotožnit s pozorovanými stimuly a odezvami, tj. O = t, e(t) : t0 ≤ t < t1 , t, r(t) : t0 ≤ t < t1 . Pro zjednodušení zápisu zavedeme pro podmnožinu Y Banachova prostoru, pro interval I reálných čísel a pro funkci y : R → Y označení yI = t, y(t) : t ∈ I . Pak můžeme psát O = e[t0,t1), r[t0,t1) . Při tomto pojetí představuje experiment vyvolání určitých stimulů e[t0,t1) a pozorování odezev r[t0,t1) objektu O. Základním předpokladem je, aby objekt O byl determinovaný1, tj. aby odezva byla stimulem jednoznačně určena. Požadujeme tedy, aby existovalo zobrazení ϕ z množiny 2R×E do množiny 2R×R takové, že r[t0,t1) = ϕ e[t0,t1) ; 1 Determinovaný objekt není totéž, co deterministický. Pozorované funkce e a r mohou být realizací nějaké náhodné funkce. 21 22 KAPITOLA 1. KONSTRUKCE MODELŮ symbol 2Y označuje potenční množinu (množinu podmnožin) množiny Y . Funkce e a r je obtížné získat (pozorovat, měřit), a pokud se to podaří, obtížně se s nimi pracuje. Jedno z nabízejících se zjednodušení spočívá v uvažování okamžitých stimulů e(t) a odezev r(t) pro t ∈ [t0, t1). Determinovanost objektu by pak znamenala, že existuje zobrazení ψ : E → R převádějící stimulus v okamžiku t na okamžitou odezvu r, r(t) = ψ e(t) . Stejný stimulus v různých časových okamžicích však často vyvolá různou odezvu, což znamená, že mezi stimulem a odezvou je nějaká zprostředkující proměnná x, která se v průběhu času mění. Tato proměnná nemusí být pozorovatelná, může, ale nemusí nějak odpovídat struktuře objektu; představuje jakousi hypotézu o uvažovaném objektu O, nějak vyjadřuje jeho stav. Nazývá se stavová proměnná. Stavovou proměnnou chápeme jako funkci času, x : R → X, kde X je opět nějaká podmnožina Banachova prostoru. Tato funkce je obecně náhodná, její hodnoty jsou dány rozložením pravděpodobnosti. V tomto textu však budeme uvažovat pouze deterministické objekty, tj. takové, že stavová proměnná má v každém čase t ∈ [t0, t1) nulový rozptyl a proto ji lze považovat za funkci klasickou (nenáhodnou). Na stavovou proměnnou x klademe dva požadavky: 1. Odezva v časovém okamžiku t ∈ [t0, t1) je jednoznačně určena stavem a stimulem v tomto čase t, tj. existuje zobrazení G : X × E → R (stimulus-state-response function) takové, že r(t) = G x(t), e(t) . (1.1) 2. Stav v nějakém časovém okamžiku je jednoznačně určen stavem v nějakém předchozím čase a stimuly, které objekt od té doby dostal, tj. existuje takové zobrazení F z množiny X × 2R×E do množiny X, že x(t + ∆t) = F x(t), e[t,t+∆t) . (1.2) Zobrazení F se nazývá přechodová funkce (state-transition function). Požadavek 1. říká, že ke znalosti objektu stačí znát jeho stav a stimuly, které na něho působí. Je splněn zejména tehdy, když jsou stavové proměnné přímo pozorovatelné. V takovém případě lze okamžitou odezvu přímo ztotožnit se stavem a rovnost (1.1) má tvar r(t) = x(t). Požadavek 2. je omezující; u skutečných objektů může stav x(t + ∆t) záviset také na historii, tj. hodnotách x(τ) pro τ < t0, nebo na budoucnosti2, tj na hodnotách x(τ) pro τ > t0. Budeme se tedy zabývat pouze neanticipativními systémy bez paměti. Rovnice (1.2) pro neznámou funkci x spolu s počáteční podmínkou x(t0) = x0 představuje model časového vývoje objektu O. Základním problémem matematického modelování je tedy nalezení (nebo konstrukce) přechodové funkce F. Speciální třídu modelů tvoří maticové modely. Jsou to modely, pro něž X = Rk a existuje matice A typu k × k taková, že přechodová funkce má tvar F x(t), e[t,t+∆t) = Ax(t). To neznamená, že by funkce F byla lineární v první proměnné. Matice A může záviset na stavu x(t). Je-li navíc ∆t > 0 pevně zvoleno, dostaneme diskrétní maticový model. Prvky 2 Taková závislost nemusí znamenat porušení kauzality, neboť přechodová funkce nevyjadřuje příčinnost ale pouze funkční závislost 1.1. STAVOVÉ PROMĚNNÉ 23 matice A pak vyjadřují stimuly, které působí v časovém intervalu [t, t + ∆t) na objekt O, který byl v čase t ve stavu x = x(t). Prvky matice A tedy závisí na stavu x a čase t, tj. A = A(x, t). Při vhodné volbě časové jednotky lze dosáhnout toho, že ∆t = 1 a rovnici (1.2) lze zapsat ve tvaru x(t + 1) = A x(t), t x(t). (1.3) Časová jednotka ∆t se nazývá projekční interval, rovnice (1.3) se nazývá projekční rovnice. Pokud závislost matice A na čase t je nekonstantní, tj. pro nějaký stav x0 ∈ Rk existují časy τ1, τ2 ∈ [t0, t1 − 1] takové, že τ1 = τ2 a A(x0, τ1) = A(x0, τ2), mluvíme o maticových modelech s externí variabilitou. Pokud matice A skutečně závisí na stavu x, tj. pro nějaký čas t ∈ [t0, t1 − 1] existují stavy x1, x2 ∈ Rk takové, že x1 = x2 a A(x1, t) = A(x2, t), mluvíme o maticových modelech s interní variabilitou. 1.1.2 Stavové proměnné v populačních modelech Populaci můžeme chápat jako objekt složený z jiných objektů. Jinak řečeno, každá populace je tvořena jedinci a každý jedinec je v nějakém stavu. Stav jedince (individua) budeme nazývat istav. Může jím být např. věk, velikost, vývojové stadium, obývaná lokalita, využitelná energie (tukové zásoby) a podobně. I-stav určuje odezvu jedince na stimulus. Ovšem v populačních modelech není zkoumaným objektem jedinec, ale populace. Stav populace nazveme p-stav. Pokud jsou splněny následující podmínky, 1. všichni jedinci jsou ovlivňováni týmž prostředím, 2. vliv populace na prostředí je součtem vlivů jedinců, pak lze p-stav vyjádřit jako rozložení i-stavů. Například, je-li jediným uvažovaným i-stavem vývojové stadium hmyzu, p-stavem v čase t může být čtyřrozměrný vektor, jehož složky jsou počty vajíček, larev, kukel a dospělců; je-li i-stavem věk, může být p-stavem v čase t integrovatelná funkce u : R+ → R+, přičemž u(a) vyjadřuje hustotu jedinců věku a, tj. takovou funkci, aby množství jedinců ve věku od a1 do a2 bylo rovno integrálu a2 a1 u(a)da (v tomto případě by však nemohlo jít o maticový model). Typickými stimuly působícími na populaci jsou rození, umírání, dospívání, migrace a podobně. Tyto stimuly jsou také projevy i-stavů. Je-li například i-stavem věk jedince, pak staří jedinci umírají s jinou pravděpodobností než mladí, v jistém věku jsou jedinci plodnější než ve stáří nebo bezprostředně po narození atd. Avšak u některých organismů plodnost nezávisí na věku, ale na velikosti (u rostlin, korýšů) nebo na věku i velikosti (u ryb); někdy přežití stejně starých jedinců závisí na jejich původu (např. zda rostlina vyrostla ze semene nebo je klonem mateřské rostliny); přechod z jednoho stadia do následujícího u některých hmyzů nezávisí na věku, ale na teplotě okolního prostředí (v tomto případě by p-stav „teplota okolního prostředí nebyl součtem nebo rozložením nějakých i-stavů) a podobně. Tyto skutečnosti ukazují, že výběr i-stavových proměnných je kritickým krokem při tvorbě modelu. Uvažujme nyní pro určitost jako stimulus proces rození. Ten lze kvantifikovat jako počet potomků za jednotku času (tedy veličinu nabývající nezáporných celočíselných hodnot), nebo 24 KAPITOLA 1. KONSTRUKCE MODELŮ i-stavová proměnná spojitá diskrétní spojitý regresní analýza analýza variance stimulus diskrétní logistická regrese kontingenční tabulky Tabulka 1.1: Statistické metody pro vyhodnocení závislosti stimulu na i-stavu celkovou biomasu potomků za jednotku času (nezáporné reálné číslo). Plodnost může být projevem i-stavu věk nebo velikost (kladné reálné číslo), ale také např. postavením v hierarchii skupiny (přirozené číslo) atd. Stimulus i i-stav tedy mohou být spojité i diskrétní veličiny. Jako vhodná i-stavová veličina určující stimulus by měla být vybrána ta, která na základě experimentů nebo pozorování vykazuje statisticky průkazný vliv na uvažovaný stimulus. Možnost volby statistické metody k vyhodnocení vlivu i-stavu na stimulus podle charakteru i-stavové veličiny a příslušného stimulu je shrnuta v tabulce 1.1. 1.2 Modely disperse Budeme se zabývat modely populace, u níž rozlišujeme dva i-stavy, z nichž jeden vyjadřuje místo výskytu jedince. Předpokládejme tedy, že populace strukturovaná podle nějakého kritéria (věku, velikosti, hmotnosti, plodnosti, stadia a podobně) je rozptýlena mezi několik lokalit (regionů, stanovišť, potravních ostrůvků a podobně). Mezi těmito lokalitami se jedinci z populace mohou přemisťovat. V takovém případě se mluví o metapopulacích (v ekologii) nebo o multiregionálních modelech (v demografii). Nejjednodušší případ disperse (migrace, šíření, rozptylu) mezi lokalitami je difúze. Při ní je množství jedinců přecházejících z jedné lokality na jinou úměrné velikosti populace na výchozí lokalitě, tj. existuje pravděpodobnost, že jedinec svou lokalitu opustí. V tomto oddílu budeme předpokládat, že tato pravděpodobnost nezávisí ani na velikosti populace, ani na nějakých vnějších vlivech. Nechť je tedy populace strukturována do s stadií a rozptýlena na p lokalit. Označme n ( ) i velikost části populace tvořené jedinci i-tého stadia na -té lokalitě. Zavedeme vektory n( ) =       n ( ) 1 n ( ) 2 ... n ( ) s       , ni =       n (1) i n (2) i ... n (p) i .       Vektor n( ) vyjadřuje strukturu populace na -té lokalitě, vektor ni vyjadřuje rozdělení jedinců i-tého stadia mezi lokalitami. Vektor popisující velikost celé populace můžeme vyjádřit dvěma způsoby — buď nejprve všechna stadia na jedné lokalitě, pak na druhé atd., nebo nejprve jedinci prvního stadia na jednotlivých lokalitách, pak jedinci druhého stadia atd. Struktura 1.2. MODELY DISPERSE 25 populace tedy může být vyjádřena buď vektorem n =      n(1) n(2) ... n(p)      =                              n (1) 1 n (1) 2 ... n (1) s n (2) 1 n (2) 2 ... n (2) s ... n (p) 1 n (p) 2 ... n (p) s                              nebo vektorem n =      n1 n2 ... ns      =                              n (1) 1 n (2) 1 ... n (p) 1 n (1) 2 n (2) 2 ... n (p) 2 ... n (1) s n (2) s ... n (p) s                              . (1.4) 1.2.1 Jednoduchý model difúze Vyjdeme ze zjednodušujících předpokladů: 1. Jednotlivé lokality jsou stejně kvalitní, tj. „demografické parametry jsou na všech stejné. 2. Pravděpodobnost emigrace z lokality závisí pouze na stadiu emigrujícího jedince, nikoliv na lokalitě. 3. Pravděpodobnost, že migrující jedinec přežije cestu závisí pouze na výchozí a cílové lokalitě, nikoliv na stadiu jedince. 4. Emigrace a přežití při migraci jsou jevy stochasticky nezávislé. Velikosti částí populace n ( ) i budeme vyjadřovat v časech t = 0, 1, 2, . . . . Pro popis struktury populace zvolíme první z možností (1.4). Dále budeme pro jednoduchost předpokládat, že k „demografické události , tj. k rození, k přechodu mezi stadii nebo k úmrtí, dochází bezprostředně po uvedených časech, nebo realističtěji řečeno, v časových intervalech (t, t + ε), kde ε je kladné malé číslo. K difúzi bude poté docházet v průběhu časových intervalů (t + ε, t + 1). Označme m ( ) i (t) velikost části populace tvořené jedinci i-tého stadia na -té lokalitě v čase t + ε, tedy bezprostředně po „demografické události . Vývoj populace tedy můžeme schematicky vyjádřit následujícím obrázkem. t0 1 2 3 n ( ) i (0) m ( ) i (0) n ( ) i (1) m ( ) i (1) n ( ) i (2) m ( ) i (2) n ( ) i (3) m ( ) i (3) Ec c c cc c c c · · · 26 KAPITOLA 1. KONSTRUKCE MODELŮ Přitom červené úsečky vyjadřují „demografické události , zelené difúzi. Nechť „demografické události na každé z lokalit popisuje matice A = (aij). To znamená, že m( ) (t) = An( ) (t), tj. m ( ) i (t) = s j=1 aijn ( ) j (t) = An( ) (t) i (1.5) pro každé = 1, 2, . . . , p a každé i = 1, 2, . . . , s, tedy m(t) =      m(1)(t) m(2)(t) ... m(p)(t)      =      An(1)(t) An(2)(t) ... An(p)(t)      =      A O . . . O O A . . . O ... ... ... ... O O . . . A           n(1)(t) n(2)(t) ... n(p)(t)      = I ⊗ A n(t); přitom O označuje nulovou a I jednotkovou čtvercovou matici řádu s. Označme dále di pravděpodobnost, že jedinec i-tého stadia opustí svou lokalitu a κ j pravděpodobnost, že jedinec, který opustil j-tou lokalitu se do konce projekčního intervalu dostane na lokalitu -tou, = j; při tomto označení musí platit p =1 =j κ j ≤ 1 (tento součet vyjadřuje pravděpodobnost, že jedinec, který opustil j-tou lokalitu, migraci přežije a dostane se na nějakou jinou lokalitu). Pravděpodobnost, že jedinec i-tého stadia opustí j-tou lokalitu a skončí na -té je tedy podle předpokladu 4. rovna diκ j. Položme κii = −1 (s tímto označením je p i=1 κij ≤ 0) a K =      κ11 κ12 . . . κ1p κ21 κ22 . . . κ2p ... ... ... ... κp1 κp2 . . . κpp      =      −1 κ12 . . . κ1p κ21 −1 . . . κ2p ... ... ... ... κp1 κp2 . . . −1      , D =      d1 0 . . . 0 0 d2 . . . 0 ... ... ... ... 0 0 . . . ds      . Po „demografické události dojde během projekčního intervalu k „dispersní události . Velikost n ( ) i subpopulace tvořené jedinci i-tého stadia na -té lokalitě na konci projekčního intervalu (neboli na začátku následujícího projekčního intervalu) před další „demografickou událostí bude sestávat z těch jedinců i-tého stadia, kteří na -té lokalitě byli na začátku uvažovaného projekčního intervalu a neemigrovali z ní (střední velikost takové subpopulace je (1−di)m ( ) i ), a z těch jedinců, kteří na -tou lokalitu během projekčního intervalu imigrovali z ostatních lokalit (střední velikost takové subpopulace je p j=1 j= dik jm (j) i ). Tedy s využitím (1.5) 1.2. MODELY DISPERSE 27 dostaneme n ( ) i (t + 1) = (1 − di)m ( ) i (t) + p j=1 j= diκ jm (j) i (t) = di p j=1 κ jm (j) i (t) + m ( ) i (t) = = p j=1 κ jdi An(j) (t) i + An( ) (t) i = p j=1 κ j DAn(j) (t) i + An( ) (t) i = =   p j=1 κ jDAn(j) (t) + An( ) (t)   i . To znamená, že n( ) (t + 1) = κ 1DAn(1) (t) + κ 2DAn(2) (t) + · · · + κ pDAn(p) (t) + An( ) (t) = = κ 1DA κ 2DA · · · κ pDA      n(1)(t) n(2)(t) ... n(p)(t)      + An( ) (t) = (K ⊗ DA)n(t) + An( ) (t) a dále n(t + 1) = (K ⊗ DA)n(t) +      A O . . . O O A . . . O ... ... ... ... O O . . . A           n(1)(t) n(2)(t) ... n(p)(t)      = (K ⊗ DA)n(t) + (I ⊗ A)n(t). Odtud vidíme, že model difúze populace lze zapsat ve tvaru n(t + 1) = (K ⊗ DA + I ⊗ A)n(t) nebo podrobněji      n(1) n(2) ... n(p)      (t + 1) = (K ⊗ DA + I ⊗ A)      n(1) n(2) ... n(p)      (t). (1.6) 1.2.2 Obecnější model difúze Pro popis struktury populace nyní zvolíme druhou z možností (1.4). Nebudeme požadovat splnění zjednodušujících předpokladů z oddílu 1.2.1 a pohyb jedinců mezi lokalitami budeme popisovat podrobněji. Lze totiž předpokládat, že migrace je proces rychlejší než „demografie . Budeme si tedy představovat, že během projekčního intervalu dojde k více „migračním událostem — přesunům z jedné lokality na jinou; počet „migračních událostí během projekčního intervalu označíme r. Budeme předpokládat, že 1 r ε. 28 KAPITOLA 1. KONSTRUKCE MODELŮ Schematicky můžeme nyní znázornit vývoj populace na -té lokalitě pro r = 3 obrázkem: Ec c cc c c c c c c c t0 1 2 n ( ) i (0) m ( ) i (0) m ( ) i (1 3 ) m ( ) i (2 3) n ( ) i (1) m ( ) i (1) m ( ) i (4 3) m ( ) i (5 3) n ( ) i (2) m ( ) i (2) m ( ) i (7 3 ) · · · Opuštění předpokladů 2.–4. vede k uvažování pravděpodobnosti, že jedinec i-tého stadia opustí j-tou lokalitu a během časového intervalu délky 1/r se dostane na lokalitu -tou. Označme tuto pravděpodobnost c ( j) i . Hodnota c ( ) i nyní vyjadřuje pravděpodobnost přežití a setrvání jedinců i-tého stadia na -té lokalitě po „demografické události . (Ve zjednodušené situaci z oddílu 1.2.1 je r = 1, c ( ) i = 1 − di a c ( j) i = dik j pro j = .) Střední množství jedinců i-tého stadia na -té lokalitě po jedné „migrační události tedy bude m ( ) i t + 1 r = p j=1 c ( j) i m (j) i (t); (1.7) čas t označuje levý krajní bod projekčního intervalu, i = 1, 2, . . . , s, = 1, 2, . . . , p. Položme Ci =       c (11) i c (12) i . . . c (1p) i c (21) i c (22) i . . . c (2p) i ... ... ... ... c (p1) i c (p2) i . . . c (pp) i       , C =      C1 O . . . O O C2 . . . O ... ... ... ... O O . . . Cs      ; nyní O označuje čtvercovou nulovou matici řádu p. Rovnosti (1.7) můžeme také přepsat ma- ticově: mi t + 1 r =       m (1) i m (2) i ... m (p) i       t + 1 r = Ci       m (1) i m (2) i ... m (p) i       (t) = Cimi(t), i = 1, 2, . . . , s, celkem tedy m t + 1 r =      m1 m2 ... ms      t + 1 r = C      m1 m2 ... ms      (t) = Cm(t). Strukturu populace po q + 1 „migračních událostech lze analogicky vyjádřit ve tvaru m t + q + 1 r = Cm t + q r , q = 0, 1, 2, . . . , r − 2. (1.8) Stejnou úvahou dostaneme, že struktura populace před další „demografickou událostí (tj. na konci projekčního intervalu) je n(t + 1) = Cm t + r − 1 r . 1.2. MODELY DISPERSE 29 S využitím (1.8) odtud dostaneme n(t + 1) = Cm t + r − 1 r = CCm t + r − 2 r = C2 m t + r − 2 r = · · · = Cr m(t). (1.9) Bez předpokladu 1. bude „demografické události na každé lokalitě popisovat jiná matice. Označme proto A( ) = a ( ) ij s i,j=1 , = 1, 2, . . . , p matici popisující rození a přežívání na -té lokalitě. Stejnou úvahou jako v případě rovnosti (1.5) odvodíme, že m ( ) i (t) = s j=1 a ( ) ij n ( ) j (t). (1.10) Tuto rovnost můžeme přepsat v maticovém tvaru                                    m (1) 1 m (2) 1 .. . m (p) 1 m (1) 2 m (2) 2 ... m (p) 2 .. . ... m (1) s m (2) s . .. m (p) s                                    (t) =                                    a (1) 11 0 . . . 0 a (1) 12 0 . . . 0 . . . . . . a (1) 1s 0 . . . 0 0 a (2) 11 . . . 0 0 a (1) 12 . . . 0 . . . . . . 0 a (1) 1s . . . 0 .. . .. . ... .. . .. . .. . ... .. . .. . .. . ... .. . 0 0 . . . a (p) 11 0 0 . . . a (p) 12 . . . . . . 0 0 . . . a (p) 1s a (1) 21 0 . . . 0 a (1) 22 0 . . . 0 . . . . . . a (1) 2s 0 . . . 0 0 a (2) 21 . . . 0 0 a (1) 22 . . . 0 . . . . . . 0 a (1) 2s . . . 0 ... ... ... ... ... ... ... ... ... ... ... ... 0 0 . . . a (p) 21 0 0 . . . a (p) 22 . . . . . . 0 0 . . . a (p) 2s .. . .. . .. . .. . .. . .. . .. . .. . .. . ... ... ... ... ... ... ... ... ... a (1) s1 0 . . . 0 a (1) s2 0 . . . 0 . . . . . . a (1) ss 0 . . . 0 0 a (2) s1 . . . 0 0 a (1) s2 . . . 0 . . . . . . 0 a (1) ss . . . 0 . .. . .. ... . .. . .. . .. ... . .. . .. . .. ... . .. 0 0 . . . a (p) s1 0 0 . . . a (p) s2 . . . . . . 0 0 . . . a (p) ss                                                                       n (1) 1 n (2) 1 .. . n (p) 1 n (1) 2 n (2) 2 ... n (p) 2 .. . ... n (1) s n (2) s . .. n (p) s                                    (t). Označme nyní Bij =       a (1) ij 0 . . . 0 0 a (2) ij . . . 0 ... ... ... ... 0 0 . . . a (p) ij       , B =      B11 B12 . . . B1s B21 B22 . . . B2s ... ... ... ... Bs1 Bs2 . . . Bss      . Rovnosti (1.10) pro i = 1, 2, . . . , s, = 1, 2, . . . , p tedy můžeme zapsat ve tvaru      m1 m2 ... ms      (t) =      B11 B12 . . . B1s B21 B22 . . . B2s ... ... ... ... Bs1 Bs2 . . . Bss           n1 n2 ... ns      (t) nebo stručně m(t) = Bn(t). Odtud a z rovnosti (1.10) nyní dostaneme model difúze n(t + 1) = Cr Bn(t) 30 KAPITOLA 1. KONSTRUKCE MODELŮ nebo podrobněji      n1 n2 ... ns      (t + 1) = Cr B      n1 n2 ... ns      (t). (1.11) 1.2.3 Příklad Uvažujme metapopulaci na dvou lokalitách strukturovanou do tří věkových tříd, tj. s = 3, p = 2. Obě lokality považujeme za stejně kvalitní, tedy specifické plodnosti i pravděpodobnosti přežití jsou na obou lokalitách stejné. Nechť plodní jsou jedinci druhé a třetí věkové třídy se specifickými fertilitami f2 a f3. Pravděpodobnost, že jedinci první, resp. druhé, věkové třídy přežijí projekční interval označíme p1, resp. p2. Jedinci třetí věkové třídy uhynou. O době migrace budeme předpokládat, že je stejná jako délka projekčního intervalu. Novorozenci nemigrují a pravděpodobnost opuštění lokality závisí pouze na věkové třídě. Náročnost cesty z první lokality na druhou může být jiná než cesty naopak; může jít např. o migraci vodních organismů proti proudu a po proudu. Pro jedince z různých věkových tříd se však neliší. Za těchto předpokladů můžeme vývoj uvažované metapopulace popisovat oběma uvedenými způsoby. V modelu popsaném v pododdílu 1.2.1 bude d1 = 0, takže A =   0 f2 f3 p1 0 0 0 p2 0   , D =   0 0 0 0 d2 0 0 0 d3   , K = −1 κ12 κ21 −1 . Projekční matice je tedy tvaru K ⊗ DA + I ⊗ A = −1 κ12 κ21 −1 ⊗   0 0 0 d2p1 0 0 0 d3p2 0   + 1 0 0 1 ⊗   0 f2 f3 p1 0 0 0 p2 0   = =         0 0 0 0 0 0 −d2p1 0 0 κ12d2p1 0 0 0 −d3p2 0 0 κ12d3p2 0 0 0 0 0 0 0 κ21d2p1 0 0 −d2p1 0 0 0 κ21d3p2 0 0 −d3p2 0         +         0 f2 f3 0 0 0 p1 0 0 0 0 0 0 p2 0 0 0 0 0 0 0 0 f2 f3 0 0 0 p1 0 0 0 0 0 0 p2 0         = =         0 f2 f3 0 0 0 (1 − d1)p1 0 0 κ12d2p1 0 0 0 (1 − d3)p2 0 0 κ12d3p2 0 0 0 0 0 f2 f3 κ21d2p1 0 0 (1 − d2)p1 0 0 0 κ21d3p2 0 0 (1 − d3)p2 0         . Při označení ˜A =   0 f2 f3 (1 − d2)p1 0 0 0 (1 − d3)p2 0   , 1.2. MODELY DISPERSE 31 M12 =   0 0 0 κ12d2p1 0 0 0 κ12d3p2 0   , M21 =   0 0 0 κ21d2p1 0 0 0 κ21d3p2 0   můžeme model (1.11) zapsat ve tvaru n(1) n(2) (t + 1) = ˜A M12 M21 ˜A n(1) n(2) (t); matice ˜A popisuje plodnosti a přežívání nemigrujících jedinců, matice M12 a M21 popisují migrace. V modelu popsaném v pododdílu 1.2.2 bude r = 1, A(1) = A(2) = A, c21 1 = c12 1 = 0, c11 1 = c22 1 = 1, c12 2 = d2κ12, c21 2 = d2κ21, c11 2 = c22 2 = 1 − d2, c12 3 = d3κ12, c21 3 = d3κ21, c11 3 = c22 3 = 1 − d3, takže B11 = 0 0 0 0 , B12 = f2 0 0 f2 , B13 = f3 0 0 f3 , B21 = p1 0 0 p1 , B22 = 0 0 0 0 , B23 = 0 0 0 0 , B31 = 0 0 0 0 , B32 = p2 0 0 p2 , B33 = 0 0 0 0 , C1 = 1 0 0 1 , C2 = 1 − d2 d2κ12 d2κ21 1 − d2 , C3 = 1 − d3 d3κ12 κ21d3 1 − d3 a projekční matice je tvaru         1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 − d2 d2κ12 0 0 0 0 d2κ21 1 − d2 0 0 0 0 0 0 1 − d3 d3κ12 0 0 0 0 d3κ21 1 − d3                 0 0 f2 0 f3 0 0 0 0 f2 0 f3 p1 0 0 0 0 0 0 p1 0 0 0 0 0 0 p2 0 0 0 0 0 0 p2 0 0         = =         0 0 f2 0 f3 0 0 0 0 f2 0 f3 (1 − d2)p1 d2κ12p1 0 0 0 0 d2κ21p1 (1 − d2)p1 0 0 0 0 0 0 (1 − d3)p2 d3κ12p2 0 0 0 0 d3κ21p2 (1 − d3)p2 0 0         . Při označení F2 = f2 0 0 f2 , F3 = f3 0 0 f3 , P1 = (1 − d2)p1 d2κ12p1 d2κ21p1 (1 − d2)p1 , P2 = (1 − d3)p2 d3κ12p1 d3κ21p2 (1 − d3)p2 32 KAPITOLA 1. KONSTRUKCE MODELŮ můžeme model (1.11) zapsat ve tvaru   n1 n2 n3   (t + 1) =   O F2 F3 P1 O O O P2 O     n1 n2 n3   (t). Matice F2, resp. F3, vyjadřuje plodnosti druhé, resp. třetí, věkové třídy, matice P1, resp. P2, popisuje přežívání migrujících i nemigrujících jedinců druhé, resp. třetí, věkové třídy; (Pi) j je pravděpodobnost úspěšné migrace jedinců (i+1)-té věkové třídy migrujících z j-té lokality na -tou, tj. v případě = j pravděpodobnost přežití a setrvání na lokalitě. Projekční matice modelu je v tomto případě blokově Leslieho typu. Pokud budeme předpokládat, že náročnost cesty mezi lokalitami nezávisí na směru, tj. κ12 = κ21 = κ, dostaneme M12 = M21 = M = κ   0 0 0 d2p1 0 0 0 d3p2 0   První model můžeme zapsat ve tvaru n(1) n(2) (t + 1) = 0 1 1 0 ⊗ M + 1 0 0 1 ⊗ ˜A n(1) n(2) (t) a druhý   n1 n2 n3   (t + 1) = M ⊗ 0 1 1 0 + ˜A ⊗ 1 0 0 1   n1 n2 n3   (t). Kapitola 2 Modely s konstantní projekční maticí 2.1 Příklad — populace strukturovaná podle plodnosti Maticový populační model n(t + 1) = An(t) je vlastně vektorová lineární diferenční rovnice neboli systém lineárních autonomních diferenčních rovnic. Její řešení ukážeme nejprve na jednoduchém příkladu — na modelu populace rozdělené na juvenilní a plodné jedince, tj. na rovnici (9), nebo rozepsané do složek (7) a (8). Místo podmínek (6) budeme uvažovat mírně slabší podmínky 0 < σ1 ≤ 1, 0 ≤ σ2 ≤ 1, 0 < γ ≤ 1, 0 < ϕ, (2.1) abychom z úvah nevyloučili „klasické Fibonacciovy králíky . Z rovnice (7) vyjádříme n2(t) = 1 ϕ n1(t + 1) − σ1(1 − γ)n1(t) (2.2) a dosadíme do (8) n2(t + 1) = σ1γn1(t) + σ2 ϕ n1(t + 1) − σ1(1 − γ)n1(t) . V rovnici (7) dosadíme t + 1 za t a dále do ní dosadíme z předchozí rovnosti: n1(t + 2) = σ1(1 − γ)n1(t + 1) + ϕn2(t + 1) = = σ1(1 − γ)n1(t + 1) + ϕσ1γn1(t) + σ2 n1(t + 1) − σ1(1 − γ)n1(t) = = σ1(1 − γ) + σ2 n1(t + 1) + ϕσ1γ − (1 − γ)σ1σ2 n1(t). První složka řešení systému (9) diferenčních rovnic prvního řádu je tedy řešením lineární diferenční rovnice druhého řádu n1(t + 2) = σ1(1 − γ) + σ2 n1(t + 1) + ϕσ1γ − (1 − γ)σ1σ2 n1(t). (2.3) Její charakteristickou rovnicí je kvadratická rovnice λ2 − σ1(1 − γ) + σ2 λ + (1 − γ)σ1σ2 − ϕσ1γ = 0, (2.4) 33 34 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ která má diskriminant D = σ1(1 − γ) + σ2 2 − 4 (1 − γ)σ1σ2 − ϕσ1γ = σ1(1 − γ) − σ2 2 + 4ϕσ1γ. (2.5) Vzhledem k předpokladům (2.1) je D > σ1(1 − γ) − σ2 2 ≥ 0, takže charakteristická rovnice (2.4) má dva reálné různé kořeny λ1 = σ1(1 − γ) + σ2 + √ D 2 , λ2 = σ1(1 − γ) + σ2 − √ D 2 . (2.6) Kořeny λ1 a λ2 zřejmě splňují nerovnosti λ1 > 0, λ1 ≥ |λ2|, λ1 − λ2 > 0. (2.7) Poznamenejme, že rovnost λ1 = |λ2|, tj. λ1 = −λ2 nastane právě tehdy, když σ1(1−γ) = −σ2, tedy vzhledem k (2.1) právě tehdy, když σ2 = 0 a γ = 1. Lineární diferenční rovnice druhého řádu (rekurentní formule) (2.3) má obecné řešení n1(t) = aλt 1 + bλt 2. Konstanty a, b získáme z počátečních podmínek. Předpokládejme, že známe počáteční hodnoty n1(0) a n2(0). Velikost složek populace nemůže být záporná a celková velikost existující populace je kladná, platí n1(0) ≥ 0, n2(0) ≥ 0, n1(0) + n2(0) > 0. (2.8) Z rovnosti (2.2) dostaneme n1(1) = σ1(1 − γ)n1(0) + ϕn2(0). Známe tedy hodnoty n1(0) a n1(1), které musí splňovat rovnosti n1(0) = a + b n1(1) = aλ1 + bλ2. Řešením tohoto systému rovnic pro neznámé parametry a, b je a = n1(1) − λ2n1(0) λ1 − λ2 , b = λ1n1(0) − n1(1) λ1 − λ2 , takže n1(t) = n1(1) − λ2n1(0) λ1 − λ2 λt 1 + λ1n1(0) − n1(1) λ1 − λ2 λt 2. Druhou složku řešení systému (9) dostaneme dosazením vypočítané první složky do rovnosti (2.2): n2(t) = n1(1) − λ2n1(0) λ1 − σ1(1 − γ) ϕ(λ1 − λ2) λt 1 + λ1n1(0) − n1(1) λ2 − σ1(1 − γ) ϕ(λ1 − λ2) λt 2. Řešení systému (9) tedy je n1(t) = σ1(1 − γ)n1(0) + ϕn2(0) − λ2n1(0) λ1 − λ2 λt 1 + λ1n1(0) − σ1(1 − γ)n1(0) − ϕn2(0) λ1 − λ2 λt 2, n2(t) = σ1(1 − γ)n1(0) + ϕn2(0) − λ2n1(0) λ1 − σ1(1 − γ) ϕ(λ1 − λ2) λt 1+ 2.1. PŘÍKLAD — POPULACE STRUKTUROVANÁ PODLE PLODNOSTI 35 + λ1n1(0) − σ1(1 − γ)n1(0) − ϕn2(0) λ2 − σ1(1 − γ) ϕ(λ1 − λ2) λt 2, kde λ1 a λ2 jsou dány rovnostmi (2.5) a (2.6). Řešení systému (9) lze také stručně zapsat ve tvaru n(t) = αλt 1w(1) + βλt 2w(2) = λt 1 αw(1) + β λ2 λ1 t w(2) , (2.9) kde α = σ1(1 − γ) − λ2 n1(0) + ϕn2(0) λ1 − λ2 , β = λ1 − σ1(1 − γ) n1(0) − ϕn2(0) λ1 − λ2 , w(1) =   w (1) 1 w (1) 2   = 1 ϕ ϕ λ1 − σ1(1 − γ) , w(2) =   w (2) 1 w (2) 2   = 1 ϕ ϕ λ2 − σ1(1 − γ) . Přímým výpočtem můžeme ověřit, že λ1, λ2 jsou vlastními hodnotami matice A a vektory w(1), w(2) jsou příslušné vlastní vektory. Z rovností (2.5), (2.6) a nerovností (2.1) dostaneme λ1 − σ1(1 − γ) = −σ1(1 − γ) + σ2 + σ2 − σ1(1 − γ) 2 + 4ϕσ1γ 2 > 0, σ1(1 − γ) − λ2 = σ1(1 − γ) − σ2 + σ1(1 − γ) − σ2 2 + 4ϕσ1γ 2 > 0. Z těchto nerovností plyne w (1) 1 > 0, w (1) 2 > 0, w (2) 1 > 0, w (2) 2 < 0. (2.10) Ze druhé z nich spolu s nerovnostmi (2.7) a (2.8) plyne α > 0. (2.11) Z vyjádření řešení (2.9) dostaneme 1 λt 1 n(t) − αw(1) = β λ2 λ1 t w(2) . (2.12) Označme dále q(t) = n1(t) n2(t) = αw (1) 1 λt 1 + βw (2) 1 λt 2 αw (1) 2 λt 1 + βw (2) 2 λt 2 = w (1) 1 + β α w (2) 1 λ2 λ1 t w (1) 2 + β α w (2) 2 λ2 λ1 t poměr velikostí složek populace (daných rovností (2.9)) v čase t. Nerovnosti (2.7), (2.10) a (2.11) ukazují, že veličina q(t) je definována korektně. Nechť λ1 = |λ2|, tj. σ2 > 0 nebo γ = 1. Podle nerovností (2.7) je λ1 > |λ2| a tedy lim t→∞ λ2 λ1 t = 0. 36 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ V tomto případě je lim t→∞ 1 λt 1 n(t) − αw(1) = 0, (2.13) tj. funkce n(t) a αλt 1w(1) jsou asymptoticky ekvivalentní, a lim t→∞ q(t) = w (1) 1 w (1) 2 . Pokud navíc λ2 = 0, což podle (2.5) a (2.6) nastane právě tehdy, když σ2(1 − γ) = ϕγ, pak n(t) = αλt 1w(1) pro t = 1, 2, 3, . . . (2.14) a q(t) = w (1) 1 w (1) 2 pro t = 1, 2, 3, . . . . Je-li λ1 = |λ2|, tj. λ1 = −λ2 podle (2.7), pak λ1 λ2 = −1, takže q(t) = w (1) 1 + β α w (2) 1 (−1)t w (1) 2 + β α w (2) 2 (−1)t , q(t + 2) = q(t). Dále 1 λt 1 n(t) − αw(1) = β(−1)t w(2) . Pro velikost vektoru na pravé straně této rovnosti vzhledem k (2.10) platí 1 λt 1 n(t) − αw(1) 1 = β w(2) 1 = β w (2) 1 − w (2) 2 . (2.15) Velikost vektoru n(t) řešení rovnice (9) tedy v každém případě podle rovností (2.13), (2.14) a (2.15) splňuje asymptotickou rovnost n(t) 1 = O(λt 1); (2.16) velikost populace se „po dostatečně dlouhém vývoji chová jako geometrická posloupnost s kvocientem λ1 . Dosud provedené výpočty můžeme shrnout: • Matice A v rovnici (9) má dvě reálné různé vlastní hodnoty λ1, λ2 takové, že λ1 > 0, λ1 ≥ |λ2|. • Vlastní vektor w(1) příslušný k vlastní hodnotě λ1 má obě složky kladné. • Řešení rovnice (9) je dáno formulí (2.9). Přitom w(1) a w(2) jsou vlastní vektory příslušné k vlastním hodnotám λ1 a λ2 matice A, parametry α > 0 a β závisí na počátečních podmínkách n1(0), n2(0). 2.2. PERRONOVA-FROBENIOVA TEORIE 37 • Řešení n(t) rovnice (9) je asymptoticky ekvivalentní s geometrickou posloupností s kvocientem λ1. • Pokud λ1 > |λ2|, pak poměr složek vektoru řešení n(t) konverguje k poměru složek vlastního vektoru w(1) příslušného k vlastní hodnotě λ1. Pokud λ1 = −λ2, pak poměr složek vektoru řešení n(t) se periodicky mění, perioda je rovna 2. Kladná vlastní hodnota λ1 matice A (dominantní vlastní hodnota) tedy představuje růstový koeficient populace. V případě populace iteroparní (σ2 > 0) nebo populace se zpožděným dospíváním (γ < 1) se poměr velikostí jednotlivých tříd v průběhu vývoje ustálí; složky normovaného vlastního vektoru příslušného k dominantní vlastní hodnotě, tj. vektoru 1 w(1) 1 w(1) = 1 w (1) 1 + w (1) 2 w(1) představují relativní zastoupení jednotlivých tříd, tedy stabilizovanou strukturu populace. Výsledky lze také ilustrovat několika konkrétními případy. Za jednotku času (délku projekčního intervalu) zvolíme dobu potřebnou k „vyprodukování jednoho potomka. Bude tedy ϕ = 1. Za počáteční hodnoty zvolíme n1(0) = 0, n2(0) = 1, tedy stejně jako v případě Fibonacciových králíků začínáme s jedním plodným párem. Platí tedy α = 1 λ1 − λ2 = −β a řešení je tvaru n1(t) = 1 λ1 − λ2 λt 1 − λt 2 , n2(t) = 1 λ1 − λ2 λt 1 λ1 − σ1(1 − γ) − λt 2 λ2 − σ1(1 − γ) , kde λ1,2 = σ1(1 − γ) + σ2 ± √ D 2 , D = σ1(1 − γ) − σ2 2 + 4σ1γ. Výsledky pro několik zvolených hodnot parametrů jsou shrnuty v tabulce 2.1. Vidíme, že celková velikost populace může neomezeně růst, klesat k nule (populace vymírá), konvergovat k nějaké hodnotě, případně této hodnoty bezprostředně dosáhnout. V případě, že populace není semelparní s bezprostředním dospíváním, struktura populace (relativní zastoupení jednotlivých složek) konverguje k nějaké hodnotě; k této hodnotě struktura konverguje monotonně nebo s tlumenými oscilacemi, případně jí dosáhne hned v prvním časovém kroku. 2.2 Perronova-Frobeniova teorie Všechny matice v tomto oddílu budou typu n×n, všechny vektory budou n-rozměrné. Symbol |A|, resp. |v|, bude označovat matici, jejíž složky jsou (|A|)ij = |aij|, resp. vektor, jehož složky jsou (|v|)i = |vi|. Dále budeme zapisovat A ≥ c . . . (∀i, j)aij ≥ c, tj. (∀i, j)(A)ij ≥ c, v ≥ c . . . (∀i)vi ≥ c, tj. (∀i)(v)i ≥ c, A ≥ B . . . (∀i, j)aij ≥ bij, tj. (∀i, j)(A)ij ≥ (B)ij, v ≥ w . . . (∀i)vi ≥ wi, tj. (∀i)(v)i ≥ (w)i 38KAPITOLA2.MODELYSKONSTANTNÍPROJEKČNÍMATICÍ σ1 σ2 γ A λ1 λ2 n(t) lim t→∞ n(t) 1 lim t→∞ n1(t) n2(t) 1 1 1 0 1 1 1 1+ √ 5 2 1− √ 5 2 n1(t) = √ 5 5 1+ √ 5 2 t − 1− √ 5 2 t n2(t) = √ 5 5 1+ √ 5 2 t+1 − 1− √ 5 2 t+1 ∞ √ 5−1 2 0 1 2 3 4 5 6 04812 abundances t 0 1 2 3 4 5 6 0.00.40.8 structure t 8 9 2 3 1 0 1 8 9 2 3 4 3 2 3 n1(t) = 2 3 t−1 2t − 1) n2(t) = 2 3 t 2t+1 − 1) ∞ 3 4 0 1 2 3 4 5 6 0123 abundances t 0 1 2 3 4 5 6 0.00.40.8 structure t 7 9 2 3 1 7 2 3 1 1 9 2 3 1 1 3 n1(t) = 3 2 1 − 1 3 t n2(t) = 1 2 1 + 1 3 t 2 3 0 1 2 3 4 5 6 0.01.0 abundances t 0 1 2 3 4 5 6 0.00.40.8 structure t 3 4 1 2 1 3 1 2 1 1 4 1 2 1 0 n1(t) =    0, t = 0 1, t ≥ 1 n2(t) =    1, t = 0 1 2 , t ≥ 1 3 2 2 0 1 2 3 4 5 6 0.00.40.8 abundances t 0 1 2 3 4 5 6 0.00.40.8 structure t 8 9 1 9 1 0 1 8 9 1 9 1 − 8 9 n1(t) = 9 17 1 − (− 8 9 )t n2(t) = 9 17 1 − (− 8 9 )t+1 18 17 1 0 1 2 3 4 5 6 0.00.40.8 abundances t 0 1 2 3 4 5 6 0.00.40.8 structure t 1 0 1 0 1 1 0 1 −1 n1(t) = 1 2 1 − (−1)t n2(t) = 1 2 1 + (−1)t 1 neexistuje 0 1 2 3 4 5 6 0.00.40.8 abundances t 0 1 2 3 4 5 6 0.00.40.8 structure t 8 9 0 1 0 1 8 9 0 2 √ 2 3 − 2 √ 2 3 n1(t) = √ 2 4 2 √ 2 3 t 1 − (−1)t n2(t) = √ 2 4 2 √ 2 3 t+1 1 + (−1)t 0 neexistuje 0 1 2 3 4 5 6 0.00.40.8 abundances t 0 1 2 3 4 5 6 0.00.40.8 structure t Tabulka2.1:Speciálnípřípadymodelu(9)ajejichřešení.Vevšechmodelechjsoupočáteční podmínkyn1(0)=0,n2(0)=1aplodnostϕ=1.Grafnalevozobrazujeprůběhvelikostí složekpopulace;velikostskupinyjuvenilníchjedincůn1jevyznačenazeleně,velikostsku- pinyplodnýchjedincůn2jevyznačenačerveně.Grafnapravozobrazujevývojrelativního zastoupeníjednotlivýchsložekpopulace;juvenilnízeleně,plodnáčerveně. 2.2. PERRONOVA-FROBENIOVA TEORIE 39 a podobně. Symbol I bude označovat jednotkovou matici. Pro matici A a vektor v dále klademe ker A = {w : Aw = 0} , v = n i=1 v2 i ; ker A je zřejmě vektorový prostor dimenze nejvýše n, tj. dim ker A ≤ n, v je eukleidovská norma vektoru v. Definice 1. Matice A se nazývá nezáporná, je-li A ≥ 0 a nazývá se kladná, je-li A > 0. Definice 2. Nezáporná matice A se nazývá – primitivní, pokud (∃k ∈ N)Ak > 0, – imprimitivní, pokud není primitivní, tj. (∀k ∈ N)(∃i, j) Ak ij = 0, – reducibilní, pokud (∃i, j)(∀k ∈ N) Ak ij = 0, – ireducibilní, pokud není reducibilní, tj. (∀i, j)(∃k ∈ N) Ak ij > 0. Poznámka 1. Přímo z definice plyne, že každá primitivní matice je ireducibilní a každá reducibilní matice je imprimitivní. Třídu nezáporných matic lze tedy rozložit na tři disjunktní části: matice reducibilní, matice primitivní a matice současně ireducibilní a imprimitivní. Tvrzení 1. Je-li A ≥ 0 a v ≥ w pak Av ≥ Aw. Je-li A > 0, v ≥ w a existuje i ∈ {1, 2, . . . , n} takové, že vi > wi pak Av > Aw. Důkaz. Plyne bezprostředně z vyjádření (Av)k = n j=1 akjvj, n j=1 akjwj = (Aw)k. Tvrzení 2. Je-li A > 0, v vlastní vektor matice A příslušný k vlastní hodnotě λ a v ≥ 0, pak v > 0 a λ > 0. Důkaz. Poněvadž v je vlastním vektorem, je v = 0 a tedy existuje i ∈ {1, 2, . . . , n} takový index, že vi > 0. Podle druhé části tvrzení 1 je λv = Av > A0 = 0. To znamená, že pro každý index j je λvj > 0. Zejména tedy λvi > 0, z čehož plyne, že λ > 0, neboť vi > 0. Dále pro libovolný index j je vj > 0, neboť λvj > 0. Tvrzení 3. Je-li A ≥ 0 primitivní a v ≥ 0 její vlastní vektor příslušný k vlastní hodnotě λ, pak v > 0, λ > 0. Důkaz. Z tvrzení 1 plyne, že λv = Av ≥ 0, takže λ ≥ 0. Poněvadž A je primitivní, existuje k ∈ N takové, že Ak > 0. Poněvadž Av = λv, je také Akv = Ak−1Av = Ak−1(λv) = λAk−1v = · · · = λkv. Tvrzení 2 nyní implikuje v > 0 a λk > 0, takže λ = 0. Tvrzení 4. Nechť matice A splňuje předpoklady (i) A ≥ 0, A = 0; 40 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ (ii) existuje číslo λ ∈ R a vektor u tak, že ATu = λu, u > 0 (vektor u je vlastní vektor matice AT příslušný k vlastní hodnotě λ, který má všechny složky kladné); (iii) existuje číslo µ ∈ R a vektor v tak, že Av = µv, v ≥ 0, v = 0 (vektor v je vlastní vektor matice A příslušný k vlastní hodnotě µ, který má všechny složky nezáporné a alespoň jednu kladnou). Pak µ = λ. Důkaz. Platí λuT v = AT u T v = uT Av = uT µv = µuT v. Z kladnosti vektoru u a z nezápornosti a nenulovosti vektoru v plyne uTv > 0. Výraz uTv lze tedy v poslední rovnosti vykrátit, takže λ = µ. Tvrzení 5. Nechť A > 0 splňuje předpoklady (ii) a (iii) tvrzení 4 (z nerovnosti A > 0 plyne i splnění předpokladu (i)) a symboly λ(= µ), v mají stejný význam jako v tvrzení 4. Je-li w vlastní vektor matice A příslušný k vlastní hodnotě λ, pak existuje číslo α ∈ R takové, že w = αv, tj. dim (ker(A − λI)) = 1. Důkaz. Buď v vektor z tvrzení 2. Pak je v > 0 a λ > 0. Položme α = min wj vj : j = 1, 2, . . . , n , i takový index, že α = wi vi . Pro každý index j tedy platí α = wi vi ≤ wj vj . Odtud plyne wj − αvj ≥ 0, wi − αvi = 0, (2.17) takže w − αv ≥ 0. Dále platí A(w − αv) = λw − αλv = λ(w − αv). To znamená, že vektor w − αv je buď vlastním vektorem příslušným k vlastní hodnotě λ, nebo platí w − αv = o. Nezáporný vlastní vektor je podle tvrzení 2 kladný a podle (2.17) má vektor w − αv alespoň jednu složku nulovou, nemůže tedy být vlastním vektorem. Nastává tedy druhá z vylučujících se možností, w − αv = o. Tvrzení 6. Nechť A ≥ 0, w je vlastní vektor příslušný k vlastní hodnotě λ. Pak A|w| i = n j=1 aij|wj| = n j=1 |aijwj| ≥ n j=1 aijwj = |(Aw)i| = |(λw)i| = |λ| |wi|, (2.18) tj. A|w| ≥ |λ| |w|. Důkaz. Nerovnost je trojúhelníková. Tvrzení 7. Nechť A ≥ 0. Pak množina SA = c ≥ 0 : ∃v(c) v(c) ≥ 0, v(c) = 1, Av(c) ≥ cv(c) je neprázdná a shora omezená. 2.2. PERRONOVA-FROBENIOVA TEORIE 41 Důkaz. Buď v(0) libovolný nezáporný vektor takový, že v(0) = 1. Podle tvrzení 1 je Av(0) ≥ A0 = 0 = 0v(0) , takže 0 ∈ SA, SA = ∅. Buď c ∈ SA a v(c) příslušný vektor, který existuje podle definice množiny SA. Nechť i je takový index, že v (c) i = max v (c) 1 , v (c) 2 , . . . , v (c) n . Pak je v (c) i > 0 a cv (c) i ≤ Av(c) i = n j=1 aijv (c) j ≤ n j=1 aijv (c) i ≤ v (c) i max    n j=1 alj : l = 1, 2, . . . , n    , tedy c ≤ max    n j=1 alj : l = 1, 2, . . . , n    a c je horní závora množiny SA. Tvrzení 8. Nechť A ≥ 0, SA je množina zavedená v tvrzení 7 a λ1 = sup SA. Pak pro každý vektor w platí A|w| ≤ λ1|w|. Důkaz. Nulový vektor splňuje uvedenou nerovnost triviálně. Připusťme, že existuje nenulový vektor w splňující nerovnost A|w| > λ1|w| a položme ε = min 1 |wi| (A|w|)i − λ1|wi| : |wi| > 0 . Pak je ε > 0 a ε|wi| ≤ (A|w|)i − λ1|wi| pro každý index i, (λ1 + ε) |wi| ≤ (A|w|)i , (λ1 + ε)|w| ≤ A|w|. Položíme-li v(λ1+ε) = 1 w |w|, dostaneme, že v(λ1+ε) ≥ 0, v(λ1+ε) = 1 a Av(λ1+ε) = 1 w A|w| ≥ 1 w (λ1 + ε)|w| = (λ1 + ε)v(λ1+ε) , takže λ1 + ε ∈ SA, což je ve sporu s definicí suprema. Tvrzení 9. Nechť A ≥ 0. Pro každou její vlastní hodnotu λ platí |λ| ≤ λ1 = sup SA. Důkaz. Buď λ vlastní hodnota matice A a w příslušný vlastní vektor. Podle tvrzení 6 a 8 je |λ||w| ≤ A|w| ≤ λ1|w|. Poněvadž w jakožto vlastní vektor je nenulový, existuje index i takový, že wi > 0. Z předchozí nerovnosti nyní dostaneme |λ|wi ≤ λ1wi a z toho dále plyne, že |λ| ≤ λ1. 42 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Tvrzení 10. Nechť A ≥ 0 a λ1 = sup SA. Pak λ1 ≥ 0, λ1 je vlastní hodnotou matice A a příslušný vlastní vektor v ≥ 0. Důkaz. Nejprve ukážeme, že množina M = {v : v ≥ 0, v = 1} je kompaktní: Z trojúhelníkové nerovnosti pro normu plyne, že pro vektory v(1), v(2) ∈ M platí v(1) − v(2) ≤ v(1) + v(2) = 1 + 1 = 2, takže množina M je ohraničená. Buď w(k) ∞ k=1 ⊆ M posloupnost vektorů konvergující k vektoru v v prostoru Rn s metrikou určenou euklidovskou normou, tj. pro každý index i platí lim k→∞ n i=1 w (k) i − vi 2 = 0, neboli lim k→∞ w (k) i = vi. Poněvadž w (k) i ≥ 0, je také vi ≥ 0, tj. v ≥ 0. Z toho, že zobrazení F : Rn → R dané předpisem F(u) = u je spojité, plyne podle Heineovy podmínky F lim k→∞ w(k) = lim k→∞ F w(k) , tj. v = lim k→∞ w(k) = 1. Celkem tedy dostáváme, že v ∈ M. Množina M s konvergentní posloupností obsahuje i její limitu, takže tato množina je také uzavřená. Hodnota λ1 = sup SA je limitou posloupnosti čísel z množiny SA, tj. existuje posloupnost {ck}∞ k=1 ⊆ SA taková, že lim k→∞ ck = λ1. K číslům ck ∈ SA existují vektory v(ck) takové, že v(ck) ≥ 0, v(ck) = 1 (2.19) a Av(ck) ≥ ckv(ck) . (2.20) Relace (2.19) říkají, že všechny vektory v(ck) jsou prvky množiny M. Z její kompaktnosti plyne, že existuje posloupnost v(ckl ) ∞ l=1 vybraná z posloupnosti vektorů v(ck) ∞ k=1 taková, že lim l→∞ v(ckl ) = v ∈ M. Z první relace (2.19) dále plyne v ≥ 0, tj. |v| = v. Z (2.20) plyne Av(ckl ) ≥ ckl v(ckl ) . Poněvadž lineární zobrazení je spojité, dostaneme limitním přechodem l → ∞ z poslední nerovnosti nerovnost Av ≥ λ1v. Z ní s využitím tvrzení 8 dostaneme Av = λ1v, což znamená, že v je vlastní vektor příslušný k vlastní hodnotě λ1. Tvrzení 11. Nechť A ≥ 0 je primitivní. Pak existuje vlastní hodnota λ1 > 0 matice A taková, že příslušný vlastní vektor v > 0, dim ker(A−λ1I) = 1 a pro každou vlastní hodnotu λ = λ1 matice A platí λ1 > |λ|. 2.2. PERRONOVA-FROBENIOVA TEORIE 43 Důkaz. Položíme λ1 = sup SA, kde SA je množina zavedená v tvrzení 7. Podle tvrzení 10 je λ1 vlastní hodnotou matice A a příslušný vlastní vektor v ≥ 0. Podle tvrzení 3 je λ1 > 0 a v > 0. Matice AT je také primitivní. Stejnou úvahou ukážeme, že existuje λ > 0 vlastní hodnota matice AT a příslušný vlastní vektor u > 0. Z tvrzení 4 dostaneme rovnost λ = λ1. Poněvadž matice A je primitivní, existuje k ∈ N takové, že Ak > 0. Úvahy lze zopakovat pro matici Ak a její vlastní hodnoty λk 1. Tím se ukáže, že matice Ak splňuje předpoklady tvrzení 5. Jsou-li nyní v(1) a v(2) dva vlastní vektory matice A příslušné k vlastní hodnotě λ1, platí Ak v(1) = λk 1v(1) , Ak v(2) = λk 1v(2) , takže podle tvrzení 5 je vektor v(2) násobkem vektoru v(1), tj. dim ker(A − λ1I) = 1. Podle tvrzení 9 nemá matice A vlastní hodnoty s absolutní hodnotou větší než λ1. Buď λ vlastní hodnota matice A taková, že |λ| = λ1 a w příslušný vlastní vektor. Z tvrzení 6 dostaneme A|w| ≥ |λ| |w| = λ1|w|, z čehož podle tvrzení 8 plyne A|w| = λ1|w|. (2.21) To znamená, že |w| ∈ ker A − λ1I , takže podle již dokázaného, je vektor |w| násobkem vektoru v a poněvadž v > 0, je také |w| > 0. (2.22) Dále pro libovolný index i platí λ1|wi| = A|w| i = n j=1 aij|wj| ≥ n j=1 aijwj = Aw i = |λwi| = |λ| |wi| = λ1|wi|. V trojúhelníkové nerovnosti tedy nastává rovnost, což znamená, že argumenty všech sčítanců jsou stejné, arg aij|wj| = arg aijwj pro všechny indexy j. Protože aij ∈ R, aij ≥ 0, tj. arg aij|wj| = 0, je také arg aijwj = 0, tj. wj ∈ R a wj ≥ 0. Dále |wj| = wj, |w| = w. Vzhledem k (2.22) je wj > 0. Nyní s využitím (2.21) dostaneme λwj = λw j = Aw j = A|w| j = λ1|w| j = λ1w j = λ1wj, takže λ = λ1. Tvrzení 12. Nechť A ≥ 0, λ1, λ jsou její vlastní hodnoty takové, že λ1 = sup SA a |λ| = λ1, arg λ = ϕ, tj. λ = eiϕλ1. Pak existují čísla ϑ1, ϑ2, . . . , ϑn ∈ R, že A = eiϕDAD−1, kde D = diag eiϑ1 , eiϑ2 , . . . , eiϑn . Důkaz. Nechť w je vlastní vektor příslušný k vlastní hodnotě λ, tj. Aw = λw. Podle tvrzení 6 a 8 je λ1|w| = |λ| |w| ≤ A|w| ≤ λ1|w|, takže A|w| = λ1|w|. (2.23) Položme ϑk =    Arg wk |wk| , wk = 0, 0, wk = 0, 44 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ tj. eiϑk = wk |wk| , pokud wk = 0. Pak eiϑk |wk| = wk pro každý index k, D|w| = w a dále AD|w| = Aw = λw = eiϕ λ1w = eiϕ λ1D|w|, tedy eiϕ λ1D|w| = AD|w|, eiϕ λ1|w| = D−1 AD|w| a s využitím (2.23) A|w| = e−iϕD−1AD|w|. Položme C = e−iϕD−1AD. Pak A|w| = C|w|. Poněvadž A|w| ≥ 0, je také C|w| ≥ 0, tedy C|w| = C|w| . Celkem s využitím trojúhelníkové nerovnosti dostaneme A|w| = C|w| = C|w| ≤ |C| |w| = A|w|. V trojúhelníkové nerovnosti nastává rovnost n j=1 clj|wj| = n j=1 |clj| |wj|, l = 1, 2, . . . , n, což znamená, že clj|wj| a |clj| |wj| mají stejné argumenty, tedy clj ∈ R, clj ≥ 0, |C| = C. Dále cij = e−iϕe−iϑi aijeiϑj , tj. |cij| = |aij| = aij, A = |C| = C = e−iϕD−1AD a odtud plyne tvrzení. Tvrzení 13. Nechť A ≥ 0, λ1 = sup SA, λ2, . . . , λd jsou všechny její různé vlastní hodnoty takové, že λ1 = |λ2| = |λ3| = · · · = |λd|, 0 = Arg λ1 < Arg λ2 < · · · < Arg λd. Pak λj = e2πi(j−1)/dλ1 a dim ker(A − λ1I) = dim ker(A − λjI) , j = 1, 2, . . . , d. Důkaz. Označme ϕj = Arg λj, tj. λj = eiϕj λ1, j = 1, 2, . . . , d. Nejprve si všimneme několika jednoduchých fakt: (i) Z tvrzení 12 plyne, že k vlastní hodnotě λj existuje regulární diagonální matice Dj taková, že A = eiϕj DjAD−1 j . (ii) Matice A a DjAD−1 j mají stejné vlastní hodnoty. Je-li totiž λ vlastní hodnotou matice A pak matice S = A − λI je singulární a tedy také matice DjSD−1 j = DjAD−1 j − λI je singulární, což znamená, že λ je také vlastní hodnotou matice DjAD−1 j . Podobně nahlédneme, že libovolná vlastní hodnota matice DjAD−1 j je také vlastní hodnotou matice A. (iii) λ je vlastní hodnotou matice A právě tehdy, když e−iϕj λ je vlastní hodnotou matice e−iϕj A, neboť matice A − λI a e−iϕj (A − λI) = e−iϕj A − e−iϕj λI jsou současně singulární nebo regulární. (iv) Je-li λj vlastní hodnotou matice A, pak λj je také vlastní hodnotou matice A, neboť matice A je reálná. Odtud dále plyne, že λj = λk pro nějaké k ∈ {1, 2, . . . , d}, tj. ϕj + ϕk = 2π, neboť λ1, λ2, . . . , λd jsou všechny vlastní hodnoty stejného modulu. 2.2. PERRONOVA-FROBENIOVA TEORIE 45 Je-li d = 1, je tvrzení triviální. Je-li d = 2, pak λ2 ∈ R — kdyby totiž λ2 měla nenulovou imaginární část, pak by také λ2 byla vlastní hodnotou různou od λ2 i λ1, což by bylo ve sporu s předpokladem, že λ1, λ2 jsou všechny vlastní hodnoty stejného modulu. To znamená, že ϕ2 = π. Buď d > 2. Je-li λ3 je vlastní hodnotou matice A, pak podle (iii) je e−iϕ2 λ3 vlastní hodnotou matice e−iϕ2 A, takže podle (i) je také vlastní hodnotou matice e−iϕ2 eiϕ2 D2AD−1 2 = D2AD−1 2 . Nyní podle (ii) je e−iϕ2 λ3 vlastní hodnotou matice A, což znamená, že existuje k ∈ {1, 2, . . . , d}, že λk = e−iϕ2 λ3, eiϕk λ1 = e−iϕ2 eiϕ3 λ1, eiϕk = ei(ϕ3−ϕ2) , ϕk = ϕ3 − ϕ2, neboť ϕk ∈ [0, 2π). To znamená, že ϕk ∈ (0, ϕ3). To je však možné jen tak, že ϕk = ϕ2, k = 2 a tedy λ3 = eiϕ2 λ2. Analogicky lze ukázat, že λ4 = eiϕ2 λ3, λ5 = eiϕ2 λ4, . . . , λd = eiϕ2 λd−1, λ1 = eiϕ2 λd. Odtud plyne, že vlastní hodnoty λ1, λ2, . . . , λd jsou vrcholy pravidelného d-úhelníku se středem 0 v komplexní rovině, tedy ϕ2 = 2π/d. Buď nyní j ∈ {2, 3, . . . , d} libovolný index. Z rovnosti Av = λ1v plyne rovnost A eiϕj v = λjv. Jsou-li tedy v(1), v(2), . . . , v(l) lineárně nezávislé vlastní vektory matice A příslušné k vlastní hodnotě λ1, pak eiϕj v(1), eiϕj v(2), . . . , eiϕj v(l) jsou vlastní vektory matice A příslušné k vlastní hodnotě λj, které jsou lineárně nezávislé. To znamená, že dim ker(A − λ1I) ≤ dim ker(A − λjI) . Analogicky z toho, že rovnost Aw = λjw implikuje rovnost A e−iϕj w = λ1w, odvodíme nerovnost dim ker(A − λjI) ≤ dim ker(A − λ1I) . Celkem tedy dostaneme, že platí rovnost dim ker(A − λ1I) = dim ker(A − λjI) . Tvrzení 14. Je-li A ≥ 0 ireducibilní, pak λ1 = sup SA > 0, příslušný vlastní vektor v > 0 a dim ker(A − λ1I) = 1. Důkaz. Nejprve ukážeme, že ireducibilní nezáporná matice A nemá nulový sloupec: Připusťme, že existuje index j takový, že aij = 0 pro všechna i = 1, 2, . . . , n. Pak pro libovolné m ∈ N je Am jj = Am−1 A jj = n k=1 (Am−1 )jkakj = 0, což je spor s ireducibilitou. 46 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Položme c = min    n j=1 aij : i = 1, 2, . . . , n    , v(c) = 1 √ n (1, 1, . . . , 1)T . Pak c > 0 a Av(c) i = 1 √ n n j=1 aij ≥ 1 √ n c = c v(c) i , tedy Av(c) ≥ cv(c), takže c ∈ SA. Odtud plyne λ1 = sup SA ≥ c > 0. Poněvadž matice A je ireducibilní, ke každé dvojici indexů i, j existuje číslo κij takové, že Aκij ij > 0. Položme m = max {κij : i = 1, . . . , n, j = 1, . . . , n, i = j}. Pak (I + A)m = I + m j=1 m j Aj > 0, a pro libovolný vlastní vektor v matice A příslušný k vlastní hodnotě λ1 platí (I + A)m v = =  I + m j=1 m j Aj   v = v + m j=1 m j λj 1v =  1 + m j=1 m j λj 1   v = (1 + λ1)m v, tedy v je současně nezáporný vlastní vektor kladné matice (I+A)m příslušný k vlastní hodnotě (1 + λ1)m. Podle tvrzení 2 je v > 0 a tedy dim ker(A − λ1I) ≥ 1. Z uvedeného výpočtu dále plyne, že ker (A − λ1I) ⊆ ker (I + A)m − (1 + λ1)m I . Prostor na pravé straně inkluze je podle tvrzení 5 jednodimenzionální a tedy dim ker (A − λ1I) ≤ 1. Věta 1. Buď A nezáporná matice. Pak existuje její vlastní hodnota λ1 ∈ R taková, že λ1 ≥ |λ| pro každou vlastní hodnotu λ matice A, a existuje nezáporný vlastní vektor v příslušný k vlastní hodnotě λ1. Je-li navíc matice A primitivní, pak λ1 > |λ| pro libovolnou vlastní hodnotu λ = λ1 matice A, příslušný vlastní vektor v > 0 a ker(A − λ1I) je jednodimenzionální. Je-li navíc matice A ireducibilní a imprimitivní, pak λ1 > 0, příslušný vlastní vektor v > 0 a existují vlastní hodnoty λ2, λ3, . . . , λd takové, že λj = e2πi(j−1)/dλ1 a ker(A − λ1I) je jednodimenzionální, j = 1, 2, . . . , d. Důkaz. První část je tvrzení 10, druhá část je tvrzení 11, třetí část je tvrzení 13 a 14. Poznámka 2. Číslo d z třetí části věty 1 je větší než 1. Tato vlastnost však nebyla dokázána. Klasifikace nezáporných matic a odpovídající vlastnosti vlastních hodnot a vlastních vektorů jsou shrnuty v obrázku 2.1. 2.3. ŘEŠENÍ PROJEKČNÍ ROVNICE 47 Reducibilní (∃i, j) (∀ ) A ij = 0 λ1 ≥ 0, w(1) ≥ 0 Ireducibilní (∀i, j) (∃ ) A ij > 0 λ1 > 0, w(1) > 0 (∃d) λj = λ1e2πi(j−1)/d, j = 1, 2, . . . , d 4 4 4 4 4 4 4 4 44 ˜ ˜ ˜ ˜ ˜ ˜ ˜ ˜ ˜˜ Imprimitivní (∀ ) (∃i, j) A ij = 0 Primitivní (∃ ) A > 0 λ1 > |λ2|, w(1) > 0 4 4 4 4 4 4 4 4 44 ˜ ˜ ˜ ˜ ˜ ˜ ˜ ˜ ˜˜ Nezáporná A ≥ 0 λ1 ∈ R Obrázek 2.1: Klasifikace nezáporných matic a charakterizace jejich vlastních hodnot a vlastních vektorů. Vlastí hodnoty λ1, . . . , λn matice A jsou uspořádané tak, že |λ1| ≥ | · · · ≥ |λn|, w(1) označuje vlastní vektor příslušný k vlastní hodnotě λ1. 2.3 Řešení projekční rovnice Budeme řešit projekční rovnici n(t + 1) = An(t) (2.24) s konstantní projekční maticí A typu k × k. Předpokládejme, že známe počáteční strukturu populace n(0). Pak platí n(1) = An(0), n(2) = An(1) = AAn(0) = A2 n(0), n(3) = An(2) = AA2 n(0) = A3 n(0), atd. Obecně n(t) = AAt−1n(0) a tedy n(t) = At n(0). (2.25) 48 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Přímým dosazením do rovnice (2.24) se lze přesvědčit, že n(t) dané rovností (2.25) je skutečně řešením projekční rovnice (2.24). Dále budeme předpokládat, že matice A v rovnici (2.24) má k různých vlastních hodnot. Z hlediska aplikací tento předpoklad není omezující. Aby totiž matice A měla násobnou vlastní hodnotu, musí její prvky splňovat určitou rovnost. Jinak řečeno, množina všech matic typu k × k majících násobné vlastní hodnoty tvoří varietu dimenze menší než k2 v prostoru Rk2 ; poněvadž k2-rozměrná míra takové variety je nulová, je (geometrická) pravděpodobnost jevu, že matice A bude mít násobné vlastní hodnoty, také nulová. Nechť různé vlastní hodnoty λ1, λ2, . . . , λk projekční matice A jsou uspořádány tak, že |λ1| ≥ |λ2| ≥ · · · ≥ |λk|; poněvadž matice A je nezáporná platí λ1 ∈ R, λ1 ≥ 0, a tedy |λ1| = λ1. Označme w(j) vlastní vektor matice A příslušný k vlastní hodnotě λj, j = 1, 2, . . . , k. Z předpokládané různosti vlastních hodnot plyne, že vektory w(1), w(2), . . . , w(k) tvoří bázi prostoru Rk. Existují tedy konstanty c1, c2, . . . , ck takové, že n(0) = c1w(1) + c2w(2) + · · · + ckw(k) . Dosazením tohoto vyjádření do rovnosti (2.25) dostaneme n(t) = At n(0) = k j=1 At cjw(j) = k j=1 cjAt−1 Aw(j) = k j=1 cjAt−1 λjw(j) = = k j=1 cjλjAt−2 Aw(j) = k j=1 cjλjAt−2 λjw(j) = k j=1 cjλ2 j At−2 w(j) = · · · = k j=1 cjλt jw(j) . Dostáváme tedy řešení rovnice (2.24) ve tvaru n(t) = k j=1 cjλt jw(j) . (2.26) Položme nyní W = w(1) , w(2) , . . . , w(k) T , Λ = diag(λ1, λ2, . . . , λk) = (δijλj) =      λ1 0 . . . 0 0 λ2 . . . 0 ... ... ... ... 0 0 . . . λk      ; W je matice, jejíž sloupce jsou vlastní vektory matice A, Λ je diagonální matice s vlastními hodnotami matice A na diagonále. Pak platí (AW)ij = Aw(j) i = λjw(j) i = λjw (j) i , (WΛ)ij = k l=1 w (l) i δljλj = λjw (j) i , 2.3. ŘEŠENÍ PROJEKČNÍ ROVNICE 49 a tedy AW = WΛ. (2.27) Poněvadž sloupce matice W tvoří bázi prostoru Rk, jsou lineárně nezávislé a tedy matice W je regulární. Z předchozí rovnosti proto dostaneme A = WΛW−1 , (2.28) dále W−1A = ΛW−1 a po transpozici AT WT−1 = WT−1 Λ; symbol WT−1 označuje matici W−1 T = WT −1 . Porovnáním s rovností (2.27) vidíme, že sloupce matice WT−1 jsou vlastní vektory transponované matice AT, která má stejné vlastní hodnoty jako původní matice A. Je-li tedy WT−1 = v(1), v(2), . . . , v(k) , pak AT v(j) = λjv(j) , neboli v(j)T A = λjv(j)T ; jinak řečeno, řádky matice W−1 jsou transponované levé vlastní vektory matice A. Platí tedy v(i)T w(j) = δij. Označme nyní c = W−1n(0). Vyjádření matice A ve tvaru (2.28) nyní dosadíme do řešení (2.25) rovnice (2.24): n(t) = At n(0) = WΛW−1 t n(0) = WΛW−1 WΛW−1 · · · WΛW−1 n(0) = WΛt W−1 n(0) = = WΛt c = w(1) , w(2) , . . . , w(k)      λt 1 0 . . . 0 0 λt 2 . . . 0 ... ... ... ... 0 0 . . . λt k           c1 c2 ... ck      = = w(1) , w(2) , . . . , w(k)      λt 1c1 λt 2c2 ... λt kck      = k j=1 cjλt jw(j) . Dostáváme tedy stejné vyjádření řešení rovnice (2.24) jako bylo v rovnosti (2.26). Navíc ale vidíme, že pro konstanty c1, c2, . . . , ck platí cj = v(j)T n(0), j = 1, 2, . . . , k. Provedené úvahy lze shrnout do věty: Věta 2. Nechť nezáporná matice A má různé vlastní hodnoty λ1, λ2, . . . , λk takové, že λ1 ≥ |λ2| ≥ · · · ≥ |λk|. Označme w(j), resp. v(j), (pravý) vlastní vektor, resp levý vlastní vektor, matice A příslušný k vlastní hodnotě λj, j = 1, 2, . . . , k. Nechť vektory v(i), w(j) jsou takové, že v(i)T w(j) = δij. Pak řešení projekční rovnice (2.24) je tvaru n(t) = k j=1 cjλt jw(j) , (2.29) kde cj = v(j)T n(0). 50 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Poznamenejme, že z předpokladu různosti vlastních hodnot plyne, že λ1 > 0. V opačném případě by totiž bylo |λ2| = |λ3| = · · · = |λk| = 0. Pro ireducibilní matici A podle PerronovyFrobeniovy věty (aplikované na matici AT) platí pro ireducibilní matici A nerovnost v(1) > 0. Má-li tedy počáteční struktura populace n(0) v takovém případě alespoň jednu složku nenulovou (tj. je-li populace na začátku sledování přítomna), pak je c1 = v(1)T n(0) > 0. Řekneme, že rovnice (2.24) je ergodická, pokud průběh jejího řešení v okolí nekonečna (tj. pro dostatečně velké t) nezávisí na počáteční podmínce n(0). Populaci, jejíž vývoj je popsán ergodickou rovnicí, nazveme také ergodickou. 2.3.1 Matice A primitivní Řešení (2.29) projekční rovnice (2.24) přepíšeme na tvar n(t) = λt 1  c1w(1) + k j=2 cj λj λ1 t w(j)   . V tomto případě je podle Perronovy-Frobeniovy věty λ1 > |λj| pro j = 2, 3, . . . , k a tedy lim t→∞ 1 λt 1 n(t) − c1w(1) = k j=2 cj lim t→∞ λj λ1 t w(j) = o. Řešení n(t) projekční rovnice (2.24) s primitivní maticí A je tedy pro libovolnou počáteční hodnotu n(0) asymptoticky ekvivalentní s funkcí c1λt 1w(1). To znamená, že pro dostatečně velké t nezávisle na počáteční struktuře n(0) populace (pokud je alespoň jedna její složka na začátku nenulová) roste velikost populace exponenciálně a relativní zastoupení jejich jednotlivých složek je úměrné složkám kladného vlastního vektoru w(1) příslušného k dominantní vlastní hodnotě λ1. Populace s primitivní projekční maticí A je tedy ergodická. Dominantní vlastní hodnotu λ1 matice A lze interpretovat jako Malthusovský koeficient růstu. Pokud tedy λ1 > 1, populace roste, pokud λ1 < 1, populace vymírá. 2.3.2 Matice A ireducibilní a imprimitivní Podle Perronovy-Frobeniovy věty je v tomto případě λ1 > 0 a existuje přirozené číslo d takové, že 1 < d ≤ k, λj = λ1e2πi(j−1)/d pro j = 1, 2, . . . , d a λ1 > |λd+1| pokud d < k. Přitom dim ker(A − λjI) = 1 pro j = 1, 2, . . . , d. Řešení (2.29) rovnice (2.24) přepíšeme na tvar n(t) = λt 1   d j=1 cje2πi(j−1)t/d w(j) + k j=d+1 cj λj λ1 t w(j)   ; při zápisu používáme konvenci p−1 j=p γj = 0. Platí tedy lim t→∞   1 λt 1 n(t) − d j=1 cje2πi(j−1)t/d w(j)   = k j=d+1 cj lim t→∞ λj λ1 t w(j) = o. (2.30) 2.3. ŘEŠENÍ PROJEKČNÍ ROVNICE 51 Vidíme, že řešení n(t) projekční rovnice (2.24) s ireducibilní a imprimitivní maticí A je tedy pro libovolnou počáteční hodnotu n(0) asymptoticky ekvivalentní s funkcí λt 1s(t), kde s = d j=1 cje2πi(j−1)t/d w(j) je d-periodická funkce. Nechť j je index takový, že 1 < j ≤ 1 2(d + 1). Pak d − j + 2 ≤ d a pro vlastní hodnotu λd−j+2 platí λd−j+2 = λ1e2πi(d−j+1)/d = λ1e2πi e−2πi(j−1)/d = λ1e−2πi(j−1)/d = λj, tj. vlastní hodnoty λj a λd−j+2 jsou komplexně sdružené. Matice A je reálná a proto pro vlastní vektor w(j) příslušný k vlastní hodnotě λj platí Aw(j) = Aw(j) = Aw(j) = λjw(j) = λjw(j) = λd−j+2w(j). Vektor w(j) je tedy vlastním vektorem matice A příslušným k vlastní hodnotě λd−j+2. Poněvadž dim ker(A − λd−j+2I) = 1, existuje konstanta αj taková, že w(d−j+2) = αjw(j). Podobně lze ukázat, že existuje konstanta βj taková, že pro levé vlastní vektory v(j) a v(d−j+2) platí v(d−j+2) = βjv(j). Dále platí 1 = v(d−j+2)T w(d−j+2) = βjv(j) T αjw(j) = αjβjv(j)T w(j) = αjβj. Poněvadž počáteční struktura populace n(0) je reálný vektor, platí cd−j+2 = v(d−j+2)T n(0) = βjv(j) T n(0) = βjv(j)T n(0) = βjcj a dále cje2πi(j−1)t/d w(j) + cd−j+2e2πi(d−j+1)t/d w(d−j+2) = = cje2πi(j−1)t/d w(j) + βjcje−2πi(j−1)t/d αjw(j) = = cje2πi(j−1)t/d w(j) + cje2πi(j−1)t/dw(j) = 2 cje2πi(j−1)t/d w(j) ; symbol označuje reálnou část komplexního čísla (vektoru). Je-li číslo d sudé, pak pro j = d 2 + 1 platí cje2πi(j−1)t/d = cd 2 +1eπit = cd 2 +1(−1)t . Položme nyní r(t) = (−1)tcd 2 +1w(d 2 +1), d sudé, o, d liché. Poznamenejme, že řeálné vlastní hodnotě λd 2 +1 matice A odpovídá reálný levý i pravý vlastní vektor. Proto je funkce r reálná. Dále platí d−1 l=0 r(t + l) = 0. 52 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Vektorovou funkci s(t) = d j=1 cje2πi(j−1)t/d w(j) nyní přepíšeme na tvar s(t) = c1w(1) + [d+1 2 ] j=2 cje2πi(j−1)t/d w(j) + cd−j+2e2πi(d−j+1)t/d w(d−j+2) + r(t) = = c1w(1) + 2 [d+1 2 ] j=2 cje2πi(j−1)t/d w(j) + r(t); symbol [γ] označuje celou část z reálného čísla γ. Z tohoto vyjádření vidíme, že funkce s(t) je reálná. Pro průměr hodnot d-periodické funkce s na intervalu délky periody platí 1 d d−1 l=0 s(t + l) = c1w(1) + 1 d d−1 l=0 2 [d+1 2 ] j=2 cje2πi(j−1)(t+l)/d w(j) + o = = c1w(1) + 2 d    [d+1 2 ] j=2 cje2πi(j−1)t/d w(j) d−1 l=0 e2πi(j−1)l/d    = = c1w(1) + 2 d    [d+1 2 ] j=2 cje2πi(j−1)t/d w(j) 1 − e2πi(j−1)d/d 1 − e2πi(j−1)/d    = c1w(1) . To vzhledem k (2.30) znamená, že pro dostatečně velký čas t nezávisle na počáteční struktuře n(0) populace (pokud je ovšem alespoň jedna z jejích složek nenulová) roste populace tak, že její velikost kolísá kolem exponenciální funkce a dlouhodobý průměr zastoupení jednotlivých složek je úměrný složkám vlastního vektoru w(1) příslušného k dominantní vlastní hodnotě λ1. Populace je opět ergodická a dominantní vlastní hodnotu λ1 matice A lze opět interpretovat jako Malthusovský koeficient růstu. 2.3.3 Matice A reducibilní Věta 3. Matice A je reducibilní právě tehdy, když její řádky a sloupce lze přeskládat tak, že ji lze blokově zapsat ve tvaru A = B1 O B12 B2 kde B1 je ireducibilní matice typu k1 × k1, B2 je matice typu (k − k1) × (k − k1) a B12 je matice typu (k − k1) × k1; přitom 1 ≤ k1 < k. Bez újmy na obecnosti lze tedy matici A přepsat v blokovém tvaru A = B1 O B12 B2 2.3. ŘEŠENÍ PROJEKČNÍ ROVNICE 53 kde B1 je ireducibilní matice typu k1 × k1. Pak A2 = B1 O B12 B2 B1 O B12 B2 = B2 1 O B12B1 + B2B12 B2 2 , A3 = B2 1 O B12B1 + B2B12 B2 2 B1 O B12 B2 = B3 1 O (B12B1 + B2B12)B1 + B2 2B12 B3 2 , atd. Obecně At =   Bt 1 O B (t) 12 Bt 2   , kde B (t) 12 je nezáporná matice typu (k − k1) × k1. Vektor n popisující strukturu populace vyjádříme jako n = q p , kde vektor q je k1-rozměrný a vektor p je (k − k1)-rozměrný. Řešení (2.25) projekční rovnice (2.24) je nyní tvaru n(t) = Bt 1 O B (t) 12 Bt 2 q(0) p(0) = Bt 1q(0) B (t) 12 q(0) + Bt 2p(0) . Modelovanou populaci tedy můžeme rozdělit na „ireducibilní část q a „zbytek p (ten je tvořen např. postreprodukčními stadii nebo subpopulací na stanovišti, na něž vedou migrační cesty z „hlavního areálu ale nikoliv zpět v případě prostorově strukturovaných modelů ap.). „Ireducibilní část populace se vyvíjí způsobem popsaným v 2.3.1 nebo v 2.3.2 podle toho, zda je matice B primitivní nebo imprimitivní. 2.3.4 Stabilizovaná struktura a reproduktivní hodnota Buď A ireducibilní matice, λ1 její dominantní vlastní hodnata a w = (w1, w2, . . . , wk)T příslušný vlastní vektor takový, že k j=1 wj = 1. Poznamenejme, že podle Perronovy-Frobeniovy věty je w > 0. Buď dále v = (v1, v2, . . . , vk)T levý vlastní vektor matice A příslušný k vlastní hodnotě λ1 takový, že vTw = 1. Opět je v > 0. Nechť n = n(t) je řešení projekční rovnice (2.24) a c1 = vTn(0). Podle 2.3.1 a 2.3.2 existuje τ ∈ N, že platí lim t→∞ τ−1 i=0 1 λi 1 n(t + i) − c1λt 1w = o; v případě primitivní matice A stačí volit τ = 1, v případě imprimitivní matice A stačí volit τ = d. Po dostatečně dlouhém vývoji populace tedy jsou průměrné velikosti jejích jednotlivých složek úměrné složkám vektoru w, tento vektor vyjadřuje stabilizovanou strukturu populace. Uvažujme nyní strukturně stabilizovanou populaci, tj. populaci po dostatečně dlouhém vývoji, jejíž průměrná struktura se dále vyvíjí podle (přibližné) rovnosti ¯n(t) = 1 τ τ−1 i=0 1 λ1 n(t + i) = c1λt 1w = λt 1vT n(0)w 54 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ a jejíž celková průměrná velikost je tedy (přibližně) rovna k j=1 ¯nj(t) = k j=1 λt 1vT n(0)wj = λt 1vT n(0) k j=1 wj = λt 1vT n(0) = λt 1 k j=1 vjnj(0). Při označení ξj = vj k l=1 vl je celková průměrná velikost populace (přibližně) rovna k j=1 ¯nj(t) = λt 1 k l=1 vl k j=1 ξjnj(0). Poněvadž k j=1 ξj = 1, představuje výraz k j=1 ξjnj(0) váženým průměr velikostí složek počáteční populace. Výsledek lze tedy interpretovat tak, že ξj vyjadřuje váhu, s jakou přispívá j-tá složka počáteční populace k velikosti populace po dostatečně dlouhém vývoji. Hodnota ξj se proto nazývá reproduktivní hodnota j-té složky populace, vektor 1 k l=1 vl      v1 v2 ... vk      se nazývá vektor reproduktivních hodnot složek populace. Pokud ve struktuře populace je jediná třída novorozenců, může být užitečné vyjadřovat reproduktivní hodnotu jednotlivých tříd populace relativně vzhledem k reproduktivní hodnotě novorozenců. Je-li tedy v takovém případě první třída třídou novorozenců, uvažujeme vektor reproduktivních hodnot ve tvaru 1 vl      v1 v2 ... vk      ; při čtení literatury je tedy potřebné dávat pozor, kterou z „variant definice vektoru reproduktivních hodnot autor používá. 2.4 Transientní dynamika V celém oddílu bude A ireducibilní matice typu k × k, (k ≥ 2), λ1 > 0 její dominantní vlastní hodnota (růstový koeficient), w = (w1, w2, . . . , wk)T příslušný (pravý) vlastní vektor takový, že k j=1 wj = 1 (stabilizovaná struktura populace) a v příslušný levý vlastní vektor takový, že vTw = 1. Vektor n = n(t) = n1(t), n2(t), . . . , nk(t) T bude řešením rovnice (2.24) v čase t. Budeme 2.4. TRANSIENTNÍ DYNAMIKA 55 dále používat označení c1 = vTn(0) a ¯n(t) = 1 d d−1 l=0 1 λi 1 n(t + l), kde d je počet vlastních hodnot matice A které mají modul rovný hodnotě λ1. Symbolem · 1 budeme označovat „taxíkářskou normu na Rk, tj. pro libovolný vektor ξ = (ξ1, ξ2, . . . , ξk)T ∈ Rk klademe ξ 1 = k j=1 |ξj|. Norma n(t) 1 vyjadřuje celkovou velikost populace v čase t; absolutní hodnoty v součtu není třeba psát, neboť vektor n(0) je nezáporný. 2.4.1 Rychlost konvergence ke stabilizované struktuře Buď µ největší z modulů vlastních hodnot matice A menších než λ1; v případě d = k klademe µ = 0. Koeficient tlumení (dumping ratio) definujeme jako = λ1 µ , pro µ = 0 klademe = ∞. Zřejmě je > 1. Porovnáním s výsledky oddílu 2.3 vidíme, že existuje konstanta K taková, že 1 λt 1 ¯n(t) − c1w ≤ K −t (2.31) pro libovolnou normu · na Rk ekvivalentní s euklidovskou. Koeficient tlumení tedy vyjadřuje rychlost konvergence ke stabilizované struktuře. V případě primitivní matice A lze nerovnost (2.31) přepsat na tvar K −t > 1 λt 1 n(t) − wc1 = 1 λ1 A t n(0) − wvT n(0) = 1 λ1 A t − wvT n(0) . Ke každému nezápornému vektoru x tedy existuje konstanta K = K(x) taková, že 1 λ1 A t − wvT x < K −t . Odtud dále plyne, že existuje kladná konstanta C taková, že 1 λ1 A t − wvT ij < C −t (2.32) pro všechna i, j = 1, 2, . . . , k. Tato vlastnost umožňuje zformulovat a dokázat jeden výsledek z teorie primitivních matic: 56 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Věta 4. Buď A primitivní matice typu k × k, λ1 > 0 její dominantní vlastní hodnota, w, resp. v, její pravý, resp. levý, vlastní vektor příslušný k dominantní vlastní hodnotě, tj. Aw = λ1w, vT A = λ1vT a nechť platí vTw = 1. Pak matice I + wvT − 1 λ1 A je regulární a řada ∞ i=1 1 λ1 A i − wvT absolutně konverguje. Přitom platí I + wvT − 1 λ1 A −1 = I + ∞ i=1 1 λ1 A i − wvT . (2.33) Důkaz. Buďte i, j ∈ N, i ≥ j libovolné indexy. Pak platí 1 λ1 A i−j wvT j = 1 λ1 A i−j−1 1 λ1 AwvT wvT j−1 = = 1 λ1 A i−j−1 wvT j = · · · = 1 λ1 A 0 wvT j = wvT j = = wvT wvT · · · wvT = w vT w vT w · · · vT w vT = = wvT , a tedy s využitím binomické věty dostaneme 1 λ1 A − wvT i = i j=0 i j 1 λ1 A i−j (−1)j wvT j = = 1 λ1 A i + i j=1 i j 1 λ1 A i−j (−1)j wvT j = = 1 λ1 A i + i j=1 i j (−1)j wvT = = 1 λ1 A i +   i j=0 i j (−1)j − 1   wvT = = 1 λ1 A i + (1 − 1)i − 1 wvT = = 1 λ1 A i − wvT . Řada ∞ i=1 1 λ1 A i − wvT konverguje absolutně, neboť podle nerovnosti (2.32) je majorizována konvergentní řadou. Z předchozího výpočtu plyne, že absolutně konverguje ke stejnému 2.4. TRANSIENTNÍ DYNAMIKA 57 součtu také geometrická řada ∞ i=1 1 λ1 A − wvT i . Platí tedy I + ∞ i=1 1 λ1 A i − wvT = I + ∞ i=1 1 λ1 A − wvT i = ∞ i=0 1 λ1 A − wvT i = = I − 1 λ1 A − wvT −1 , což je dokazovaná rovnost (2.33). 2.4.2 Vzdálenost od stabilizované struktury Označme x(t) = n(t) n(t) 1 = 1 k j=1 nj(t) n(t). Keyfitzova vzdálenost populace od stabilizované struktury je definována jako ∆ n(t), w = 1 2 x(t) − w 1 = 1 2 k j=1 |xj(t) − wj| ; jedná se o standardní vzdálenost pravděpodobnostních vektorů. Poněvadž všechna xj(t) jsou nezáporná a všechna wj jsou kladná, je |xj(t) − wj| ≤ |xj(t)| + |wj| = xj(t) + wj a tedy 0 ≤ ∆ n(t), w ≤ 1 2 k j=1 xj(t) + wj = 1. Rovnost ∆ n(t), w = 0 přitom nastane právě tehdy, když xj(t) = wj pro všechny indexy j = 1, 2, . . . , k, tedy právě tehdy, když x(t) = w. Keyfitzova vzdálenost od vektoru w závisí pouze na hodnotě n(t), nikoliv na trajektorii, po níž se struktura populace n(t) ke stabilizované struktuře w přibližuje. Pro jemnější hodnocení odchylky populace od stabilizované struktury je potřebné tuto trajektorii zohlednit. Položme proto s A, n(0), t = t i=0 1 λi 1 ¯n(i) − c1w . Tento vektor kumuluje rozdíly aktuální „průměrné struktury populace od struktury stabilizované. Cohenova kumulativní vzdálenost počáteční struktury populace n(0) od struktury stabilizované je definována jako D A, n(0) = k j=1 lim t→∞ sj A, n(0), t = lim t→∞ s A, n(0), t 1 . 58 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Nechť matice A je primitivní, tedy d = 1, ¯n(t) = n(t) = Atn(0) podle (2.25). V tomto případě tedy s využitím Věty 4 dostaneme lim t→∞ s A, n(0), t = ∞ i=0 1 λi 1 n(i) − wc1 = ∞ i=0 1 λ1 A i n(0) − wvT n(0) = = ∞ i=0 1 λ1 A i − wvT n(0) = I − wvT + ∞ i=1 1 λ1 A i − wvT n(0) = = I + wvT − 1 λ1 A −1 − wvT n(0). Označme Z = I + wvT − 1 λ1 A −1 . Pak Cohenovu vzdálenost n(0) od stabilizované struktury můžeme vyjádřit vztahem D A, n(0) = Z − wvT n(0) 1 . 2.4.3 Populační moment Předpokládejme, že populace v čase t = 0 má stabilizovanou strukturu, ale nikoliv stabilizovanou velikost, tj. populace roste nebo vymírá. V čase t = 0 se nějakým vnějším zásahem změní projekční matice tak, aby její dominantní vlastní hodnota byla rovna 1, tj. aby se stabilizovala i velikost populace. Označme Aold projekční matici populace v čase t < 0 a Anew projekční matici populace v čase t ≥ 0. Populační moment (population momentum) definujeme jako M = lim t→∞ n(t) 1 n(0) 1 , (2.34) pokud tato limita existuje. Poněvadž struktura populace n(t) v čase t ≥ 0 závisí pouze na matici Anew a na počáteční struktuře populace n(0), populační moment závisí na n(0) a Anew, M = M n(0), Anew . V případě M > 1 velikost populace po stabilizaci vzroste, v případě M < 1 se zmenší. Dominantní vlastní hodnota matice Anew je rovna 1. Pravý, resp. levý, vlastní vektor matice Anew příslušný k vlastní hodnotě 1 označíme wnew, resp. vnew. Nechť wnew 1 = 1, vT newwnew = 1 a matice Anew je primitivní. V tomto případě podle 2.3.1 je lim t→∞ n(t) = vT newn(0)wnew. Odtud zejména plyne, že limita v definici populačního momentu (2.34) existuje. Platí tedy M = vT newn(0)wnew 1 n(0) 1 = vT newn(0) wnew 1 n(0) 1 = vT newn(0) n(0) 1 . Nechť matice Aold je také primitivní a wold je vlastní vektor příslušný k dominantní vlastní hodnotě matice Aold takový, že wold 1 = 1. V takovém případě je n(0) = c1wold a M = vT newc1wold c1wold 1 = vT newwold. Jsou-li tedy obě matice Aold, Anew primitivní, je populační moment těmito maticemi jednoznačně určen, M = M(Aold, Anew). 2.5. ANALÝZA CITLIVOSTI A PRUŽNOSTI 59 2.5 Analýza citlivosti a pružnosti Vývoj populace podle rovnice (2.24) je charakterizován demografickými charakteristikami — např. růstovým koeficientem λ1, stabilizovanou věkovou strukturou w, reproduktivní hodnotou složek populace a podobně. Tyto charakteristiky závisí na projekční matici A. Určitou charakteristiku však jednotlivé složky matice A neovlivňují stejnou měrou: charakteristika je více či méně citlivá na změny nějaké složky projekční matice; charakteristika reaguje na změny složky matice A více či méně pružně. Citlivost charakteristiky χ = χ(A) na složku aij projekční matice A (sensitivity of χ to aij) je definována jako ∂χ ∂aij , pružnost charakteristiky χ = χ(A) vzhledem ke složce aij projekční matice A (elasticity of χ with respect to aij) je definována jako aij χ ∂χ ∂aij = ∂ ln χ ∂ ln aij . 2.5.1 Citlivost a pružnost růstového koeficientu Nechť λ je vlastní hodnota matice A, w, resp. v je pravý, resp. levý, vlastní vektor příslušný k vlastní hodnotě λ. Rovnost Aw = λw zderivujeme podle aij, vynásobíme zleva vektorem vT a upravíme pomocí vztahu vTA = λvT. Tímto způsobem dostaneme Aw = λw, ∂A ∂aij w + A ∂w ∂aij = ∂λ ∂aij w + λ ∂w ∂aij vT ∂A ∂aij + vT A ∂w ∂aij = vT ∂λ ∂aij w + λvT ∂w ∂aij viwj = ∂λ ∂aij vT w, tedy ∂λ ∂aij = viwj vTw . (2.35) V případě, že matice A je ireducibilní, λ1 je dominantní vlastní hodnota projekční matice A (Malthusovský koeficient růstu), w a v jsou pravý a levý vlastní vektor příslušný k dominantní vlastní hodnotě, které splňují rovnost vTw = 1, dostaneme ∂λ1 ∂aij = viwj. Nyní můžeme položit sij = viwj a definovat matici citlivosti S růstového koeficientu λ1 na složky projekční matice A jako S = (sij) = ∂λ1 ∂aij = (viwj) = vwT . Matice citlivosti vyjadřuje vliv změn populačních prametrů na růstový koeficient. A to včetně změn těch parametrů. které se v reálné populaci měnit nemohou, neboť jsou nutně nulové 60 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ (např. nelze přeskočit některé vývojové stadium hmyzu). Citlivost tedy vyjadřuje, co by se stalo, kdyby se jistý parametr změnil nebo mohl změnit. I tento hypotetický výsledek může být v některých situacích zajímavý (např. jaký vliv na evoluční zdatnost populace by měla mutace způsobující přechod ze stadia larvy přímo v dospělce bez stadia kukly). Pružnost eij růstového koeficientu λ1 vzhledem ke složce aij je nyní dána rovností eij = aij λ1 ∂λ1 ∂aij = 1 λ1 aijsij = 1 λ1 aijviwj, matice pružnosti růstového koeficientu je definována jako E = (eij) = 1 λ1 A ◦ vwT , kde ◦ označuje Hadamardův součin matic (součin „po složkách ). Lemma 1 (Eulerova věta o homogenních funkcích). Je-li funkce f = f(x1, x2, . . . , xm) homogenní řádu κ, tj. f(cx1, cx2, . . . , cxm) = cκ f(x1, x2, . . . , xm) (2.36) pro libovolnou konstantu c, pak m i=1 xi ∂f(x1, x2, . . . , xm) ∂xi = κf(x1, x2, . . . , xm). Důkaz. Rovnost (2.36) zderivujeme podle c, tj. m i=1 ∂f(cx1, cx2, . . . , cxm) ∂cxi ∂cxi ∂c = κcκ−1 f(x1, x2, . . . , xm). Poněvadž ∂cxi ∂c = xi, stačí položit c = 1 abychom dostali dokazovanou rovnost. Pro růstový koeficient λ1 = λ1(a11, . . . , aij, . . . , akk) platí Aw = λ1w a také cAw = cλ1w pro libovolnou konstantu c. To znamená, že c-násobek vlastní hodnoty λ1 je vlastní hodnotou matice cA, cλ1(a11, . . . , aij, . . . , akk) = λ1(ca11, . . . , caij, . . . , cakk), jinak řečeno, růstový koeficient λ1 je homogenní funkcí řádu 1 složek projekční matice A. Podle Eulerovy věty o homogenních funkcích tedy platí k i,j=1 eij = 1 λ1 k i,j=1 aij ∂λ1 ∂aij = 1 λ1 λ1 = 1. Z tohoto důvodu bývá pružnost růstového koeficientu eij vzhledem ke složce aij interpretována jako relativní příspěvek složky aij projekční matice k růstovému koeficientu. 2.6. ANALÝZA VĚKOVĚ STRUKTUROVANÉ POPULACE 61 2.6 Analýza věkově strukturované populace Projekční maticí věkově strukturované populace je Leslieho matice A =          F1 F2 F3 . . . Fk−1 Fk P1 0 0 . . . 0 0 0 P2 0 . . . 0 0 ... ... ... ... ... ... 0 0 0 . . . 0 0 0 0 0 . . . Pk−1 0          . Parametry P1, P2, . . . , Pk−1, F1, F2, . . . , Fk jsou nezáporné; Fj označuje specifickou fertilitu jedinců j-té věkové třídy, Pj pravděpodobnost přežití projekčního intervalu jedincem j-té věkové třídy. Budeme předpokládat, že 0 < Pj < 1 pro j = 1, 2, . . . , k−1 (je možné dosáhnout maximálního věku k − 1, tj. být v k-té věkové třídě, a v libovolném věku je možné uhynout). Pro zjednodušení zápisu zavedeme symbol Pk a položíme Pk = 0. Nejprve (v pododdílech 2.6.1 a 2.6.2) provedeme pravděpodobnostní úvahy, které nevychází z teorie nezáporných matic. Pomocí nich odvodíme některé standardní demografické charakteristiky populace. Ty bývají používány nejen v demografii, ale i v pojistné matematice. Některé z nich umožní zpřehlednit vyjadřévání. Pokud budeme předpokládat, že Fk > 0 (i nejstarší jedinci jsou plodní), pak je matice A ireducibilní. Kdyby se v populaci vyskytovali i jedinci v postreprodukčním věku, byla by taková matice A projekční maticí reproduktivní části populace. 2.6.1 Čistá míra reprodukce Označme mj(t0) počet novorozených potomků rodičů z j-té věkové třídy v jistém čase t0, tj. mj(t0) = Fjnj(t0 − 1). Pak mj(t0 + i) je počet jedinců věku i, kteří se narodili v čase t0 a jejichž rodiče měli v čase t0 věk j. Úmrtí v různých časových okamžicích považujeme za stochasticky nezávislé jevy. Za tohoto předpokladu je klasická pravděpodobnost, že se jedinec dožije věku j svých rodičů v době svého narození, rovna mj(t0 + j) mj(t0) = j−1 q=1 Pq. Odtud mj(t0 + j) = mj(t0) j−1 q=1 Pq =  Fj j−1 q=1 Pq   nj(t0 − 1). Hodnota Fj j−1 q=1 Pq (2.37) tedy vyjadřuje počet potomků jednoho jedince věku j, kteří se dožili alespoň věku svého rodiče v době narození. Tato hodnota se nazývá fertilitní funkce. Součet R0 = k j=1 Fj j−1 q=1 Pq. (2.38) 62 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ se nazývá čistá míra reprodukce (net reproduction rate). Vyjadřuje poměr, v jakém nahradí generace potomků generaci svých rodičů; je-li tedy tento poměr menší než 1, nestačí následující generace nahradit generaci předcházející a populace vymírá. Fertilitní funkci (2.37) lze interpretovat i jiným způsobem. Výraz j−1 q=1 Pq vyjadřuje pravděpodobnost, že se jedinec dožije věku j −1, vstoupí do j-té věkové třídy. Jinak řečeno, v čase t představuje podíl živých jedinců mezi všemi, kteří se narodili v čase t − j. Fertlitní funkce násobená počtem jedinců narozených v čase t−j vyjadřuje množství potomků, kteří se těmto jedincům v čase t narodili. To znamená, že fertilitní funkce (2.37) představuje očekávaný (střední) počet potomků jedince v čase j po jeho narození. Při této interprataci čistá míra reprodukce (2.38) je očekávaný počet potomků jedince během jeho celého života. 2.6.2 Očekávaná doba dožití Nechť náhodná veličina T vyjadřuje celkovou dobu života nějakého jedince z uvažované populace. Uvažujme kohortu (skupinu jedinců, kteří se narodili ve stejném čase t0) o počáteční velikosti N0. Velikost kohorty v čase t označíme N(t); je tedy N(t0) = N0. Úmrtí v různých časových okamžicích považujeme za stochasticky nezávislé jevy. Klasická pravděpodobnost jevu, že jedinec z kohorty bude žít ještě ve věku j je proto dána výrazem N(t0 + j) N0 = P1P2 · · · Pj = j q=1 Pq. Odtud N(t0 + j) = N0 j q=1 Pq. Klasická pravděpodobnost, že doba života jedince je právě j, tj. že se jedinec dožije věku j a věku j + 1 se nedožije, je rovna P(T = j) = N(t0 + j) − N(t0 + j + 1) N0 = N0 j q=1 Pq − N0 j+1 q=1 Pq N0 = (1 − Pj+1) j q=1 Pq pro j = 1, 2, . . . , k − 1. Dále P(t = k) = 0, poněvadž k je věk, kterého již není možné dosáhnout. Střední délka života (očekávaná doba dožití, life expectancy) je definována jako střední hodnota náhodné veličiny T a je tedy dána formulí E T = k−1 j=1 j(1 − Pj+1) j q=1 Pq = k j=1 (j − 1)(1 − Pj) j−1 q=1 Pq = = k j=1 (j − 1) j−1 q=1 Pq − k−1 j=1 (j − 1) j q=1 Pq = k−1 j=1 j j q=1 Pq − k−1 j=1 (j − 1) j q=1 Pq = = k−1 j=1 j q=1 Pq = k j=2 j−1 q=1 Pq. (2.39) 2.6. ANALÝZA VĚKOVĚ STRUKTUROVANÉ POPULACE 63 Určíme ještě rozptyl náhodné veličiny T. K tomu nejprve vypočítáme E T2 = k−1 j=1 j2 P(T = j) = k−1 j=1 j2 (1 − Pj+1) j q=1 Pq = k j=1 (j − 1)2 (1 − Pj) j−1 q=1 Pq = = k j=1 (j − 1)2 j−1 q=1 Pq − k j=1 (j − 1)2 j q=1 Pq = k−1 j=1 j2 j q=1 Pq − k j=1 (j − 1)2 j q=1 Pq = = k−1 j=1 j2 − (j − 1)2 j q=1 Pq = k−1 j=1 (2j − 1) j q=1 Pq = 2 k−1 j=1 j j q=1 Pq − k−1 j=1 j q=1 Pq. Rozptyl doby života tedy je var T = E T2 − (E T)2 = 2 k−1 j=1 j j q=1 Pq − k−1 j=1 j q=1 Pq −   k−1 j=1 j q=1 Pq   2 . Označme dále Ta délku života, který má před sebou jedinec ve věku a. Pak je P(Ta = j) = N(t0 + a + j) − N(t0 + a + j + 1) N(t0 + a) = a+j q=1 Pq − a+j+1 q=1 Pq a q=1 Pq = (1 − Pa+j+1) a+j q=a+1 Pq pro j = 0, 1, . . . , k − a − 1 a střední délka života ve věku a (očekávaná doba dožití ve věku a) označovaná symbolem ea je dána výrazem ea = E Ta = k−a−1 j=0 j(1 − Pa+j+1) a+j q=a+1 Pq = k−a j=2 a+j−1 q=a+1 Pq; použili jsme analogické úpravy jako při odvození (2.39). Poznamenejme, že T0 = T a tedy očekávaná doba dožití je střední délkou života při narození. V demografických studiích se kromě střední délky života také udávají dvě další charakteristiky přežití. Pravděpodobná délka života ve věku a označovaná a je doba, po jejímž uplynytí zůstane na živu polovina jedinců z původního rozsahu. Přesněji řečeno, a splňuje nerovnosti N(t0 + a + a) ≤ 1 2N(t0 + a) < N(t0 + a + a − 1), neboli a+ a q=1 Pq ≤ 1 2 a q=1 Pq < a+ a−1 q=1 Pq, po úpravě a+ a q=a+1 Pq ≤ 1 2 < a+ a−1 q=a+1 Pq, takže pravděpodobnou délku života ve věku a můžeme vyjádřit formulí a = min    j : j = 0, 1, . . . , k − a, a+j q=a+1 Pq ≤ 1 2    . Normální délka života ve věku a označovaná a je doba, po jejímž uplynutí je úmrtí nejpravděpodobnější, tj. a = arg max {P(Ta = j) : j = 0, 1, . . . , k − a − 1}. 64 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ 2.6.3 Růstový koeficient populace Najdeme vlastní hodnoty matice A. Označme ∆k = det(A − λI) = F1 − λ F2 F3 . . . Fk−1 Fk P1 −λ 0 . . . 0 0 0 P2 −λ . . . 0 0 ... ... ... ... ... ... 0 0 0 . . . −λ 0 0 0 0 . . . Pk−1 −λ . Determinant ∆k rozvineme podle posledního sloupce, ∆k = (−1)k+1 Fk k−1 q=1 Pq − λ∆k−1. (2.40) Odtud je vidět, že pro λ = 0 je ∆k = (−1)k+1 Fk k−1 q=1 Pq = 0, pokud předpokádáme Fk = 0, a tedy matice A v souladu s Perronovou-Frobeniovou větou nemá nulové vlastní hodnoty. Rovnost (2.40) lze považovat za lineární diferenční rovnici (rekurentní formuli) prvního řádu pro neznámou posloupnost {∆k}. Jejím řešením je ∆k = (−λ)k−1 ∆1 − (−1)k k p=2 λk−p Fp p−1 q=1 Pq. Poněvadž ∆1 = F1 − λ, dostaneme ∆k = (−λ)k + (−λ)k−1 F1 − (−1)k k p=2 λk−p Fp p−1 q=1 Pq = (−λ)k  1 − k p=1 λ−p Fp p−1 q=1 Pq   . Vlastní hodnoty matice A tedy jsou řešením rovnice k p=1  Fp p−1 q=1 Pq   λ−p = 1. (2.41) Levou stranu této rovnice můžeme považovat za funkci proměnné λ, f(λ) = k p=1  Fp p−1 q=1 Pq   λ−p . Při tomto označení je f(1) = R0 podle (2.38). Dále lim λ→0+ f(λ) = ∞, lim λ→∞ f(λ) = 0 a f (λ) = − k p=1  pFp p−1 q=1 Pq   λ−p−1 < 0 2.6. ANALÝZA VĚKOVĚ STRUKTUROVANÉ POPULACE 65 f λ 1 R0 λ1 1 f λ R0 1 1 λ1 Obrázek 2.2: Grafické řešení rovnice (2.41) — charakteristické rovnice Leslieho matice. Vlevo: čistá míra reprodukce R0 < 1 (vymírající populace), vpravo: R0 > 1 (rostoucí populace). pro λ > 0. To znamená, že na intervalu (0, ∞) funkce klesá od nekonečna k nule, takže rovnice (2.41) má jediné kladné řešení, označme ho λ1. Hodnota λ1 je dominantní vlastní hodnotou matice A, tedy Malthusovským koeficientem růstu populace. Pokud R0 = f(1) > 1, pak λ1 > 1; pokud R0 = f(1) < 1, pak λ1 < 1. Situace je znázorněna na obrázku 2.2. Můžeme tedy formulovat závěr: Je-li R0 = k p=1 Fp p−1 q=1 Pq > 1, pak populace roste, je-li R0 < 1, pak populace vymírá. Pokud R0 = 1 a populace je strukturně stabilizovaná, pak se její velikost nemění. Tento výsledek, k němuž jsme dospěli s využitím Perronovy-Frobeniovy teorie, je stejný, jako závěr pravděpodobnostní úvahy provedené v oddíle 2.6.1. Délka generace γ je definována jako doba, po jejímž uplynutí jsou rodiče vystřídáni potomky stejně starými, jako byli rodiče při jejich narození. Tedy poměr velikosti generace potomků a generace rodičů je roven R0. Tento poměr je však u strukturně stabilizované populace s růstovým koeficientem λ1 roven hodnotě λγ 1. Z rovnosti R0 = λγ 1 pro délku generace dostaneme vyjádření γ = ln R0 ln λ1 . 2.6.4 Stabilizovaná věková struktura Stabilizovanou věkovou strukturu populace, tj. vlastní vektor w = (w1, w2, . . . , wk)T příslušný k dominantní vlastní hodnotě λ1 matice A, získáme řešením homogenní soustavy lineárních rovnic          F1 F2 F3 . . . Fk−1 Fk P1 0 0 . . . 0 0 0 P2 0 . . . 0 0 ... ... ... ... ... ... 0 0 0 . . . 0 0 0 0 0 . . . Pk−1 0                   w1 w2 w3 ... wk−1 wk          = λ1          w1 w2 w3 ... wk−1 wk          . 66 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Poněvadž dim ker(A−λ1I) = 1, musí být jedna z rovnic lineární kombinací ostatních. Druhá až k-tá rovnice této soustavy jsou P1w1 = λ1w2, P2w2 = λ1w3, P3w3 = λ1w4, ... Pk−1wk−1 = λ1wk. Tyto rovnice jsou lineárně nezávislé a jejich řešení je w2 = w1 λ1 P1 w3 = w2 λ1 P2 = w1 λ2 1 P1P2 w4 = w3 λ1 P3 = w1 λ3 1 P1P2P3 ... wj = w1λ1−j 1 j−1 q=1 Pq ... wk = w1λ1−k 1 k−1 q=1 Pq. Odtud je vidět, že pokud λ1 ≥ 1, pak w1 ≥ w2 ≥ · · · ≥ wk. Tyto nerovnosti jsou nutnou podmínkou k tomu, aby strukturně stabilizovaná populace nevymírala. Dostáváme tak závěr: je-li některá věková třída ve strukturně stabilizované populaci početnější než věková třída mladších jedinců, pak populace vymírá. Z požadavku, aby složky vektoru w vyjadřovaly relativní zastoupení věkových tříd ve strukturně stabilizované populaci, tj. w 1 = 1, plyne podmínka 1 = k p=1 wp = w1 k p=1 λ1−p 1 p−1 q=1 Pq, tedy wj = λ1−j 1 j−1 q=1 Pq k p=1 λ1−p 1 p−1 q=1 Pq = j−1 q=1 Pq k p=1 λj−p 1 p−1 q=1 Pq . 2.6.5 Reproduktivní hodnota věkových tříd Reproduktivní hodnotu věkových tříd, tj. levý vlastní vektor v = (v1, v2, . . . , vk)T příslušný k dominantní vlastní hodnotě λ1 matice A, získáme řešením homogenní soustavy lineárních rovnic ATv = λ1v:          F1 P1 0 . . . 0 0 F2 0 P2 . . . 0 0 F3 0 0 . . . 0 0 ... ... ... ... ... ... Fk−1 0 0 . . . 0 Pk−1 Fk 0 0 . . . 0 0                   v1 v2 v3 ... vk−1 vk          = λ1          v1 v2 v3 ... vk−1 vk          . (2.42) 2.6. ANALÝZA VĚKOVĚ STRUKTUROVANÉ POPULACE 67 Tuto soustavu přepíšeme na tvar F1v1 + P1v2 = λ1v1, F2v1 + P2v3 = λ1v2, ... Fk−2v1 + Pk−2vk−1 = λ1vk−2, Fk−1v1 + Pk−1vk = λ1vk−1, Fkv1 = λ1vk. Z poslední až druhé rovnice postupně vypočítáme vk = Fk λ1 v1, vk−1 = 1 λ1 (Fk−1v1 + Pk−1vk) = Fk−1λ−1 1 + FkPk−1λ−2 1 v1, vk−2 = 1 λ1 (Fk−2v1 + Pk−2vk−1) = Fk−2λ−1 1 + Fk−1Pk−2λ−2 1 + FkPk−1Pk−2λ−3 1 v1 = = k p=k−2  Fp p−1 q=k−2 Pq   λk−3−p 1 v1, ... vi = k p=i  Fp p−1 q=i Pq   λi−1−p 1 v1 ... v2 = k p=2  Fp p−1 q=2 Pq   λ1−p 1 v1. Dosazením vypočítaného v2 do první rovnice dostaneme v1 = 1 λ1  F1 + P1 k p=2  Fp p−1 q=2 Pq   λ1−p 1   v1 = v1 k p=1  Fp p−1 q=1 Pq   λ−p 1 ; poněvadž vlastní hodnota λ1 splňuje charakteristickou rovnici (2.41), vidíme, že první rovnice soustavy (2.42) je závislá na druhé až k-té (což odpovídá tomu, že dim ker(A − λ1I) = 1). Bývá vhodné volit v1 = 1, tj. vyjadřovat reproduktivní hodnotu věkové třídy relativně k reproduktivní hodnotě novorozenců. V takovém případě je vi = k p=i  Fp p−1 q=i Pq   λi−1−p 1 = k p=i  Fp p−1 q=1 Pq   λi−1−p 1 i−1 q=1 Pq . Porovnáním s (2.37) vidíme, že reproduktivní hodnota i-té věkové třídy je součtem fertilitních funkcí do nejvzdálenějšího možného konce života jedinců této věkové třídy „diskontovaných růstovým koeficientem a pravděpodobností přežívání. 68 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Předpokládejme nyní, že jedinci z uvažované populace jsou plodní až od jistého věku, tj. že existuje index m, 1 < m < k, takový, že F1 = F2 = · · · = Fm = 0 < Fm+1. Dále předpokládejme, že populace nevymírá, tj. λ1 ≥ 1. V takovém případě pro i ≤ m platí vi = k p=m+1  Fp p−1 q=i Pq   λi−1−p 1 v1 < k p=m+1  Fp p−1 q=i+1 Pq   λi−p 1 v1 = vi+1. Z tohoto výsledku můžeme zformulovat závěr: V nevymírající populaci se stabilizovanou věkovou strukturou reproduktivní hodnota nedospělých jedinců s věkem roste. 2.6.6 Citlivost růstového koeficientu na plodnost a přežívání Nechť λ1 je dominantní vlastní hodnota Leslieho matice A (kladné řešení charakteristické rovnice (2.41)), w a v příslušný pravý a levý vlastní vektor. Podle rovnosti (2.35) platí ∂λ1 ∂Fj = v1wj vTw , ∂λ1 ∂Pj = vj+1wj vTw . (2.43) Z první rovnosti a s využitím výsledků pododdílu 2.6.4 dostaneme ∂λ1 ∂Fj ∂λ1 ∂Fj+1 = wj wj+1 = λ1−j 1 j−1 q=1 Pq λ−j 1 j q=1 Pq = λ1 Pj . (2.44) Pokud λ1 ≥ 1 (populace nevymírá), pak je pravá strana poslední rovnosti větší než 1, tedy ∂λ1 ∂Fj > ∂λ1 ∂Fj+1 . V nevymírající populaci citlivost růstového koeficientu na plodnost klesá s věkem. S využitím rovnosti (2.44) a výsledků pododdílu 2.6.5 z druhé rovnosti (2.43) plyne ∂λ1 ∂Pj ∂λ1 ∂Pj+1 = vj+1wj vj+2wj+1 = vj+1 vj+2 λ1 Pj = k p=j+1 Fp p−1 q=j+1 Pq λj−p 1 k p=j+2 Fp p−1 q=j+2 Pq λj+1−p 1 λ1 Pj = = Fj+1 + k p=j+2 Fp p−1 q=j+1 Pq λj−p 1 k p=j+2 Fp p−1 q=j+2 Pq λj+1−p 1 λ1 Pj = Fj+1 Pj+1 + k p=j+2 Fp p−1 q=j+2 Pq λj−p 1 k p=j+2 Fp p−1 q=j+2 Pq λj−p 1 Pj+1 Pj ≥ Pj+1 Pj , (2.45) v případě Fj+1 > 0 je poslední nerovnost ostrá. Předpokládejme, že pravděpodobobnost přežití nedospělých jedinců roste s věkem, tj. největší úmrtnost mají novorozenci a úmrtnost s věkem až do dosažení plodnosti klesá, takže P1 < P2 < · · · < Pm. V takovém případě je ∂λ1 ∂Pj > ∂λ1 ∂Pj+1 pro j = 1, 2, . . . , m − 1. 2.6. ANALÝZA VĚKOVĚ STRUKTUROVANÉ POPULACE 69 Pokud úmrtnost nedospělých jedinců klesá s věkem, pak citlivost růstového koeficientu na pravděpodobnost přežití u nedospělých jedinců s věkem klesá. Poněvadž λ1 > 0, lze nerovnost (2.45) přepsat na tvar Pj λ1 ∂λ1 ∂Pj ≥ Pj+1 λ1 ∂λ1 ∂Pj+1 , rovnost nastane právě tehdy, když Fj+1 = 0. Porovnáním s definicí pružnosti v oddílu 2.5 vidíme, že pružnost růstového koeficientu vzhledem k pravděpodobobnosti přežití s věkem neroste, u nedospělých jedinců tato pružnost na věku nezávisí. Uvažujme ještě jeden důsledek rovností (2.43), a to ∂λ1 ∂Pj ∂λ1 ∂Fj = vj+1wj v1wj = vj+1 v1 . Růstový koeficient je citlivější na přežívání nějaké věkové třídy než na její plodnost právě tehdy, když reproduktivní hodnota následující věkové třídy je větší než reproduktivní hodnota novorozenců. Volíme-li v1 = 1, lze uvést další interpretaci reproduktivní hodnoty: reproduktivní hodnota věkové třídy vyjadřuje poměr citlivosti růstového koeficientu na přežívání a plodnost věkové třídy předchozí. 2.6.7 Očekávaný věk při úmrtí Uvažujme populaci, která je strukturně stabilizovaná, tj. populaci, jejíž vyvoj je popsán rov- ností n(t) = cλ1w, kde λ1 je dominantní vlastní hodnota matice A a w je příslušný vlastní vektor. Nechť V je náhodná veličina vyjadřující věk nějakého jedince z populace, který zemřel během časového intervalu (t, t + 1). Počet všech jedinců, kteří měli v čase t věk j − 1 a zemřeli během intervalu (t, t + 1), je podle výsledků pododdílu 2.6.4 roven nj(t) − nj+1(t + 1) = cλt 1wj − cλt+1 1 wj+1 = cλt 1  λ1−j 1 j−1 q=1 Pq − λ1λ 1−(j+1) 1 j q=1 Pq   w1 = = cw1λt+1 1 λ−j 1 (1 − Pj) j−1 q=1 Pq pro j = 1, 2, . . . , k −1. Počet všech jedinců, kteří měli v čase t věk k −1, a tedy všichni během intervalu (t, t + 1) zemřeli, je nk(t) = cλt 1wk = cw1λt+1−k 1 k−1 q=1 Pq. 70 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Při zavedené konvenci Pk = 0 je počet všech jedinců, kteří měli v čase t věk j − 1 a zemřeli během časového intervalu (t, t + 1) roven cw1λt+1 1 λ−j 1 (1 − Pj) j−1 q=1 Pq pro j = 1, 2, . . . , k. Počet všech jedinců, kteří zemřeli během intervalu (t, t + 1), je tedy roven k p=1 cw1λt+1 1 λ−p 1 (1 − Pp) p−1 q=1 Pq = cw1λt+1 1 k p=1 λ−p 1 (1 − Pp) p−1 q=1 Pq. Považujeme-li pravděpodobnost za klasickou, můžeme psát P(V = j − 1) = cw1λt+1 1 λ−j 1 (1 − Pj) j−1 q=1 Pq cw1λt+1 1 k p=1 λ−p 1 (1 − Pp) p−1 q=1 Pq = λ−j 1 (1 − Pj) j−1 q=1 Pq k p=1 λ−p 1 (1 − Pp) p−1 q=1 Pq . Pravděpodobnost, že jedinec ze strukturně stabilizované populace uhynulý během projekčního inervalu má určitý věk tedy nezávisí na čase. Střední hodnota věku při úmrtí je E V = k j=1 (j − 1)λ−j 1 (1 − Pj) j−1 q=1 Pq k p=1 λ−p 1 (1 − Pp) p−1 q=1 Pq = k j=1 jλ−j 1 (1 − Pj) j−1 q=1 Pq k j=1 λ−j 1 (1 − Pj) j−1 q=1 Pq − 1. (2.46) Pro zjednodušení zápisu označme na chvíli αj = (1−Pj) j−1 q=1 Pq, βj = j λ−j 1 αj , γj = λ−j 1 αj . Pak E V = k j=1 jλ−j 1 αj k j=1 λ−j 1 αj − 1, takže ∂ E V ∂λ1 = − k j=1 j2λ−j−1 1 αj k j=1 λ−j 1 αj − k j=1 jλ−j 1 αj − k j=1 jλ−j−1 1 αj k j=1 λ−j 1 αj 2 = = k j=1 jλ−j 1 αj 2 − k j=1 j2λ−j 1 αj k j=1 λ−j 1 αj λ1 k j=1 λ−j 1 αj 2 = = k j=1 βjγj 2 − k j=1 β2 j k j=1 γ2 j λ1 k j=1 γ2 j 2 < 0 2.7. UDÁLOSTI V ŽIVOTNÍM CYKLU 71 podle Cauchyovy-Buňakovského-Schwartzovy nerovnosti. Tedy s klesajícím růstovým koeficientem roste očekávaný věk při úmrtí. Upravme nyní jmenovatel zlomku v rovnosti (2.46): k j=1 λ−j 1 (1 − Pj) j−1 q=1 Pq = k j=1 λ−j 1 j−1 q=1 Pq − k j=1 λ−j 1 j q=1 Pq = k j=1 λ−j 1 j−1 q=1 Pq − k−1 j=1 λ−j 1 j q=1 Pq = = 1 λ1 + k j=2 λ−j 1 j−1 q=1 Pq − k j=2 λ1−j 1 j−1 q=1 Pq = 1 λ1 + k j=2 λ−j 1 (1 − λ1) j−1 q=1 Pq. Očekávaný věk při úmrtí tedy můžeme vyjádřit jako E V = λ1 k j=1 (j − 1)λ−j 1 (1 − Pj) j−1 q=1 Pq 1 + k j=2 λ1−j 1 (1 − λ1) j−1 q=1 Pq . (2.47) Poněvadž E V je klesající funkcí proměnné λ1, porovnáním výrazů (2.39) a (2.47) vidíme, že E V = E T právě tehdy, když λ1 = 1. Očekávaný věk při úmrtí a střední délka života jsou stejné jedině v populaci se stabilizovanou velikostí. V rostoucí populaci je střední věk při úmrtí nižší než střední délka života, ve vymírající populaci naopak vyšší. 2.7 Události v životním cyklu Za životní cyklus jedince považujeme období od narození (vylíhnutí, vyklíčení) do smrti. Události, které ho v té době potkají, mohou být dosažení dospělosti, přeměna do dalšího vývojového stadia, vykvetení, ztráta plodnosti, emigrace, imigrace a podobně. Událost v životním cyklu je vyjádřena jako přechod do jiného i-stavu v průběhu projekčního intervalu, je tedy charakterizována nenulovou složkou v projekční matici A. V té jsou ovšem zahrnuty také „objevení se nových jedinců. Abychom odlišili události v životním cyklu od rození, provedeme dekompozici projekční matice, tj. vyjádříme ji ve tvaru součtu dvou matic, z nichž jedna vyjadřuje přechody mezi jednotlivými i-stavy, druhá reprodukci: A = T + F. Přitom složka tij matice T je pravděpodobnost, že jedinec z j-té třídy přejde během projekčního intervalu do třídy i-té; složka fij matice F je střední počet potomků jedince z j-té třídy, kteří se během projekčního intervalu objeví ve třídě i-té. Poznamenejme, že tříd, v nichž mohou být „noví jedinci, může být více. Mohou to být novorozenci na různých lokalitách u metapopulací, dormantní a klíčící semena u rostlin a podobně. Část populace tvořenou jedinci, kteří „se objevili (narodili, vylíhli, . . . ) ve stejném okamžiku, nazýváme kohorta. Složení kohorty vzniklé v čase t je dáno součinem Fn(t−1), kohorta se vyvíjí jako populace s projekční maticí T. Ke k třídám, do nichž je populace strukturovaná, přidáme k+1-ní třídu, v níž jsou jedinci mrtví. Zavedeme parametry mj = 1 − k i=1 tij, j = 1, 2, . . . , k 72 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ vyjadřující pravděpodobnost, že jedinec z j-té třídy během projekčního intervalu zemře. Pak lze životní cyklus interpretovat jako Markovův řetězec s maticí pravděpodobností přechodu P = T O mT 1 . O matici T je rozumné předpokládat: • Dominantní vlastní hodnota λ nezáporné matice T je menší než 1, tj. celková velikost populace, v níž „nevznikají noví jedinci se v průběhu času zmenšuje. • Ke každému indexu j ∈ {1, 2, . . . , k} existuje konečná posloupnost indexů i1, i2, . . . , is taková, že ti1,j > 0, ti2,i1 > 0, . . . , tis,is−1 > 0, mis > 0, tj. jedinec z jakékoliv třídy se v konečném čase může dostat do třídy zemřelých, neexistuje nějaká třída nesmrtelných. Poněvadž T 2 = λ < 1, platí lim t→∞ Tt 2 ≤ lim t→∞ λt = 0, lim t→∞ t p=0 Tp 2 ≤ lim t→∞ t p=0 T p 2 = lim t→∞ t p=0 λp = lim t→∞ 1 − λt+1 1 − λ = 1 1 − λ , takže lim t→∞ Tt = O a řada ∞ t=0 Tt konverguje. Označme dále N = ∞ t=0 Tt = (I − T)−1 ; matice N se nazývá fundamentální matice uvažovaného Markovova řetězce. Poněvadž P2 = T O mT 1 T O mT 1 = T2 O mT(T + I) 1 , P3 = T2 O mT(T + I) 1 T O mT 1 = T3 O mT(T2 + T + I) 1 , . . . . . . , Pt =   Tt O mT t−1 p=0 Tp 1   , vidíme, že (Tt)ij vyjadřuje pravděpodobnost, že jedinec, který byl na počátku ve třídě j bude v čase t ve třídě i. Dále z uvedeného výpočtu plyne, že lim t→∞ Pt = O O mTN 1 , takže stav k + 1 je jediným absorbujícím stavem Markovova řetězce. 2.7. UDÁLOSTI V ŽIVOTNÍM CYKLU 73 2.7.1 Čas strávený v jedné třídě Uvažujme populaci, která je strukturovaná do k tříd. Předpokládejme, že do jednotlivých tříd vstupují jedinci na začátku projekčního intervalu. Pokud považujeme délku projekčního intervalu za jednotkový čas, lze celkový čas, po který je jedinec příslušníkem i-té třídy vyjádřit jako celkový počet přechodů z nějaké třídy j-té do třídy i-té; přitom připouštíme i možnost j = i, tedy setrvání jedince v i-té třídě. Zavedeme náhodnou veličinu µij(t), který indikuje, zda jedinec, který byl na počátku ve třídě j, je v čase t ve třídě i, tj. µij(t) = 1, jedinec, který byl v čase 0 ve třídě j, je v čase t ve třídě i, 0, jinak. Její střední hodnota je E µij(t) = 1(Tt )ij + 0 1 − (Tt )ij = (Tt )ij. Náhodná veličina νij vyjadřující celkový čas, po který je ve třídě i jedinec, který byl na počátku ve třídě j, je dána součtem všech jeho „pobytů ve třídě j , νij = ∞ t=0 µij(t). Očekávaný čas, po který je ve třídě i jedinec, který byl na počátku ve třídě j je tedy E νij = E ∞ t=0 µij(t) = ∞ t=0 E µij(t) = ∞ t=0 (Tt )ij = (N)ij. Dostáváme tak interpretaci fundamentální matice: N = (E νij)k i,j=1. 2.7.2 Očekávaná doba dožití Za dobu dožití jedince ze třídy j považujeme čas, za který se jedinec z této třídy dostane do absorpční třídy k + 1 mrtvých jedinců. Označme jako τj náhodnou veličinu, která tento čas vyjadřuje. Pak pravděpodobnost P(τj > t), že jedinec, který byl na počátku ve třídě j, bude v čase t ještě naživu, je vlastně pravděpodobností jevu, že tento jedinec nebude ve třídě k +1, tj. že bude v nějaké jiné třídě. Tato pravděpodobnost je dána součtem P(τj > t) = k i=1 (Tt )ij. Pravděpodobnostní funkci náhodné veličiny τj lze vyjádřit jako P(τj = t) = P(τj > t − 1) − P(τj > t) = k i=1 (Tt−1 − Tt )ij, neboť podle principu inkluze a exkluze platí 1 = P(τj > t − 1 ∨ τj ≤ t) = P(τj > t − 1) + P(τj ≤ t) − P(τj > t − 1 ∧ τj ≤ t) = = P(τj > t − 1) + 1 − P(τj > t) − P(τj = t). 74 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Tuto pravděpodobnostní funkci můžeme také přepsat ve tvaru P(τj = t) = k i=1 (Tt−1 − Tt )ij = 1T (Tt−1 − Tt ) T j = (Tt−1 − Tt )T 1 j . Ze známé pravděpodobností funkce můžeme vypočítat střední hodnotu i rozptyl náhodné veličiny. Očekávaná doba dožití (specifická střední délka života) jedince j-té třídy je střední hodnota náhodné veličiny τj, kterou můžeme vyjádřit jako E τj = ∞ t=0 t P(τj = t) = ∞ t=1 t k i=1 (Tt−1 − Tt )ij = k i=1 ∞ t=1 tTt−1 − ∞ t=1 tTt ij = = k i=1 ∞ t=0 (t + 1)Tt − tTt ij = k i=1 ∞ t=0 (Tt )ij = k i=1 (I − T)−1 ij = (NT 1)j. Vektor specifických středních délek života tedy je (E τ1, E τ2, . . . , E τk)T = NT1, neboli (E τ1, E τ2, . . . , E τk) = 1T N. Nyní vypočítáme rozptyl náhodné veličiny τj. Poněvadž N = ∞ t=0 Tt = ∞ t=0 (t + 1)Tt − tTt = ∞ t=1 tTt−1 − tTt = ∞ t=1 tTt−1 (I − T) = ∞ t=1 tTt−1 N−1 , můžeme vyjádřit N2 = ∞ t=1 tTt−1 a dále ∞ t=0 tTt = ∞ t=1 (tTt−1 − tTt−1 + tTt ) = N2 − ∞ t=1 tTt−1 (I − T) = N2 − ∞ t=1 tTt−1 N−1 = N2 − N, takže E τ2 j = ∞ t=0 t2 P(τj = t) = ∞ t=0 t2 k i=1 (Tt−1 − Tt )ij = k i=1 ∞ t=0 t2 (Tt−1 − Tt ) ij = = k i=1 ∞ t=1 t2 Tt−1 − ∞ t=0 t2 Tt = k i=1 ∞ t=0 (t + 1)2 Tt − t2 Tt ij = = k i=1 ∞ t=0 (2t + 1)Tt ij = k i=1 2(N2 − N) + N ij = k i=1 (2N2 − N)ij = (2N2 − N)T 1 j . Pro rozptyl náhodné veličiny τj tedy dostáváme var τj = E τ2 j − (E τj)2 = (2N2 − N)T 1 j − (NT 1)j 2 a řádkový vektor rozptylů dob dožití můžeme psát ve tvaru (var τ1, var τ2, . . . , var τk) = 1T (2N2 − N) − 1T N ◦ 1T N. Pokud jsou ve třídě i novorozenci (tedy jedinci věku 0), je očekávaný věk při úmrtí těchto jedinců roven E τi − 1. 2.7. UDÁLOSTI V ŽIVOTNÍM CYKLU 75 2.7.3 Věkově specifická plodnost Připomeňme, že složka (Ta)lj vyjadřuje pravděpodobnost, že jedinec, který byl na počátku (tj. v čase 0) ve třídě j bude v čase a ve třídě l; je-li l = k +1, je to pravděpodobnost, že tento jedinec je ve třídě l a žije. To znamená, že k l=1 (Ta)lj vyjadřuje pravděpodobnost, že jedinec, který byl na počátku ve třídě j bude v čase a ještě naživu, tj. v nějaké jiné třídě než k + 1-ní. Označme nyní σlj(a) podmíněnou pravděpodobnost, že jedinec, který byl na počátku ve třídě j bude v čase a ve třídě l za podmínky, že dosud žije. Tato pravděpodobnost je dána vztahem σlj(a) = (Ta)lj k l=1 (Ta)lj = Ta (diag 1T Ta )−1 lj . Nechť i-tá třída je třídou „nově vzniklých jedinců (novorozenců, semen ap.), tj. fil = 0 pro nějaký index l. Uvažujme jedince, který byl na počátku ve třídě j, dožije se věku o a většího a v něm bude ve třídě l. Očekávaný počet jeho potomků, kterým přispěje v čase a do třídy i je dán součinem filσlj(a). Očekávaný počet potomků ve třídě i všech jedinců, kteří byli na počátku ve třídě j a v čase a jsou ještě naživu, je ϕij(a) = k l=1 filσlj(a) = FTa (diag 1T Ta )−1 ij ; ϕij(a) se nazývá věkově specifická plodnost (age-specific fertility) jedinců třídy j. Věkově specifické plodnosti můžeme zapsat do matice Φ(a) = ϕij(a) k i,j=1 = FΣ(a) = FTa (diag 1T Ta )−1 , kde Σ(a) = σij(a) k i,j=1 . Pokud je populace strukturovaná podle věku a jediná třída novorozenců je první, pak ϕ11(a) je věkově specifická plodnost jedince ve věku a v obvyklém demografickém smyslu. 2.7.4 Čistá míra reprodukce a délka generace Čistá míra reprodukce (net reproductive rate) R0 je definována jako očekávaný počet potomků jedince během jeho života. Střední hodnota doby, kterou jedinec, který byl na počátku ve třídě j, stráví ve třídě l je E νlj = (N)lj. Očekávaný počet nových jedinců „vzniklých ve třídě i, které vyprodukuje během doby strávené ve třídě l jedinec, který byl na počátku ve třídě j, je dán součinem fil E νlj. Očekávaný počet všech potomků jedince, který byl na počátku ve třídě j, objevivších se ve třídě i je tedy k l=1 fil E νlj = (FN)ij. Pokud tedy je první třída jedinou třídou, v níž jsou novorozenci (obecně: nově vzniklí jedinci), můžeme položit R0 = (FN)11. 76 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Pokud ale jsou novorozenci ve více třídách, je situace komplikovanější. Označme y(0) složení populace na počátku. Její vývoj jako generace, tj. bez potomků, je projektován maticí T, y(t + 1) = Ty(t), kde y(t) označuje složení generace v čase t. Tedy y(t) = Tt y(0). Očekávané složení potomků uvažované generace v čase t je Fy(t), takže všichni očekávaní potomci iniciální generace, tj. následující generace, má složení ∞ t=0 Fy(t) = ∞ t=0 FTt y(0) = F ∞ t=0 Tt y(0) = FNy(0). Tento výsledek lze interpretovat tak, že matice FN projektuje jednu generaci na následující. To nás opravňuje vzít za čistou míru reprodukce R0 dominantní vlastní hodnotu matice FN. Délka generace T je definována jako doba, po jejímž uplynutí je populace právě R0-krát větší, tj. n(T) 1 = R0 n(0) 1 . Má-li populace stabilizovanou strukturu, tj. n(t) = λt 1n(0), kde λ1 je dominantní vlastní číslo projekční matice A a vektor n(0) je úměrný příslušnému vlastnímu vektoru, pak n(T) 1 = λT 1 n(0) 1, tedy R0 = λT 1 , neboli T = log R0 log λ1 . 2.7.5 Rozložení věku v jednotlivých třídách stabilizované populace Uvažujme populaci, která se vyvíjí dostatečně dlouho, takže její složení v čase t je dáno vztahem n(t) = αλt 1w, kde λ1 je dominantní vlastní číslo projekční matice A, vektor w je příslušný vlastní vektor a α > 0 je vhodná konstanta. Složení kohorty vzniklé v čase t je tedy Fn(t − 1) = αλt−1 1 Fw a složení kohorty vzniklé v čase t − a je αλt−a−1 1 Fw. Kohorta vzniklá v čase t − a má v čase t složení Ta αλt−a−1 1 Fw = (αλt−1 1 )λ−a 1 Ta Fw. Označme x(a) = λ−a 1 Ta Fw. Tento vektor je úměrný složení kohorty stáří a, tj. složení, jaké má v čase t kohorta vzniklá v čase t − a. V poli (matici s neomezeným počtem sloupců) X = x(0), x(1), x(2), . . . = Fw, λ−1 1 TFw, λ−2 1 T2 Fw, . . . 2.7. UDÁLOSTI V ŽIVOTNÍM CYKLU 77 odpovídají sloupce věku a řádky třídě. Tedy zlomek x(a) j x(a) 1 = λ−a 1 TaFw j λ−a 1 TaFw 1 = (TaFw)j TaFw 1 , a = 0, 1, 2, . . . , vyjadřuje relativní zastoupení jedinců třídy j mezi všemi jedinci věku a, tedy pravděpodobnostní funkci náhodné veličiny „třída, z níž je jedinec věku a ; zlomek x(a) j ∞ s=0 (x(s))j = λ−a 1 TaFw j ∞ s=0 λ−s 1 TsFw j = (TaFw)j ∞ s=0 λa−s 1 TsFw j , j = 1, 2, . . . , k, vyjadřuje relativní zastoupení jedinců věku a mezi všemi jedinci třídy j, tj. pravděpodobnostní funkci náhodné veličiny „věk jedince z j-té třídy . Střední hodnotu věku jedinců j-té třídy je nyní dána výrazem βj = ∞ a=0 a (TaFw)j ∞ s=0 λa−s 1 TsFw j = ∞ a=0 aλ−a 1 (TaFw)j ∞ a=0 λ−a 1 (TaFw)j . Připomeňme, že mj = 1− k i=1 tij vyjadřuje pravděpodobnost, že jedinec z j-té třídy během projekčního intervalu zemře. Očekávaný počet jedinců třídy j, kteří zemřou, je tedy mjnj(t) = mjαλt 1wj = αλt 1mjwj a celkový počet umírajících jedinců je k j=1 mjnj(t) = αλt 1 k j=1 mjwj = αλt 1mT w. To znamená, že pravděpodobnost, že umírající jedinec je z třídy j je dána podílem mjwj mTw . Průměrný věk γ jedinců umírajících během projekčního intervalu nyní vyjádříme jako vážený průměr středních hodnot věků umírajících jedinců ve všech třídách, γ = k j=1 mjwj mTw βj = 1 mTw k j=1 mjwj ∞ a=0 aλ−a 1 (TaFw)j ∞ a=0 λ−a 1 (TaFw)j . 78 KAPITOLA 2. MODELY S KONSTANTNÍ PROJEKČNÍ MATICÍ Kapitola 3 Identifikace parametrů modelu 3.1 Inversní metody časových řad V tomto oddílu budeme předpokládat, že z pozorování nebo experimentu známe strukturu populace v T + 1 po sobě následujících časových okamžicích, tedy že máme vektory n(0), n(1), n(2), . . . , n(T). Přitom je populace složena z k tříd, tj. vektory n(t), t = 0, 1, . . . , T jsou k-rozměrné. U prvních dvou metod — regresní a maximální věrohodnosti — obecně nebudeme nic předpokládat o projekční matici A kromě její nezápornosti. Budeme se tedy snažit z pozorovaných dat identifikovat (odhadnout) všech k2 prvků této matice. To může být užitečné například v situacích, kdy nevíme, zda některé věkové skupiny jsou plodné či nikoliv, zda některé stadium může přežít projekční interval a podobně. Těmito metodami pak ale odhadujeme i parametry, které musí být nulové, například ty, které vyjadřují nemožný přechod mezi stadii; z kukly již nemůže vzniknout vajíčko. To zvyšuje nároky na výpočetní kapacitu, zanáší do výpočtu chyby a také může vytvářet obtížně řešitelné problémy při interpretaci výsledků. Třetí metoda využívající kvadratického programování naopak vyžaduje znalost tvaru projekční matice. Je tedy potřebné „a priori vědět, které prvky projekční matice vyjadřují pravděpodobnost možných jevů a jsou tedy z intervalu (0, 1], které vyjadřují plodnost a mohou být větší než 1 a podobně. 3.1.1 Regresní metody Budeme předpokládat, že pozorované složení populace v časovém okamžiku t + 1 je projekcí jejího složení v okamžiku t, a navíc se na něm projevují nějaké náhodné vlivy nebo chyby pozorování. Tedy ni(t + 1) = k j=1 aijnj(t) + εi(t), i = 1, 2, . . . , k, t = 0, 1, 2, . . . , T − 1. 79 80 KAPITOLA 3. IDENTIFIKACE PARAMETRŮ MODELU Hodnota εi(t) je přitom realizací nějaké náhodné veličiny; je rozumné předpokládat, že její střední hodnota je 0. Předchozí rovnosti můžeme pro každé i přepsat maticově,      ni(1) ni(2) ... ni(T)      =      n1(0) n2(0) . . . nk(0) n1(1) n2(1) . . . nk(1) ... ... ... ... n1(T − 1) n2(T − 1) . . . nk(T − 1)           ai1 ai2 ... aik      +      εi(0) εi(1) ... εi(T − 1)      , i = 1, 2, . . . , k. Označíme N =      n1(0) n2(0) . . . nk(0) n1(1) n2(1) . . . nk(1) ... ... ... ... n1(T − 1) n2(T − 1) . . . nk(T − 1)      a přepíšeme všechny rovnosti jako jednu rovnost maticovou,                     n1(1) ... n1(T) n2(1) ... n2(T) ... nk(1) ... nk(T)                     =      N O . . . O O N . . . O ... ... ... ... O O . . . N                          a11 ... a1k a21 ... a2k ... ak1 ... akk                     +                     ε1(0) ... ε1(T − 1) ε2(0) ... ε2(T − 1) ... εk(0) ... εk(T − 1)                     . Vektor na levé straně rovnosti označíme y, vektor chyb (vektor za znakem +) označíme ε. Předchozí rovnost nyní můžeme přepsat do tvaru y = (I ⊗ N) vec(AT ) + ε, nebo při označení X = I ⊗ N, β = vec(AT) ještě stručněji y = Xβ + ε. (3.1) Vektor y a matici X známe z pozorování. Vektor β je vektorem neznámých parametrů, složek projekční matice A, který chceme identifikovat. Tvar rovnosti sugeruje, že vektor parametrů β bychom mohli odhadnout metodami lineární regrese. „Klasickou metodou nejmenších čtverců tak dostaneme odhad parametrů ve tvaru ˆβ = (XT X)−1 XT y. (3.2) Tato formule byla ovšem odvozena za předpokladu, že nezávisle proměnné (tj. složky „matice plánu X) jsou nenáhodné veličiny (nejsou zatíženy chybou). To však v rovnosti (3.1) neplatí, matice X obsahuje tytéž složky, které jsou také složkami vektoru y. Z tohoto důvodu 3.1. INVERSNÍ METODY ČASOVÝCH ŘAD 81 není korektní parametry β odhadovat výrazem na pravé straně rovnosti (3.2), ale metodami orthogonální regrese (total least squares). Další potíž spočívá v tom, že některé parametry mohou vyjít jako záporné, nebo že parametry, které by měly vyjadřovat pravděpodobnosti, mohou vyjít větší než 1. V takovém případě nahradíme nerealistické hodnoty nulami, respektive jedničkami. 3.1.2 Metoda maximální věrohodnosti Stejně jako u regresních metod budeme předpokládat, že pozorované složení populace v čase t + 1 je projekcí jejího složení v čase t a náhodné odchylky. Nyní však budeme předpokládat, že náhodná odchylka je multiplikativní, ni(t + 1) = eδi(t) k j=1 aijnj(t), i = 1, 2, . . . , k, t = 0, 1, 2, . . . , T − 1. (3.3) O chybách budeme předpokládat, že v jednom každém časovém okamžiku jsou realizací náhodného vektoru z k-rozměrného normálního rozdělení se střední hodnotou o a varianční maticí Σ, δ(t) =      δ1(t) δ2(t) ... δk(t)      ∼ N(o, Σ) (3.4) a že jsou v jednotlivých časových okamžicích nezávislé, tj. δ(t1) a δ(t2) jsou nezávislé náhodné vektory pro t1 = t2. Poznamenejme, že varianční matice Σ nemusí být diagonální; např. podmínky, které jsou dobré pro mladé jedince, mohou být dobré i pro staré nebo naopak. Dále budeme předpokládat, že všechny pozorované hodnoty jsou kladné. Můžeme je tedy zlogaritmovat, tj. položit mi(t) = ln ni(t), i = 1, 2, . . . , k, t = 0, 1, 2, . . . , T, neboli m(t) = ln n(t), t = 0, 1, 2, . . . , T. Poněvadž podle (3.3) platí mi(t) = ln ni(t) = δi(t − 1) + ln k j=1 aijnj(t − 1), je při daných hodnotách n(t − 1) vektor m(t) realizací náhodného vektoru z k-rozměrného normálního rozdělení se střední hodnotou µ(t), µ(t) i = ln k j=1 aijnj(t − 1), t = 1, 2, . . . , T (3.5) 82 KAPITOLA 3. IDENTIFIKACE PARAMETRŮ MODELU a varianční maticí Σ. Z předpokládané nezávislosti chyb v různých časových okamžicích nyní plyne, že věrohodnostní funkce je tvaru L(A, Σ) = P m(1), m(2), . . . , m(T)|A, Σ, m(0) = T t=1 P m(t)|A, Σ, m(t − 1) = = T t=1 1 (2π)k det Σ exp − 1 2 m(t) − µ(t) T Σ−1 m(t) − µ(t) . Odtud dostaneme − ln L(A, Σ) = T t=1 1 2 ln(2π)k + ln det Σ + m(t) − µ(t) T Σ−1 m(t) − µ(t) = = Tk 2 ln 2π + 1 2 − T ln det Σ−1 + T t=1 m(t) − µ(t) T Σ−1 m(t) − µ(t) . Poněvadž první člen, tj. výraz 1 2Tk ln 2π, nezávisí na matici A ani na matici Σ, můžeme maximálně věrohodný odhad parametrů, tj. složek matic A a Σ−1, vypočítat jako (ˆA, ˆΣ−1 ) = argmin T t=1 m(t) − µ(t) T Σ−1 m(t) − µ(t) , kde složky vektorů µ(t), t = 1, 2, . . . , T jsou dány rovnostmi (3.5) (a závisí tedy na matici A). Minimum hledáme nějakou iterační metodou, jako výchozí aproximaci můžeme použít odhad (3.2). Pokud jsou pozorované hodnoty n(t), t = 0, 1, . . . , T získány opakovaným měřením, můžeme spočítat jejich výběrový rozptyl, varianční matici Σ nahradit výběrovou kovarianční maticí S a odhadovat pouze složky matice A, ˆA = argmin T t=1 m(t) − µ(t) T S−1 m(t) − µ(t) . Ještě poznamenejme, že předpoklad o nezávislosti náhodných odchylek od deterministického modelu v jednotlivých časových okamžicích je dost silný, ve skutečné populaci nemusí být splněn; např. pokud byla populace v jednom časovém okamžiku díky příznivým podmínkám v dobrém fyziologickém stavu, může být její růst do dalšího okamžiku větší než obvykle, i když podmínky se nezávisle změní na nějaké méně příznivé. 3.1.3 Metoda kvadratického programování Tuto metodu nejprve ukážeme na konkrétním příkladu. Uvažujme populaci strukturovanou do tří tříd (stadií), přičemž v první třídě jsou novorozenci a ve třetí jsou plodní jedinci. Vývoj populace je tedy popsán projekční rovnicí   n1 n2 n3   (t + 1) =   0 0 F P1 Q2 0 0 P2 Q3     n1 n2 n3   (t). (3.6) 3.1. INVERSNÍ METODY ČASOVÝCH ŘAD 83 Tuto rovnici můžeme přepsat v jednotlivých složkách jako systém n1(t + 1) = F n3(t) = n3(t)F n2(t + 1) = P1n1(t)+Q2n2(t) = n1(t)P1+n2(t)Q2 n3(t + 1) = P2n2(t)+Q3n3(t) = n2(t)P2 +n3(t)Q3 nebo v jiném maticovém tvaru   n1 n2 n3   (t + 1) =   0 0 0 n3(t) 0 n1(t) n2(t) 0 0 0 0 0 n2(t) 0 n3(t)         P1 Q2 P2 F Q3       . (3.7) Matici na pravé straně této rovnice označíme M(t), vektor označíme p a rovnici zapíšeme stručně jako n(t + 1) = M(t)p, t = 0, 1, 2, . . . , T − 1. Tyto rovnice můžeme zapsat jako jednu      n(1) n(2) ... n(T)      =      M(0) M(1) ... M(T − 1)      p; 3T-rozměrný vektor na levé straně označíme z, matici typu 3T × 6 na levé straně označíme M a dostaneme z = Mp. Složky vektoru z a složky matice M jsou měřené hodnoty, vektor p je tvořen parametry, které chceme odhadnout. Pokud by se populace vyvíjela přesně podle modelu (3.6) a pozorování by nebyla zatížena chybou, platilo by podle předchozí rovnosti z−Mp = o, neboli z − Mp = 0. Proto za odhady parametrů p vezmeme takové, které minimalizují normu vektoru z − Mp; konkrétně použijeme normu euklidovskou. Parametry jsou nezáporné, P1, P2, Q2, Q3 vyjadřují pravděpodobnosti, součet pravděpodobnosti P2 přežití a přechodu ze druhé do třetí třídy a Q2 přežití ve druhé třídě nemůže převýšit hodnotu 1. Parametry tedy musí splňovat nerovnosti 0 ≤ F, 0 ≤ P1 ≤ 1, 0 ≤ P2 ≤ 1, 0 ≤ Q2 ≤ 1, 0 ≤ Q3 ≤ 1, P2 + Q2 ≤ 1. Tyto nerovnosti přepíšeme v maticovém tvaru                 0 0 0 −1 0 −1 0 0 0 0 1 0 0 0 0 0 0 −1 0 0 0 0 1 0 0 0 −1 0 0 0 0 1 0 0 0 0 0 0 0 −1 0 0 0 0 1 0 1 1 0 0                       P1 Q2 P2 F Q3       ≤                 0 0 1 0 1 0 1 0 1 1                 . (3.8) 84 KAPITOLA 3. IDENTIFIKACE PARAMETRŮ MODELU Matici na levé straně označíme C, vektor na pravé straně označíme b a dostaneme podmínky ve tvaru Cp ≤ b. Celkem tak dostáváme, že odhad parametrů p můžeme hledat tak, že najdeme minimum normy vektoru z − Mp za podmínky (3.8), stručně ˆp = argmin z − Mp 2 2 : Cp ≤ b . V obecném případě uvažujeme populaci vyvíjející se podle rovnice n(t + 1) = An(t), (3.9) která odpovídá rovnici (3.6) z úvodního příkladu a kterou můžeme přepsat ve tvaru n(t + 1) = n(t)T ⊗ I vec A, (3.10) kde I je jednotková matice řádu k (viz výpočet na str. 18). Označíme-li N(t) = n(t)T ⊗ I, můžeme tuto rovnost zapsat stručněji, n(t + 1) = N(t) vec A. (3.11) Předpokládáme, že známe strukturu matice A, takže víme, že mezi jejími složkami je právě l nenulových a zbývajících k2 − l je nulových; v úvodním příkladu bylo l = 5, k = 3 a tedy k2 − l = 4. Vektor vec A tedy obsahuje pouze l nenulových prvků. Je-li j-tá složka vektoru vec A rovna nule, pak každý prvek z j-tého sloupce matice N(t) je při násobení v předchozí rovnosti násoben nulou. To znamená, že j-tý sloupec matice N(t) k výsledku ničím nepřispívá, je zbytečný. Tato úvaha vede k tomu, že můžeme snížit dimenzi problému. Vektor vec A nahradíme vektorem p, který obsahuje nenulové prvky vektoru vec A, matici N(t) nahradíme maticí M(t), která vznikne z matice N(t) tak, že v ní vynecháme všechny sloupce, které odpovídají nulovým prvkům vektoru vec A; v úvodním příkladu se jednalo o rovnici (3.7). Rovnosti (3.11) tedy přepíšeme ve tvaru n(t + 1) = M(t)p, t = 0, 1, 2, . . . , T − 1, (3.12) nebo souhrnně      n(1) n(2) ... n(T)      =      M(0) M(1) ... M(T − 1)      p. Vektor na levé straně označíme z, matici na pravé straně označíme M a dostaneme z = Mp. (3.13) Vektor z i matice M jsou složeny z pozorovaných hodnot, vektor p je l-ticí parametrů, které chceme odhadnout. Pokud by se vývoj populace přesně řídil modelem (3.9), pak by podle rovnosti (3.13) platilo z − Mp = o. Tato úvaha vede k tomu, že za odhady parametrů p vezmeme takové hodnoty, aby norma vektoru z − Mp byla co nejmenší. 3.2. PARAMETRY POPULACE SE STABILIZOVANOU VĚKOVOU STRUKTUROU 85 Platí z − Mp 2 2 = (z − Mp)T (z − Mp) = zT z − zT Mp − pT MT z + pT MT Mp = = zT z − 2zT Mp + pT MT Mp = zT z + 2 1 2 pT MT Mp − zT Mp . Hodnota zTz nezávisí na parametrech, proto stačí minimalizovat výraz v závorce. Označíme G = MT M, q = MT z. Pak je matice G typu l × l a je symetrická. Hledáme vektor p s nezápornými složkami tak, aby 1 2 pT Gp − qT p → min . (3.14) Kromě nezápornosti musí složky vektoru p splňovat i další podmínky — pravděpodobnosti nemohou překročit hodnotu 1, součet všech pravděpodobností vyjadřujících přechod z nějakého stadia do jiných také nemůže být větší než 1 a podobně. Všechna taková omezení jsou lineární, můžeme je tedy podobně jako nerovnosti (3.8) z úvodního příkladu obecně zapsat ve tvaru Cp ≤ b, (3.15) kde C je vhodná matice a b je vhodný vektor; matice C má l sloupců, počet jejích a řádků je roven dimenzi vektoru b. Úloha (3.14), (3.15) je úlohou kvadratického programování v základním tvaru. 3.2 Parametry populace se stabilizovanou věkovou strukturou Uvažujme věkově strukturovanou populaci, tj. populaci která se vyvíjí podle Leslieho modelu po dostatečně dlouhou dobu. Taková populace má v nějakém čase t strukturu n(t) = N(t)w, kde N(t) = k j=1 nj(t) je celková velikost populace a w je normovaný vlastní vektor příslušný k dominantní vlastní hodnotě λ (růstovému koeficientu); norma v tomto případě používáme součtovou, w 1 = k j=1 wj = 1. 3.2.1 Odhad růstového koeficientu Předpokládejme, že máme změřenu velikost populace v časech ti, i = 1, 2, . . . , m, tj. že známe hodnoty N(t1), N(t2), . . . , N(tm) a všechny tyto hodnoty jsou kladné. Při stabilizované věkové struktuře platí N(ti) = N(t1)λti−t1 , 86 KAPITOLA 3. IDENTIFIKACE PARAMETRŮ MODELU tedy po zlogaritmování ln N(ti) = ti ln λ + ln N(t1) − t1 ln λ. Tuto rovnost lze považovat za zobecněný lineární regresní model a hodnotu růstového koeficientu odhadnout jako ˆλ = exp    m m i=1 ti ln N(ti) − m i=1 ti m i=1 ln N(ti) m m i=1 t2 i − m i=1 ti 2    . (3.16) Pokud bychom měli nepřerušenou časovou řadu pozorování velikosti populace, tj. znali všechny hodnoty N(0), N(1), . . . , N(T), lze výpočet odhadu růstového koeficientu zjednodušit na tvar ˆλ = exp    (T + 1) T t=0 t ln N(t) − T t=0 t T t=0 ln N(t) (T + 1) T t=0 t2 − T t=0 t 2    = = exp 6 T(T + 1)(T + 2) T t=0 (2t − T) ln N(t) . 3.2.2 Odhady pravděpodobností přežití a fertilit Předpokládejme nyní, že navíc máme změřené velikosti jednotlivých věkových tříd v jednom časovém okamžiku, tj. známe hodnoty n1(t), n2(t), . . . , nk(t) pro nějaký čas t. Z tohoto měření můžeme odhadnout složky vektoru w, ˆwi = ni(t) N(t) = ni(t) k j=1 nj(t) , i = 1, 2, . . . , k. (3.17) V populaci se stabilizovanou věkovou strukturou platí Piwi = λwi+1, tj. Pi = λ wi+1 wi , i = 1, 2, . . . , k − 1. Za odhad pravděpodobností Pi tedy můžeme vzít ˆPi = ˆλ ˆwi+1 ˆwi = ˆλ ni+1(t) ni(t) , i = 1, 2, . . . , k; (3.18) odhad ˆλ růstového koeficientu je přitom dán rovností (3.16). K odhadu věkově specifických plodností F1,F2, . . . , Fk využijeme charakteristickou rovnici (2.41) Leslieho matice. Dosadíme do ní odhady růstového koeficientu i pravděpodobností přežití, 1 = k j=1 Fj ˆλ−j j−1 q=1 ˆPq = k j=1 Fj ˆλ−j j−1 q=1 ˆλ nq+1(t) nq(t) = 1 ˆλn1(t) k j=1 Fjnj(t). 3.2. PARAMETRY POPULACE SE STABILIZOVANOU VĚKOVOU STRUKTUROU 87 Dostaneme tak rovnost ˆλn1(t) = k j=1 Fjnj(t). (3.19) Označíme Φ součet plodností a zavedeme relativní věkově specifické plodnosti vztahy fj = Fj Φ = Fj j q=1 Fj , j = 1, 2, . . . , k. Pokud bychom znali hodnoty fj například z nějaké teorie nebo z dalšího pozorování, můžeme dosadit do rovnosti (3.19), ˆλn1(t) = Φ k j=1 fjnj(t) a součet plodností odhadovat vztahem ˆΦ = ˆλn1(t) k j=1 fjnj(t) . (3.20) Nejjednodušší předpoklad o plodnostech je ten, že jsou na začátku života nulové, v plodném věku, tj. ve věku od menarche po menopauzu (ve věkových kategoriích m, m+1, . . . , M) jsou konstantní a v postreproduktivním období jsou opět nulové. Můžeme tedy předpokládat f1 = f2 = · · · = fm−1 = fM+1 = fM+2 = · · · = fk = 0, fm = fm+1 = · · · = fM = 1 M − m + 1 . Za tohoto předpokladu odhadneme sumární plodnost Φ vztahem ˆΦ = ˆλn1(t) M − m + 1 M j=m nj(t) . (3.21) Také lze například předpokládat, že plodnost od m-té věkové kategorie narůstá, dosahuje maxima Fmax v p-té věkové třídě, a pak klesne až na nulovou hodnotu po menopauze. Pokud nárůst a pokles budeme považovat za lineární, specifické plodnosti vyjádříme ve tvaru Fj =    Fmax j − m + 1 p − m + 1 , m ≤ j ≤ p, Fmax M − j + 1 M − p + 1 , p < j ≤ M, 0, jinak. V tomto případě dostaneme z rovnosti (3.19) pro maximální plodnost odhad ˆFmax = ˆλn1(t) (p − m + 1)(M − p + 1) (M − p + 1) p j=m (j − m + 1)nj(t) + (p − m + 1) M j=p+1 (M − j + 1)nj(t) . (3.22) 88 KAPITOLA 3. IDENTIFIKACE PARAMETRŮ MODELU V populaci se stabilizovanou věkovou strukturou platí k j=1 Fjn(t) = n1(t + 1) = λn1(t). Pokud tedy budeme znát hodnotu n1(t+1), tj. počet novorozenců v čase t+1, lze v rovnostech (3.20), (3.21), (3.22) výraz ˆλn1(t) nahradit výrazem n1(t + 1). Kapitola 4 Modely s externí variabilitou 4.1 Sezónní variabilita Budeme se zabývat situací, kdy životní cyklus populace je tvořen několika fázemi navazujícími na sebe v průběhu času. V jedné fázi je populace strukturována do několika tříd, přitom se počty tříd mohou v jednotlivých fázích životního cyklu lišit. Dobu trvání životního cyklu budeme považovat za jednotkovou, doby trvání jednotlivých fází mohou být různé. Nechť konkrétně je životní cyklus rozdělen na m fází. Předpokládejme, že počáteční (nultá) fáze prvního cyklu začíná v čase t = 0 a i-tá fáze prvního cyklu trvá od času τi do času τi+1, kde τ0, τ1, . . . , τm jsou reálná čísla taková, že 0 = τ0 < τ1 < τ2 < · · · < τm−1 < τm = 1. Předpokládejme dále, že v i-té fázi je populace strukturována do ki tříd, k0, k1, . . . , km−1 jsou přirozená čísla. Velikost populace v i-té fázi t-tého cyklu je tedy vyjádřena ki-rozměrným vektorem n(t + τi). Nechť nakonec nezáporná matice Bi typu ki+1 × ki projektuje velikost populace v i-té fázi na její velikost v i + 1-ní fázi (popisuje přechod populace z i-té fáze do následující), i = 0, 1, . . . , m − 2, nezáporná matice Bm−1 typu k0 × km−1 projektuje velikost populace v poslední fázi na její velikost v počáteční fázi následujícího cyklu. Vývoj populace budeme tedy modelovat rovnicemi n(t + τh+1) = Bhn(t + τh), h = 0, 1, . . . , m − 1, t = 0, 1, 2, . . . . (4.1) Podle této rovnice platí n(t + 1 + τh) = Bh−1n(t + 1 + τh−1) = = Bh−1Bh−2n(t + 1 + τh−2) = Bh−1Bh−2Bh−3n(t + 1 + τh−3) = · · · · · · = Bh−1Bh−2 · · · B0n(t + 1 + τ0) = Bh−1Bh−2 · · · B0n(t + τm) = = Bh−1Bh−2 · · · B0Bm−1n(t + τm−1) = Bh−1Bh−2 · · · B0Bm−1Bm−2n(t + τm−2) = · · · · · · = Bh−1Bh−2 · · · B0Bm−1Bm−2 · · · Bhn(t + τh). Pro h ∈ {0, 1, . . . , m − 1} nyní položíme Ah = Bh−1Bh−2 · · · B0Bm−1Bm−2 · · · Bh, zejména A0 = Bm−1Bm−2 · · · B0. 89 90 KAPITOLA 4. MODELY S EXTERNÍ VARIABILITOU Abychom zjednodušili zápis výpočtů, označíme ještě Am = A0. Každá z matic Ah je čtvercová řádu kh. Předchozí výsledek můžeme nyní zapsat ve tvaru n(t + τh + 1) = Ahn(t + τh) a z tohoto zápisu je vidět, že n(t + τh) = At hn(τh) = At hBh−1Bh−2 · · · B1B0n(0) = Bh−1Bh−2 · · · B1B0At 0n(0). (4.2) Budeme ještě používat označení Dh = Bh−1Bh−2 · · · B0Bm−1Bm−2 · · · Bh+1, h = 0, 1, . . . , m1. Matice Dh je typu kh × kh+1 a platí Ah = DhBh, Ah+1 = BhDh, h = 0, 1, . . . , m − 1. (4.3) Tvrzení 15. Všechny matice A0, A1, . . . , Am−1 mají stejné nenulové vlastní hodnoty. Důkaz: Nechť h ∈ {0, 1, 2, . . . , m − 1} je libovolné číslo. Označme pro stručnost r = kh, s = kh+1. Matice Bh je typu s × r, matice Dh je typu r × s. Poněvadž podle (4.3) je Ah O Bh O I Dh O I = DhBh O Bh O I Dh O I = DhBh DhBhDh Bh BhDh = = I Dh O I O O Bh BhDh = I Dh O I O O Bh Ah+1 a matice I Dh O I je regulární, platí I Dh O I −1 Ah O Bh O I Dh O I = O O Bh Ah+1 , což znamená, že matice P1 = Ah O Bh O a P2 = O O Bh Ah+1 jsou podobné a tedy mají stejná vlastní čísla1. Charakteristický polynom první matice je det(P1 − λI) = det Ah − λIr O Bh −λIs = det (Ah − λI) (−λ)s , charakteristický polynom druhé matice je det(P2 − λI) = det −λIr O Bh Ah+1 − λIs = (−λ)r det (Ah+1 − λIs) . To znamená, že vlastní hodnoty matice P1 jsou stejné, jako vlastní hodnoty matice Ah plus s nul; vlastní hodnoty matice P2 jsou stejné, jako vlastní hodnoty matice Ah+1 plus r nul. A poněvadž matice P1 a P2 mají stejné vlastní hodnoty, mají matice Ah a Ah+1 stejné nenulové vlastní hodnoty. 1 Matice M a N jsou podobné, pokud existuje regulární matice P taková, že P−1 MP = N. Číslo λ je vlastní hodnotou matice M právě tehdy, když M − λI = S a matice S je singulární. Vynásobením této rovnosti maticí P−1 zleva a maticí P zprava dostaneme ekvivalentní rovnost N − λI = P−1 SP. Přitom matice P−1 SP je singulární, což znamená, že λ je vlastním číslem matice N. 4.1. SEZÓNNÍ VARIABILITA 91 Tvrzení 16. Nechť všechny matice A0, A1, . . . , Am−1 jsou primitivní a λ je jejich společná dominantní vlastní hodnota. Pak pro každé h ∈ {0, 1, . . . , m − 1} platí lim t→∞ n(t + τh) λt n(τh) 1 = wh, kde wh je (pravý) vlastní vektor matice Ah příslušný k vlastní hodnotě λ. Důkaz: Tvrzení plyne z první rovnosti (4.2) a z 2.3.1. Vývoj populace směřuje ke stavu, že se její struktura (relativní zastoupení jednotlivých tříd) v jednotlivých fázích nemění; struktura populace se periodicky mění (s periodou délky populačního cyklu). Předpokládejme, že matice A0 má vlastní hodnotu λ, která je větší než absolutní hodnota všech ostatních vlastních hodnot. V takovém případě je λ společná dominantní vlastní hodnota matic Ah, h = 0, 1, . . . , m − 1, tj. λ je růstový koeficient populace. Označme wh pravý normovaný vlastní vektor matice Ah příslušný k této dominantní vlastní hodnotě λ. Pak podle (4.3) platí λwh = Ahwh = DhBhwh. Vynásobením této rovnosti maticí Bh zleva a s novým využitím vztahů (4.3) dostaneme λ (Bhwh) = BhDhBhwh = Ah+1 (Bhwh) . To znamená, že vektor Bhwh je vlastním vektorem matice Ah+1 příslušný k dominantní vlastní hodnotě λ. Platí tedy wh+1 = Bhwh Bhwh 1 , h = 0, 1, . . . , m − 1; přitom klademe wm = w0. Známe-li tedy normovaný pravý vlastní vektor matice A0 příslušný k vlastní hodnotě λ, pak můžeme snadno spočítat normované pravé vlastní vektory matic A1, A2, . . . , Am−1. Nechť vh, h = 0, 1, 2, . . . , m − 1 je levý vlastní vektor matice Ah takový, že (vh)1 = 1 (reproduktivní hodnoty vyjadřujeme relativně k reproduktivní hodnotě první třídy). Pak platí vT h+1Ah+1 = vT h+1BhDh = λvT h+1. Vynásobením této rovnosti zprava maticí Bh dostaneme vT h+1Bh DhBh = λ vT h+1Bh , což znamená, že vektor BTvh+1 je levým vlastním vektorem matice DhBh = Ah příslušným k vlastní hodnotě λ. Známe-li tedy levý vlastní vektor vm = v0 matice A0, pak můžeme spočítat levé vlastní vektory matic Am−1, Am−2, . . . , A1 pomocí vztahů vh−1 = BT h−1vh BT h−1vh 1 , h = m, m − 1, . . . , 2. Upozorněme ještě na skutečnost, že společné dominantní vlastní číslo matic Ah (tj. rychlost růstu populace se sezónní variabilitou) v případě, že všechny matice Bh, h = 0, 1, . . . , m − 1 92 KAPITOLA 4. MODELY S EXTERNÍ VARIABILITOU jsou čtvercové (tj. populace je v každé fázi členěna do stejných tříd), nemusí nijak souviset s vlastními hodnotami jednotlivých matic Bh. Například pro m = 2 uvažujme matice B0 = 0,2 0,2 0,9 0 , B1 = 0,1 3 0,2 0 . Dominantní vlastní hodnota matice B0, resp. matice B1, je 0,5359, resp. 0,8262. Z toho by se mohlo zdát, že populace vymírá. Ale dominantní vlastní hodnota matice A0 = 0,1 3 0,2 0 0,2 0,2 0,9 0 = 2,72 0,02 0,04 0,04 je rovna 2,7203, takže populace dosti rychle roste. Tento příklad není nějak umělý, může popisovat populaci, která se v nepříznivém období roku (například období sucha) soustřeďuje na přežívání, v příznivém období na rozmnožování. Uvedený jev tedy může sloužit k analýze strategie dormance semen nebo spor. Na závěr ještě vyšetříme citlivost růstového koeficientu λ na složky matice Bh. Poněvadž podle (4.3) je ∂(Ah)pq ∂(Bh)ij = ∂ ∂(Bh)ij kh+1 l=1 (Dh)pl(Bh)lq = δqj(Dh)pi = δqj(DT h )ip, platí podle řetězového pravidla pro derivování složené funkce ∂λ ∂(Bh)ij = kh p=1 kh q=1 ∂λ ∂(Ah)pq ∂(Ah)pq ∂(Bh)ij = kh p=1 kh q=1 ∂λ ∂(Ah)pq δqj(DT h )ip = = kh p=1 ∂λ ∂(Ah)pj (DT h )ip = DT h ∂λ ∂Ah ij . Označíme-li S(Ah) = ∂λ ∂(Ah)ij = (vh)i(wh)j vT h wh matici citlivosti růstového koeficientu λ na složkách matice Ah (sr. 2.5.1) a S(Bh) = ∂λ ∂(Bh)ij matici citlivosti růstového koeficientu λ na složkách matice Bh, můžeme psát S(Bh) = DT h S(Ah). Matici pružnosti E(Bh) růstového koeficientu λ vzhledem ke složkám matice Bh můžeme zapsat ve tvaru E(Bh) = (Bh)ij λ ∂λ ∂(Bh)ij = 1 λ Bh ◦ S(Bh). 4.2. PERIODICKÁ VARIABILITA 93 4.2 Periodická variabilita Představme si populaci, jež je strukturovaná do k tříd a vyvíjí se v prostředí, které se periodicky mění. To může například být způsobeno měnícím se počasím v průběhu roku a podobně. Takovou populaci můžeme modelovat rovnicí n(t + 1) = A(t)n(t), (4.4) kde o časově závislé projekční matici A(t) předpokládáme, že je pro všechna t = 0, 1, 2, . . . nezáporná a má periodu m, tj. A(t + m) = A(t) a m je kladné celé číslo. Změníme časové měřítko tak, aby délka periody byla jednotková, tj. zavedeme novou nezávisle proměnnou s = 1 m t a položíme ν(s) = n(ms). Pak pro h ∈ {0, 1, 2, . . . , m − 1} a s nezáporné celé číslo platí ν s + h + 1 m = n(ms + h + 1) = A(ms + h)n(ms + h) = = A(h)n m s + h m = A(h)ν s + h m . Model tedy můžeme zapsat ve tvaru ν s + h + 1 m = A(h)ν s + h m , h = 0, 1, . . . , m − 1, s = 0, 1, 2, . . . , což je model (4.1) s τh = h m , Bh = A(h). Model (4.4) můžeme považovat za speciální případ modelu se sezónní externí variabilitou. Alternativu k uvedenému přístupu k modelům s externí periodickou variabilitou představuje využití Fourierovy analýzy. Prvky aij(t) matice A(t) v modelu (4.4) jsou periodické funkce s periodou m. Můžeme je tedy vyjádřit ve tvaru Fourierovy řady aij(t) = cij + ∞ l=1 bij cos 2π m t + ϕij . O koeficientech budeme předpokládat, že |aij| ≥ |bij|, i, j = 1, 2, . . . , k; (4.5) tento předpoklad zaručí, že matice A(t) je nezáporná pro všechna t. Je-li nerovnost v podmínce (4.5) ostrá, pak matice A(t) je primitivní, resp. ireducibilní, pro všechna t právě tehdy, když A(0) je primitivní, resp. ireducibilní. 94 KAPITOLA 4. MODELY S EXTERNÍ VARIABILITOU 4.3 Aperiodická variabilita Uvažujme model vývoje populace strukturované do k tříd s časově závislou projekční maticí A(t) n(t + 1) = A(t)n(t); (4.6) přitom matice A(t) je pro každé t = 0, 1, 2, . . . nezáporná. Řešením této rovnice s počáteční hodnotou n(0) je n(t) = A(t − 1)A(t − 2) · · · A(1)A(0)n(0). (4.7) Abychom mohli model (4.6) nějak analyzovat, potřebujeme zavést několik dalších pojmů z teorie nezáporných matic. Hilbertova projektivní pseudometrika d je definována pro každou dvojici (x, y) nezáporných vektorů se stejným nosičem, tj. takových, že xi = 0 právě tehdy, když yi = 0. Pseudometrika d je definována vztahem d(x, y) =    ln max xi yi : yi > 0 min xi yi : yi > 0 , x = o, 0, x = o. Tvrzení 17. Pseudometrika d má vlastnosti: 1. d(x, y) ≥ 0 pro všechny přípustné vektory x, y. 2. d(x, y) = d(y, x) pro všechny přípustné vektory x, y. 3. d(x, z) ≤ d(x, y) + d(y, z) pro všechny přípustné vektory x, y. 4. d(x, y) = 0 právě tehdy, když x = ay pro nějaké a > 0. 5. d(x, y) = d(ax, by) pro všechny přípustné vektory x, y a kladná čísla a, b. 6. Pro každou nezápornou čtvercovou matici A a všechny přípustné vektory x, y platí d(Ax, Ay) ≤ d(x, y). Je-li d(x, y) > 0, je tato nerovnost ostrá. Důkaz: 1. Tvrzení plyne přímo z definice zobrazení d. 2. Pro nulové vektory je symetrie zřejmá z definice. Pokud y = o, platí xi yi : yi > 0 = xi yi : xi > 0 neboť vektory x a y mají stejné nosiče. Dále max xi yi : yi > 0 = 1 min yi xi : yi > 0 , min xi yi : yi > 0 = 1 max yi xi : yi > 0 a z toho již plyne tvrzení. 4.3. APERIODICKÁ VARIABILITA 95 3. Je-li některý z vektorů x, y, z nulový, jsou nulové všechny a nerovnost je triviální. V opačném případě d(x, y) + d(y, z) = ln max xi yi : yi > 0 min xi yi : yi > 0 + ln max yi zi : yi > 0 min yi zi : zi > 0 = = ln max xi yi : yi > 0 max yi zi : zi > 0 min xi yi : yi > 0 min yi zi : zi > 0 ≥ ln max xi yi yi zi : yizi > 0 min xi yi yi zi : yizi > 0 = d(x, z). 4. Je-li d(x, y) = 0, pak max xi yi : yi = 0 = min xi yi : yi = 0 , a z toho dále plyne, že všechny hodnoty xi yi jsou stejné, a tedy rovny nějaké konstantě a. Opačná implikace je zřejmá. 5. max axi byi : byi = 0 = a b max xi yi : yi = 0 a podobně pro minimum. 6. Nechť i je libovolný index takový, že (Ay)i = 0. Pak platí (Ax)i (Ay)i = j aijxj p aipyp = j aijyj p aipyp xj yj , a poněvadž j aijyj p aipyp = 1, aijyj p aipyp ≥ 0 pro všechna i, j, je hodnota (Ax)i (Ay)i váženým průměrem hodnot z množiny xj yj : yj = 0 . To znamená, že min xj yj : yj = 0 ≤ (Ax)i (Ay)i ≤ max xj yj : yj = 0 (4.8) pro všechny indexy i takové, že (Ay)i = 0. Odtud dále plyne, že min xj yj : yj = 0 ≤ min (Ax)j (Ay)j : (Ax)j = 0 ≤ ≤ max (Ax)j (Ay)j : (Ax)j = 0 ≤ max xj yj : yj = 0 , což je ekvivalentní s dokazovanou nerovností. Pokud je d(x, y) > 0, pak podle již dokázané vlastnosti 4. je alespoň jedna z nerovností v (4.8) ostrá. 96 KAPITOLA 4. MODELY S EXTERNÍ VARIABILITOU Vlastnosti 1., 2., a 3. jsou axiomy pseudometriky. Vlastnosti 4. a 5. říkají, že pseudometrika d nerozlišuje (stotožňuje) vektory, které se liší pouze velikostí, nikoliv směrem. Z vlastnosti 6. plyne, že pro všechny přípustné vektory x, y a nezáporné matice A platí nerovnost d(Ax, Ay) d(x, y) ≤ 1, neboli, že násobení nezápornou maticí A nezvětšuje Hilbertovu pseudovzdálenost vektorů. To nám dále umožňuje pro nezápornou matici A definovat Birkhoffův kontrakční koeficient τ(A) = sup d(Ax, Ay) d(x, y) : d(x, y) > 0 . Tvrzení 18. Koeficient τ(A) má vlastnosti: 1. 0 ≤ d(Ax, Ay) ≤ τ(A)d(x, y), τ(A) ≤ 1 pro všechny přípustné vektory x, y a všechny nezáporné matice A. 2. Pro nezáporné čtvercové matice A, B platí τ(AB) ≤ τ(A)τ(B). 3. Pro nezápornou nenulovou matici A platí, že τ(A) = 0 právě tehdy, když A = cwvT, kde c je nějaká konstanta a w, v jsou levý a pravý vlastní vektor matice A příslušné k její dominantní vlastní hodnotě. 4. Je-li A > 0, pak τ(A) < 1. Důkaz: 1. Plyne přímo z definice koeficientu τ a z vlastností pseudometriky d. 2. Poněvadž podle tvrzení 17.6. platí d(ABx, ABy) ≤ d(Bx, By), je τ(AB) = sup d(ABx, ABy) d(x, y) : d(x, y) ≤ sup d(Bx, By) d(x, y) : d(x, y) = τ(B). Z toho, že τ(A) ≤ 1 dále dostaneme τ(AB) ≤ τ(B) ≤ τ(B)τ(A). 3. Nechť τ(A) = 0. Pak pro všechny vektory x splňující podmínku d(x, y) > 0 platí d(Ax, Ay) = 0, což vzhledem k vlastnosti 17.4 pseudometriky znamená, že existuje číslo a > 0 takové, že Ax = aAy. Násobení maticí A tedy zobrazuje všechny vektory se stejným nosičem do jednorozměrného prostoru a z toho dále plyne, že hodnost matice A je 1. Všechny sloupce matice A jsou tedy násobkem nějakého nezáporného a nenulového vektoru q, tj. A = qpT. Přitom p = o, neboť A = O. Nechť nyní w je vlastní vektor matice A příslušný k dominantní vlastní hodnotě λ. Pak λw = Aw = qpT w = pT w q, tedy q = λ pTw w. Dále platí λvT = vT A = vT qpT = vT λ pTw wpT = λvTw pTw pT , tedy pT = pTw vTw vT . 4.3. APERIODICKÁ VARIABILITA 97 Odtud dostáváme, že A = λ vTw wvT. Nechť A = cwvT. Pak pro libovolné vektory x, y platí Ax = cwvT x = cvT x w, Ay = cwvT y = cvT y w, takže Ax = cwvT x = cvT x w = cvTx cvTy Ay, což podle tvrzení 17.4 znamená, že d(Ax, Ay) = 0. 4. Např. J. E. Carroll, Birkhoff’s contraction coefficient. Linear Algebra and its Applications 389 (2004) 227-234. Vrátíme se nyní k rovnici (4.6) a jejímu řešení (4.7). Položme Ht = A(t − 1)A(t − 2) · · · A(1)A(0). Řekneme, že posloupnost matic {Ht}∞ t=0 je slabě ergodická, pokud lim t→∞ τ(Ht) = 0. Pro každé dva nezáporné vektory x, y se stejným nosičem podle tvrzení 18.1 platí 0 ≤ d(Htx, Hty) ≤ τ(Ht)d(x, y), takže z věty o třech posloupnostech plyne lim t→∞ d(Htx, Hty) = 0 pro slabě ergodickou posloupnost matic {Ht}. Pokud je tedy posloupnost matic {Ht} slabě ergodická, pak řešení rovnice (4.6) mají pro libovolné počáteční podmínky asymptoticky ekvivalentní směr. Z vlastnosti 18.3 můžeme usoudit, že slabě ergodická posloupnost matic {Ht} je asymptoticky ekvivalentní s posloupností matic λtwtvT t , kde λt je dominantní vlastní hodnota matice Ht a wt, resp. vt, je příslušné levý, resp. pravý, vlastní vektor. Řešení rovnice (4.6) je tedy asymptoticky ekvivalentní s posloupností vektorů λtwtvT t n(0) = λtvT t n(0) wt . Slabě ergodická poslupnost matic je tedy jistým zobecněním pojmu primitivní matice. Přes- něji: Pokud je matice A(t) v rovnici (4.6) konstantní a primitivní, tj. A(t) = A pro všechna t ≥ 0 a existuje t0 ≤ 0 takové, že At0 > 0, pak je posloupnost matic {Ht} = At slabě ergodická. Důkaz: Podle tvrzení 18.1 a 18.2 je 0 ≤ τ(At ) = τ At−[t/t0]t0 A[t/t0]t0 ≤ ≤ τ At−[t/t0]t0 τ A[t/t0]t0 ≤ τ At−[t/t0]t0 τ(At0 )[t/t0] ≤ τ(At0 )[t/t0] . Poněvadž podle 18.4 je τ(At0 ) < 1, je lim t→∞ τ(At0 )[t/t0] = 0, takže také lim t→∞ τ(At) = 0. 98 KAPITOLA 4. MODELY S EXTERNÍ VARIABILITOU Nechť jsou všechny matice A(t), t = 0, 1, 2, . . . , v rovnici (4.6) primitivní a mají stejnou incidenční matici, tj. pro všechna t, s ≥ 0 a všechny dvojice indexů i, j platí, že aij(t) = 0 právě tehdy, když aij(s) = 0. Pokud existuje konstanta K taková, že lim sup t→∞ max {aij(t) : aij(0) > 0} min {aij(t) : aij(0) > 0} < K, (4.9) pak je posloupnost matic {Ht} slabě ergodická. Podmínka (4.9) říká, že variabilita prostředí není taková, že by některý koeficient projekční matice „téměř vymizel . Kapitola 5 Modely s interní variabilitou 5.1 Příklad — populace strukturovaná podle plodnosti Uvažujme opět model (9) popisující vývoj populace, v níž lze jedince rozlišit na juvenilní a dospělé (plodné). Projekční matice populace je tvaru A = σ1(1 − γ) ϕ σ1γ σ2 , kde σ1 . . . pravděpodobnost přežití juvenilních jedinců do dalšího období; σ2 . . . pravděpodobnost přežití plodných jedinců do dalšího období; γ . . . pravděpodobnost, že juvenilní jedinec během období dospěje; ϕ . . . střední počet potomků plodného jedince za jedno období. Charakteristicky polynom matice A je na levé straně rovnosti (2.4), takže dominantní vlastní hodnota (závisející na všech parametrech) je podle (2.6) rovna λ1 = λ1(σ1, σ2, γ, ϕ) = 1 2 σ1(1 − γ) + σ2 + σ1(1 − γ) − σ2 2 + 4σ1γϕ . Odtud je vidět, že λ1(σ1, σ2, γ, ϕ) ≥ 1 právě tehdy, když ϕ ≥ (1 − σ2) 1 − σ1(1 − γ) σ1γ . Zejména λ1(σ1, σ2, 1, ϕ) = 1 2 σ2 + σ2 2 + 4σ1ϕ . a λ1(σ1, σ2, 1, ϕ) ≥ 1 právě tehdy, když ϕσ1 ≥ 1 − σ2. Jinak řečeno, populace s bezprostředním dospíváním nevymírá právě tehdy, když plodnost násobená pravděpodobností přežití juvenilních jedinců není menší než úmrtnost dospělých. Každý z ekologických (demografických) parametrů modelu může záviset na velikosti populace nebo na jejím složení (relativním nebo absolutním zastoupením jednotlivých tříd). Velká populace, tj. velká vnitrodruhová konkurence, může omezovat přežití, rychlost dospívání i 99 100 KAPITOLA 5. MODELY S INTERNÍ VARIABILITOU plodnost: σ1 = Σ1e−s11n1−s12n2 , (5.1) σ2 = Σ2e−s21n1−s22n2 , (5.2) γ = Γe−g1n1−g2n2 , (5.3) ϕ = Φe−f1n1−f2n2 . (5.4) Parametry Σ1, Σ2, Γ, Φ lze interpretovat jako pravděpodobnosti přežití juvenilních jedinců, přežití plodných jedinců, maturace během projekčního intervalu a specifickou plodnost dospělých jedinců (v tomto pořadí) za předpokladu, že se neprojevuje vnitrodruhová konkurence. Parametry sij, gi, fi, i, j = 1, 2 udávají „velikost vlivu vnitrodruhové konkurence na přežití, dobu dospívání a plodnost. Všechny parametry jsou nezáporné; pro pravděpodobnosti Γ, Σ1, Σ2 platí: 0 < γ ≤ 1, tj. juvenilní jedinec dospěje v konečném čase, 0 < Σ1 ≤ 1, tj. v reálné populaci nemohou všichni jedinci zemřít před dosažením plodnosti, 0 ≤ Σ2 < 1, tj. dospělí jedinci nemohou neumírat. Parametry σ1, σ2, γ a ϕ budeme chápat jako funkce vektoru n = (n1, n2)T. Označme λ0 1 = λ1 σ1(0), σ2(0), γ(0), ϕ(0) , λ∞ 1 = lim n →∞ λ1 σ1(n), σ2(n), γ(n), ϕ(n) , pokud poslední limita existuje. Platí: Je-li 0 ≤ λ∞ 1 < 1 < λ0 1, pak 0 < inf { n(t) : t ∈ N0} ≤ sup { n(t) : t ∈ N0} < ∞, tj. populace dlouhodobě přežívá a její velikost je omezená. Jsou-li funkce σ1, σ2, γ a ϕ dány rovnostmi (5.1)–(5.4), pak lim n →∞ σ1(n) = 0 pokud min {s11, s12} > 0, lim n →∞ σ2(n) = 0 pokud min {s21, s22} > 0, lim n →∞ γ(n) = 0 pokud min {g1, g2} > 0, lim n →∞ ϕ(n) = 0 pokud min {f1, f2} > 0, dále lim ϕ→0 λ1(σ1, σ2, γ, ϕ) = σ1(1 − γ), lim γ→0 λ1(σ1, σ2, γ, ϕ) = σ1, lim σ1→0 λ1(σ1, σ2, γ, ϕ) = σ2, lim σ2→0 λ1(σ1, σ2, γ, ϕ) = 1 2 σ1(1 − γ) + σ2 1(1 − γ)2 + 4σ1γϕ a funkce λ1 je spojitou funkcí svých proměnných. Odtud plyne: • pokud plodnost závisí na velikosti populace podle vztahu (5.4) s min{f1, f2} > 0 a ostatní parametry modelu jsou konstantní, pak velikost populace nemůže růst neomezeně (neboť σ1(1 − γ) < 1) — stabilizace populace zmenšením plodnosti při velké populační hustotě; 5.2. KONSTRUKCE MODELŮ 101 • pokud převrácená hodnota doby dospívání závisí na velikosti populace podle vztahu (5.3) s min{g1, g2} > 0 a ostatní parametry jsou konstantní, přičemž přežívání juvenilních jedinců není jisté (σ1 < 1), pak populace nemůže růst neomezeně — stabilizace populace odložením reprodukce při velké populační hustotě; • pokud pravděpodobnost přežití juvenilních jedinců závisí na velikosti populace podle vztahu (5.1) s min{s11, s12} > 0 a ostatní parametry jsou konstantní, pak populace nemůže růst neomezeně — stabilizace populace zvětšením úmrtnosti juvenilních jedinců (nebo infanticidou) při velké populační hustotě; • i když pravděpodobnost přežití dospělých jedinců závisí na velikosti populace podle vztahu (5.2) s min{s21, s22} > 0, může populace růst neomezeně; k tomu například dojde, když plodnost je velká, ϕ > 1 − σ1(1 − γ) σ1γ . Stejné úvahy se stejnými závěry lze provést i v případě, že parametry σ1, σ2, γ a ϕ závisí na velikosti populace jiným způsobem, než podle vztahů (5.1)–(5.4), ale stále mají vlastnost lim n →∞ σ1(n) = 0, lim n →∞ σ2(n) = 0, lim n →∞ γ(n) = 0, lim n →∞ ϕ(n) = 0. Na obrázku 5.1 vlevo je znázorněna dynamika populace, jejíž ekologické (demografické) charakteristiky σ1, σ2, γ, ϕ nezávisejí na populační hustotě. Na obrázcích 5.1 vpravo, 5.2, 5.3 vlevo jsou příklady populace se stejnými hodnotami parametrů Σ1, Σ2, Γ a Φ takových, že právě jeden z ekologických (demografických) parametrů závisí na velikosti (hustotě) populace. U populace na obr. 5.3 vlevo se projevuje vliv vnitrodruhové konkurence na přežití dospělých jedinců (např. vnitrodruhová agresivita rostoucí s populační hustotou); tento vliv však nezajistí regulaci velikosti populace. Vliv vnitrodruhové konkurence na přežití dospělých stabilizuje velikost populace při nižší plodnosti, obr. 5.3 vpravo. Na obrázcích 5.4–5.7 jsou ukázky různých invariantních množin a příslušné dynamiky populace. U populací na obrázcích 5.4, 5.5, 5.7 se jedná o stabilizaci populace omezením plodnosti při velkých populačních hustotách, u populace na obrázku 5.6 se jedná o stabilizace populace odložením reprodukce při vyšších populačních hustotách. 5.2 Konstrukce modelů Obecný model růstu strukturované populace s interní variabilitou je tvaru n(t + 1) = A n(t) n(t). (5.5) Čtvercová matice A = A(n) řádu k je pro každý vektor n ∈ ¯Rk + nezáporná. Pokud lze projekční matici A dekomponovat na součet matice přechodů mezi třídami a matice plodností, A(n) = T(n) + F(n), 102 KAPITOLA 5. MODELY S INTERNÍ VARIABILITOU 0 5 10 15 20 050000150000 time juveniles 0 5 10 15 20 010003000 time adults 0 5 10 15 20 0.51.01.5 time juveniles 0 5 10 15 20 0.000.040.08 time adults Obrázek 5.1: Vývoj populace strukturované podle plodnosti. Použité parametry: Σ1 = 0.5, Σ2 = 0.1, Γ = 0.1, Φ = 50, n1(0) = 1, n2(0) = 0. Vlevo: Ekologické parametry nezávisí na velikosti populace, tj. f1 = f2 = 0, g1 = g2 = 0, s11 = s12 = 0, s21 = s22 = 0 a v důsledku λ1 = 1.8658. Vpravo: Plodnost ϕ závisí na velikosti populace, ostatní parametry nikoliv, tj. f1 = f2 = 1, g1 = g2 = 0, s11 = s12 = 0, s21 = s22 = 0 a v důsledku λ0 1 = 1.8658, λ∞ 1 = 0.45. Stabilizace populace omezením plodnosti. 5.2. KONSTRUKCE MODELŮ 103 0 5 10 15 20 0.61.01.41.8 time juveniles 0 5 10 15 20 0.0000.0100.020 time adults 0 5 10 15 200.20.61.0 time juveniles 0 5 10 15 20 0.0000.0100.020 time adults Obrázek 5.2: Vývoj populace strukturované podle plodnosti. Použité parametry: Σ1 = 0.5, Σ2 = 0.1, Γ = 0.1, Φ = 50, n1(0) = 1, n2(0) = 0. Vlevo: Dospívání γ závisí na velikosti populace, ostatní parametry nikoliv, tj. f1 = f2 = 0, g1 = g2 = 1, s11 = s12 = 0, s21 = s22 = 0 a v důsledku λ0 1 = 1.8658, λ∞ 1 = 0.5. Stabilizace populace odložením reprodukce. Vpravo: Přežití juvenilních σ1 závisí na velikosti populace, ostatní parametry nikoliv, tj. f1 = f2 = 0, g1 = g2 = 0, s11 = s12 = 1, s21 = s22 = 0 a v důsledku λ0 1 = 1.8658, λ∞ 1 = 0.1. Stabilizace populace zvětšením úmrtnosti juvenilních (infanticidou). 104 KAPITOLA 5. MODELY S INTERNÍ VARIABILITOU 0 5 10 15 20 04000080000 time juveniles 0 5 10 15 20 050015002500 time adults 0 5 10 15 200.50.70.9 time juveniles 0 5 10 15 20 0.000.020.04 time adults Obrázek 5.3: Vývoj populace strukturované podle plodnosti. Použité parametry: Σ1 = 0.5, Σ2 = 0.1, Γ = 0.1, n1(0) = 1, n2(0) = 0. Vlevo: Přežití dospělých σ2 závisí na velikosti populace, ostatní parametry nikoliv a fertilita je velká, tj. Φ = 50, f1 = f2 = 0, g1 = g2 = 0, s11 = s12 = 0, s21 = s22 = 1 a v důsledku λ0 1 = 1.8658, λ∞ 1 = 1.8221. Zpomalení růstu zvětšením úmrtnosti dospělých. Vpravo: Přežití dospělých σ2 závisí na velikosti populace, ostatní parametry nikoliv a fertilita je malá, tj. Φ = 10.5, f1 = f2 = 0, g1 = g2 = 0, s11 = s12 = 0, s21 = s22 = 1 a v důsledku λ0 1 = 1.0204, λ∞ 1 = 0.9837. Stabilizace populace zvětšením úmrtnosti dospělých. 5.2. KONSTRUKCE MODELŮ 105 0.0 0.5 1.0 1.5 0.000.020.040.060.08 Attractor juveniles adults 0 20 40 60 80 0.00.51.01.52.0 Population projection time juveniles 0 20 40 60 80 0.000.040.08 time adults Obrázek 5.4: Rovnovážný bod. Použité parametry: Σ1 = 0.5, Σ2 = 0.1, Γ = 0.1, Φ = 50, f1 = f2 = 1, g1 = g2 = 0, s11 = s12 = s21 = s22 = 0, n1(0) = 1, n2(0) = 0 106 KAPITOLA 5. MODELY S INTERNÍ VARIABILITOU 0 2 4 6 8 0.00.10.20.30.4 Attractor juveniles adults 0 20 40 60 80 051015 Population projection time juveniles 0 20 40 60 80 0.00.20.40.60.8 time adults Obrázek 5.5: Cyklus periody 4. Použité parametry: Σ1 = 0.5, Σ2 = 0.1, Γ = 0.1, Φ = 500, f1 = f2 = 1, g1 = g2 = 0, s11 = s12 = s21 = s22 = 0, n1(0) = 1, n2(0) = 0 5.2. KONSTRUKCE MODELŮ 107 0 1 2 3 4 5 6 0.0000.0050.0100.015 Attractor juveniles adults 0 50 100 150 200 02468 Population projection time juveniles 0 50 100 150 200 0.0000.010 time adults Obrázek 5.6: Invariantní smyčka. Použité parametry: Σ1 = 0.5, Σ2 = 0.1, Γ = 0.1, Φ = 300, f1 = f2 = 0, g1 = g2 = 1, s11 = s12 = s21 = s22 = 0, n1(0) = 1, n2(0) = 0 108 KAPITOLA 5. MODELY S INTERNÍ VARIABILITOU 0 5 10 15 0.00.20.40.60.8 Attractor juveniles adults 0 50 100 150 200 01020304050 Population projection time juveniles 0 50 100 150 200 0.01.02.0 time adults Obrázek 5.7: Podivný atraktor. Použité parametry: Σ1 = 0.5, Σ2 = 0.1, Γ = 0.1, Φ = 1800, f1 = f2 = 1, g1 = g2 = 0, s11 = s12 = s21 = s22 = 0, n1(0) = 1, n2(0) = 0 5.2. KONSTRUKCE MODELŮ 109 musí prvky tij = tij(n) a fij = fij(n) matic T a A splňovat nerovnosti fij(n) ≥ 0, tij(n) ≥ 0, i, j = 1, 2, . . . , k, k i=1 tij(n) ≤ 1, j = 1, 2, . . . , k pro všechny nezáporné vektory n. Jako vhodný tvar funkcí tij navrhli Fujiwara a Caswell tij(n) = exp(αi + βi T n) 1 + k p=1 exp(αp + βp T n) , i = 1, 2, . . . , k. (5.6) Parametry αi určují pravděpodobnosti přechodu do i-té třídy nebo setrvání v ní při nulové velikosti populace (při tak malé populaci, že se neprojeví vnitrodruhová konkurence ani kooperace). Vektor βi určuje vliv velikostí jednotlivých tříd populace na přechod do i-té třídy nebo přežívání v ní. Pokud budeme ještě uvažovat (k+1)-ní třídu (uhynulé jedince) a položíme tk+1,j(n) = 1 1 + k p=1 exp(αp + βp T n) , představují funkce tij hustotu mnohorozměrného logistického rozdělení pravděpodobonosti. Snadno ověříme, že k+1 i=1 tij(n) = 1. Často je užitečné uvažovat poněkud specifičtější model, konkrétně takový, že všechny složky matice A závisí na váženém součtu velikostí jednotlivých tříd populace N(t) = k i=1 cini(t) = cT n(t), c ≥ 0; veličina N ve speciálním případě c = 1 vyjadřuje celkovou velikost populace. Ve vyjádření pravděpodobností přechodu rovnostmi (5.6) bude u těchto modelů βi = c pro všechna i = 1, 2, . . . , k. Nechť aij = aij(N) je diferencovatelná funkce a označme gij(N) = Naij(N). Pokud aij(N) > 0 pro nějaké N > 0, řekneme, že vliv N na aij je depensující. Pokud aij(N) ≤ 0 a gij(N) ≥ 0 pro všechna N ≥ 0, řekneme, že vliv N na aij je kompensující; je-li přitom lim N→∞ gij(N) = ∞, mluvíme o nedostatečné kompensaci, je-li lim N→∞ gij(N) = 0, mluvíme o nadměrné kompensaci. Často používané závislosti jsou aij(N) = αij 1 1 + βN , Beverton-Holt, kompensující aij(N) = αije−γN , Ricker, nadměrně kompensující; parametry αij, β, γ jsou kladné. 110 KAPITOLA 5. MODELY S INTERNÍ VARIABILITOU 5.3 Asymptotické vlastnosti Označme n(t; n0) řešení projekční rovnice (5.5) s počáteční podmínkou n(0) = n0. Definice 3. XXX • Nechť x0 ∈ ¯Rk +. Množina Ω(x0) = x ∈ ¯Rk + : (∃ {ti}∞ i=0) lim i→∞ ti = ∞, lim i→∞ n(ti; x0) = x se nazývá ω-limitní množina bodu x0 vzhledem k rovnici (5.5). • Nechť M ⊆ ¯Rk +. Množina Ω(M) = x∈M Ω(x) se nazývá ω-limitní množina množiny M vzhledem k rovnici (5.5). • Nechť S ⊆ ¯Rk + je taková množina, že x ∈ S ⇒ (∀t ∈ N) n(t; x) ∈ S. Pak se S nazývá invariantní množina rovnice (5.5). • Nechť S ⊆ ¯Rk + je taková invariantní množina rovnice (5.5), že (∀Q ⊆ S) (S Q = ∅) ⇒ (∃x ∈ Q)(∃t ∈ N)n(t; x) ∈ Q, tj. každá její vlastní podmnožina není invariantní množinou rovnice (5.5). Pak se množina S nazývá minimální invariantní množina rovnice (5.5). Poznámka 3. Množina S je minimální invariantní množinou rovnice (5.5) právě tehdy, když S = Ω(S). Definice 4 (Typy invariantních množin). Minimální invariantní množina S ⊆ ¯Rk + rovnice (5.5) se nazývá stacionární bod (rovnovážný bod, equilibrium), pokud množina S je jednoprvková; cyklus délky (periody) p, pokud p je celé číslo větší než 1 a množina S má p prvků; invariantní smyčka, pokud množina S je uzavřenou křivkou v prostoru Rk; podivná, pokud není žádného z předchozích typů. Vektor ˆn je stacionární bod rovnice (5.5) právě tehdy, když pro každý čas t je n(t; ˆn) = ˆn, tj. vektor ˆn je řešením rovnice A(ˆn)ˆn = ˆn. (5.7) Odtud je vidět, že nulový vektor je stacionárním bodem rovnice (5.5). Tento stacionární bod nazýváme triviální. Netriviální stacionární body vyjadřují stálou velikost i složení populace — velikosti jednotlivých tříd se v průběhu času nemění; populace je v dynamické rovnováze se svým prostředím. 5.3. ASYMPTOTICKÉ VLASTNOSTI 111 Cyklus délky p je množina S = {x0, x1, . . . , xp−1} taková, že A(xi)xi = xi+1 mod p. Zejména pro p = 2 platí S = {x0, x1} a x1 = A(x0)x0, x0 = A(x1)x1, stacionární bod x0 je tedy nenulovým řešením rovnice A A(x0)x0 A(x0)x0 = x0. Analogicky lze hledat cykly délky větší než 2. Definice 5 (Stabilita stacionárních bodů). Stacionární bod ˆn rovnice (5.5) se nazývá stabilní, pokud (∀ε > 0)(∃δ > 0)(∀t ≥ 1) n0 − ˆn < δ ⇒ n(t; n0) − ˆn < ε; asymptoticky stabilní, pokud (∃δ > 0) n0 − ˆn < δ ⇒ lim t→∞ n(t; n0) = ˆn; globálně asymptoticky stabilní, pokud (∀n0 = o) lim t→∞ n(t; n0) = ˆn, repulsivní, pokud (∃ε > 0) n0 = ˆn ⇒ lim inf t→∞ n(t; n0) − ˆn > ε. Poznamenejme, že v definici asymptotické stability stacionárního bodu nepožadujeme stabilitu (na rozdíl od Persidského definice asymptotické stability v teorii obyčejných diferenciálních rovnic). Nechť ˆn je stacionární bod rovnice (5.5) a n je její řešení. Označme x(t) = n(t) − ˆn odchylku řešení od stacionárního bodu. Pak n(t) = ˆn + x(t), a ˆn + x(t + 1) = n(t + 1) = A ˆn + x(t) ˆn + x(t) . (5.8) Předpokládejme dále, že odchylka x od stacionárního bodu je „malá (tj. norma vektoru x je „malá ) a že složky matice A jsou dvakrát spojitě diferencovatelné. Položme B = A(ˆn) + ∂A ∂n1 (ˆn) , ∂A ∂n2 (ˆn) , . . . , ∂A ∂nk (ˆn) ˆn. (5.9) Z rovnosti (5.8) s využitím Taylorovy věty a rovnosti (5.7) nyní dostaneme ˆni + xi(t + 1) = k j=1 aij ˆn + x(t) ˆnj + xj(t) = = k j=1 aij(ˆn) + k l=1 ∂aij ∂nl (ˆn) xl(t) + O x(t) 2 ˆnj + xj(t) = = k j=1 aij(ˆn) ˆnj + xj(t) + k j=1 k l=1 ∂aij ∂nl (ˆn) xl(t) ˆnj + xj(t) + O x(t) 2 = = A(ˆn)ˆn i + k j=1 aij(ˆn)xj(t) + k l=1 xl(t) k j=1 ∂aij ∂nl (ˆn) ˆnj + O x(t) 2 = = ˆni + k l=1 ail(ˆn)xl(t) + k l=1 xl(t) ∂A ∂nl (ˆn) ˆn i + O x(t) 2 = = ˆni + k l=1 ail(ˆn) + ∂A ∂nl (ˆn) ˆn i xl(t) + O x(t) 2 = ˆni + Bx(t) i + O x(t) 2 112 KAPITOLA 5. MODELY S INTERNÍ VARIABILITOU pro libovolné i = 1, 2, . . . , k. Odtud x(t + 1) = Bx(t) + O x(t) 2 . To znamená, že odchylka x od stacionárního bodu ˆn se přibližně vyvíjí jako řešení lineárního homogenního systému diferenčních rovnic s maticí B (která nemusí být nezáporná), tj. x(t) ≈ y(t), kde y je řešení úlohy y(t + 1) = By(t), y(0) = x(0), pokud x(t) je „malá . Označme λB = max |λ| : (∃v ∈ Rk )Bv = λv . (5.10) Je-li λB < 1, pak lim t→∞ y(t) = o, tedy y zůstává „malý a „malá je i odchylka x od stacionárního stavu; v důsledku toho je lim t→∞ x(t) = o, což pro řešení rovnice (5.5) s počáteční podmínkou blízko stacionárního bodu ˆn znamená, že lim t→∞ n(t) = ˆn. Je-li λB > 1, pak lim t→∞ y(t) = ∞. Provedené úvahy můžeme zformulovat jako větu. Věta 5. Nechť ˆn je stacionární bod rovnice (5.5) a matice A je v okolí bodu ˆn dvakrát spojitě diferencovatelná. Definujme matici B rovností (5.9) a číslo λB rovností (5.10). Pak platí: je-li λB < 1, pak je stacionární bod ˆn asymptoticky stabilní, je-li λB > 1, pak je stacionární bod ˆn repulsivní. V případě, že projekční matice A závisí na váženém součtu složek vektoru n, A = A(N) = A(cTn), jsou její první parciální derivace ve stacionárním bodě rovny ∂A ∂nl (ˆn) = clA ( ˆN), kde ˆN = cT ˆn. To znamená, že matici B můžeme zapsat ve tvaru B = A( ˆN) + c1A ( ˆN), c1A ( ˆN), . . . , ckA ( ˆN) ˆn = A( ˆN) + cT ⊗ A ( ˆN) ˆn. Definice 6 (Klasifikace minimálních invariantních množin). Minimální invariantní množina S ⊆ ¯Rk + rovnice (5.5) se nazývá stabilní, pokud ke každému okolí V množiny S existuje okolí U množiny S takové, že z n0 ∈ U plyne n(t; n0) ∈ V pro všechna t ∈ N; atraktor, pokud existuje okolí U množiny S takové, že z n0 ∈ U plyne lim t→∞ inf { n(t; n0) − x : x ∈ S} = 0; okolí U se nazývá oblast přitažení atraktoru; 5.3. ASYMPTOTICKÉ VLASTNOSTI 113 globální atraktor, pokud množina S je atraktor a celá množina U = ¯Rk + {o} je jeho oblastí přitažení; repelor, pokud existuje okolí U množiny S takové, že ke každému počátečnímu stavu n0 ∈ S existuje čas t0 takový, že {n(t; n0) : t ≥ t0} ∩ U = ∅. 114 KAPITOLA 5. MODELY S INTERNÍ VARIABILITOU Kapitola 6 Modely dvojpohlavní populace 6.1 Populace strukturovaná podle plodnosti Uvažujme populaci, v níž jsou jedinci dvou pohlaví a jedinci každého pohlaví jsou rozlišeni na juvenilní a plodné. Juvenilní jedinci mohou přežívat a dospívat, plodní jedinci produkují gamety a mohou přežívat. Jedná se tedy o „dvojpohlavní analogii populace, jejíž model je sestaven na str. 5. Za začátek životního cyklu takové populace budeme považovat spojení samičí a samčí gamety, kterým vznikne zatím bezpohlavní zygota, ze které se vyvine buď mladá samička nebo sameček. Za jednotku času zvolíme dobu, za niž ze zygoty vznikne jedinec. Označme n1 = n1(t) množství životaschopných zygot v čase t, tj. takových zygot, z nichž vznikne živý jedinec. Dále označme n2 = n2(t) a n3 = n3(t), resp. n4 = n4(t) a n5 = n5(t), množství juvenilních a plodných samic, resp. samců. Analogicky jako v případě nepohlavní populace označíme γf nebo γm pravděpodobnosti, že juvenilní samička nebo sameček dospěje během časového intervalu jednotkové délky, a σ1f , σ1m, σ2f a σ2m označíme pravděpodobnosti, že jednotkový čas přežije juvenilní samička, juvenilní sameček, plodná samice a plodný samec (v tomto pořadí). Primární poměr pohlaví, tj. relativní zastoupení zygot, z nichž se vyvinou samičky, označíme . Nakonec označíme symbolem B počet zygot, které během jednotkového času vyprodukují plodní jedinci. Tento počet bude určitě nějak záviset na aktuálním množství f samic a m samců aktuálně přítomných v populaci, tedy B = B(f, m); v uvažovaném modelu je f = n3, m = n5. Životní cyklus uvažované populace je znázorněn na obr. 6.1. Nejedná se ovšem o graf životního cyklu ve vlastním slova smyslu; do uzlu N1 vede „hrana vycházející ze dvou uzlů N3 a N5. Vývoj modelované populace lze popsat rovnostmi n1(t + 1) = B n3(t), n5(t) , n2(t + 1) = n1(t) + σ1f (1 − γf )n2(t), n4(t + 1) = (1 − )n1(t) + σ1m(1 − γm)n4(t), n3(t + 1) = σ1f γf n2(t) + σ2f n3(t), n5(t + 1) = σ1mγmn4(t) + σ2mn5(t), 115 116 KAPITOLA 6. MODELY DVOJPOHLAVNÍ POPULACE   1   2   3   4   5E E d d d ds        ©   c   T  '  ' ' σ2m σ2fσ1f (1 − γf ) σ1m(1 − γm) σ1f γf σ1mγm 1 − B juvenilní dospělí samci samice zygoty Obrázek 6.1: Schematické znázornění životního cyklu dvojpohlavní populace strukturované na juvenilní (neplodné) a dospělé (plodné) jedince. Uzel N1 representuje zygoty, které jsou produkovány společně dospělými samicemi N3 a samci N5, ale které ještě nemají určité pohlaví. Parametr označuje primární poměr pohlaví, parametry γi pravděpodobnost maturace během projekčního intervalu, parametry σ1i a σ2i pravděpodobnosti přežití juvenilních a dospělých jedinců příslušného pohlaví. nebo v maticovém tvaru       n1 n2 n3 n4 n5       (t + 1) = =       B n3(t), n5(t) 0 0 0 0       +       0 0 0 0 0 σ1f (1 − γf ) 0 0 0 0 σ1f γf σ2f 0 0 1 − 0 0 σ1m(1 − γm) 0 0 0 0 σ1mγm σ2m             n1 n2 n3 n4 n5       (t). (6.1) Opět se nejedná o maticový model v obvyklém smyslu; model je porušen (perturbován) přičítaným vektorem, který závisí na struktuře populace a tato závislost je nelineární. Funkce B = B(f, m) vyjadřující množství „vyprodukovaných zygot bývá nazývána funkce rození (birth function) nebo funkce manželství (mariage function). Tvar této funkce je potřebné specifikovat. Funkce rození B = B(f, m) přiřadí danému množství plodných samic f a samců m množství jimi vyprodukovaných zygot. Měla by mít následující vlastnosti: (i) B : [0, ∞) × [0, ∞) → [0, ∞), tj. množství vyprodukovaných zygot je nezáporné číslo. (ii) B(0, n) = 0 = B(n, 0) pro jakékoliv n ≥ 0, tj. pokud v populaci chybí plodné samice nebo samci, nebudou žádné nové zygoty. (iii) B(αf, αm) = αB(f, m) pro každé α > 0, tj. funkce je homogenní prvního řádu. Tato vlastnost vyjadřuje, že počet „vyprodukovaných zygot se mění ve stejném poměru, 6.1. POPULACE STRUKTUROVANÁ PODLE PLODNOSTI 117 v jakém se změní množství dospělých samic i samců; například, pokud se zdvojnásobí množství plodných samic i samců, pak se zdvojnásobí i počet „vyprodukovaných zygot (pro α = 2). (iv) Funkce B je neklesající v každém svém argumentu, tj. zvětší-li se množství plodných samic (nebo samců), množství „vyprodukovaných zygot se nezmenší. První dvě vlastnosti jsou přirozené, třetí a čtvrtá mohou být předmětem diskuse. Např. zdvojnásobí-li se počet plodných samic i samců, může být výsledné množství vyprodukovaných zygot větší než dvojnásobné (při větší populační hustotě může být větší šance, že se samice se samcem setkají, projevuje se Alleeho efekt), nebo menší než dvojnásobné (při velké populační hustotě mohou být spotřebovány zdroje prostředí a na jedince zbývá méně energie pro produkci gamet, projevuje se vnitrodruhová konkurence). Podobné námitky lze mít i proti neklesání funkce B. Můžeme uvažovat různé strategie oplodňování. Jedna extrémní možnost je tvorba trvalých párů, tj. že jeden samec během jednotkového času oplodní nejvýše jednu samici a jedna samice je oplodněna nejvýše jedním samcem. Navíc těchto párů je maximální možný počet — pokud není samic méně než samců, pak všichni samci realizují své spermie, pokud není samic více než samců, pak všechny samice jsou oplodněné. Druhá krajnost je ta, že všechny dospělé samice jsou oplodněny, pokud je v populaci alespoň jeden plodný samec, nebo že všichni samci realizují své gamety, pokud je v populaci alespoň jedna samice (to je možné například u dvoudomých rostlin). Samozřejmě může nastat také nějaká možnost mezi těmito krajnostmi. Tyto úvahy vedou k závěru, že hodnoty funkce B(f, m) mohou být úměrné některému z následujících výrazů: min {f, m} maximální možné množství párů, f, m > 0 0, m = 0 dominance samic, m, f > 0 0, f = 0 dominance samců, qf + (1 − q)m, fm = 0 0, fm = 0 vážený aritmetický průměr, pro váhu q platí 0 < q < 1, √ fm geometrický průměr, 2fm f + m harmonický průměr. Všechny takové funkce splňují vlastnosti (i)–(iv). Předpokládejme nyní, že známe funkci B. Pro libovolné hodnoty f ≥ 0 a m ≥ 0 položíme F3(f, m) = B(f, m) 2f , F5(f, m) = B(f, m) 2m . 118 KAPITOLA 6. MODELY DVOJPOHLAVNÍ POPULACE Pak B(f, m) = fF3(f, m) + mF5(f, m) a model (6.1) můžeme přepsat jako maticový model tvaru       n1 n2 n3 n4 n5       (t + 1) = =       0 0 F3 n3(t), n5(t) 0 F5 n3(t), n5(t) σ1f (1 − γf ) 0 0 0 0 σ1f γf σ2f 0 0 1 − 0 0 σ1m(1 − γm) 0 0 0 0 σ1mγm σ2m             n1 n2 n3 n4 n5       (t). 6.2 Věkově strukturovaná dvojpohlavní populace Uvažujme populaci tvořenou jedinci dvou pohlaví, kteří jsou charakterizováni svým věkem a kteří tvoří poměrně stabilní páry. Páry mohou vznikat nebo se rozpadat, páry mohou plodit potomky, jedinci přežívají nebo umírají. Za časovou jednotku (délku projekčního intervalu) zvolíme takový čas, během kterého může u každého jedince dojít k nejvýše jedné z událostí: vytvoření páru s jedincem opačného pohlaví, rozpad páru v němž byl zapojen, úmrtí. Nechť k označuje věk vyjádřený v této časové jednotce, který nemůže samec ani samice překročit. Označme fi = fi(t), resp. mj = mj(t), množství nespárovaných samic věkové třídy i, resp. nespárovaných samců věkové třídy j, v čase t. Dále označme cij = cij(t) množství párů, v nichž samice je z věkové třídy i a samec z věkové třídy j; takové páry budeme stručně nazývat „páry typu (i, j) . Při uvedeném označení je celkový počet samic z věkové třídy i, resp. samců z věkové třídy j, roven fi + k j=1 cij, resp. mj + k i=1 cij. Střední množství potomků, které vyprodukuje pár typu (i, j) během projekčního intervalu označíme Fij. Dále budeme předpokládat, že mezi novorozenci je konstantní podíl samic a že novorozenci nejsou spárovaní. Z těchto předpokladů plyne, že množství novorozených samic a samců a množství párů tvořených novorozenci je dáno rovnostmi f1(t + 1) = k i,j=1 Fijcij(t), m1(t + 1) = (1 − ) k i,j=1 Fijcij(t), c11(t) = 0 (6.2) pro každé t = 0, 1, 2, . . . . Označme dále P (f) i , resp. P (m) j , pravděpodobnost, že samice věkové třídy i, resp. samec věkové třídy j, přežije projekční interval, Dij pravděpodobnost, že se pár typu (i, j) během projekčního intervalu rozpadne a oba partneři přežijí (tuto pravděpodobnost můžeme nazvat rozvodovost, divorse rate), Mij množství párů typu (i, j), které se během projekčního intervalu vytvoří. Budeme předpokládat, že pravděpodobnosti přežití i rozvodovost nezávisí na 6.2. VĚKOVĚ STRUKTUROVANÁ DVOJPOHLAVNÍ POPULACE 119 velikosti ani struktuře populace. Naopak množství nově vzniklých párů závisí přinejmenším na množství a věkovém složení nespárovaných samic a samců, tj. Mij = Mij(f, m). Tuto funkci nazýváme funkce partnerství, mating function. Množství nespárovaných samic věkové třídy i + 1 v čase t + 1 je tvořeno těmi, které měly v čase t věk ze třídy i a přežily projekční interval zmenšené o ty, které během projekčního intervalu vytvořily pár s nějakým samcem. K nim přibudou samice, které byly spárovány s nějakým samcem a tento pár se během projekčního intervalu rozpadl nebo jim partner uhynul. Tedy fi+1(t + 1) = P (f) i fi(t) − k j=1 Mij f(t), m(t) + k j=1 Dij + (1 − P (m) j )P (f) i cij(t) (6.3) pro všechna t = 0, 1, 2, . . . , i = 1, 2, . . . , k. Podobně pro množství nespárovaných samců platí mj+1(t + 1) = P (m) j mj(t) − k i=1 Mij f(t), m(t) + k i=1 Dij + (1 − P (f) i )P (m) j cij(t) (6.4) pro všechna t = 0, 1, 2, . . . , j = 1, 2, . . . , k. Z párů typu (i, j), v nichž oba partneři přežijí a které se nerozpadnou po uplynutí projekčního intervalu budou páry typu (i+1, j+1). K nim se přidají páry vzniklé během projekčního intervalu ze samice věkové třídy i a samce věkové třídy j. Tedy ci+1,j+1(t + 1) = P (f) i P (m) j − Dij cij(t) + Mij f(t), m(t) (6.5) pro všechna t = 0, 1, 2, . . . , i, j = 1, 2, . . . , k. Model vývoje uvažované populace je zapsán rovnicemi (6.2), (6.3), (6.4), (6.5). Přitom primární poměr pohlaví , párově specifické plodnosti Fij, věkově specifické koeficienty přežívání P (f) i a Pm j a rozvodovost Dij splňují nerovnosti 0 < < 1, Fij ≥ 0, 0 ≤ P (f) i ≤ 1, 0 ≤ P (m) j ≤ 1, 0 ≤ Dij ≤ P (f) i P (m) j pro všechna i, j = 1, 2, . . . , k. Podmíněná pravděpodobnost, že pár typu (i, j) se během projekčního intervalu rozpadne za podmínky, že oba partneři přežijí, je dána výrazem dij = Dij P (f) i P (m) j , pokud P (f) i P (m) j > 0; jinak položíme dij = 0. Pro podmíněnou rozvodovost dij platí 0 ≤ dij ≤ 1 pro všechna i, j = 1, 2, . . . , k. Rovnice (6.3), (6.4), (6.5) můžeme při tomto označení přepsat na tvar fi+1(t + 1) = P (f) i  fi(t) + k j=1 1 − P (m) j (1 − dij) cij(t)   − k j=1 Mij f(t), m(t) , mj+1(t + 1) = P (m) j mj(t) + k i=1 1 − P (f) i (1 − dij) cij(t) − k i=1 Mij f(t), m(t) , ci+1,j+1 = P (f) i P (m) j (1 − dij)cij(t) + Mij f(t), m(t) . 120 KAPITOLA 6. MODELY DVOJPOHLAVNÍ POPULACE 6.2.1 Funkce partnerství Funkce partnerství M = M(f, m) musí být nezáporná, tj. M : [0, ∞)2k → [0, ∞). Dále by měla mít následující vlastnosti: P1) Pokud v populaci není nespárovaná samice věkové třídy i nebo samec věkové třídy j, pak pár typu (i, j) nevznikne. Tedy Mij(f1, . . . , fi−1, 0, fi+1, . . . , fk, m1, . . . , mk) = 0, Mij(f1, . . . , fk, m1, . . . , mj−1, 0, mj+1, . . . , mk) = 0 pro libovolné nezáporné vektory f, m. P2) Celkový počet nově spárovaných samic věkové třídy i nemůže být větší, než byl celkový počet nespárovaných samic této věkové třídy a podobně pro samce. Tedy k j=1 Mij(f, m) ≤ fi, k i=1 Mij(f, m) ≤ mj pro libovolné nezáporné vektory f, m a všechny indexy i, j = 1, 2, . . . , k. P3) Funkce M je homogenní řádu p > 0, tedy Mij(αf, αm) = αp Mij(f, m) pro libovolné nezáporné vektory f, m, kladné číslo α a pro všechny indexy i, j = 1, 2, . . . , k. Pokud se v populaci projevuje vnitrodruhová konkurence, je p > 1, pokud se v populaci projevuje Alleho efekt, je p < 1; sr. diskusi k vlastnosti (iii) u populace strukturované podle plodnosti v 6.1. P4) Pokud se zvětší počet nespárovaných samic věkové třídy i a samců věkové třídy j, nezmenší se počet nově vznikajících párů typu (i, j). Tedy pro všechny nezáporné vektory f, m, ˜f, ˜m takové, že f ≤ ˜f a m ≤ ˜m a pro všechny indexy i, j = 1, 2, . . . , k platí Mij(f, m) ≤ Mij(f1, . . . , fi−1, ˜fi, fi+1, . . . , fk, m1, . . . , mj−1, ˜mj, mj+1, . . . , mk), P5) Na „manželském trhu je konkurence. Pokud počet nespárovaných samic věkové třídy i a samců věkové třídy j se nezmění, ale přibudou nějací nespárovaní jedinci jiných věkových tříd, může se zmenšit počet nově vznikajících párů typu (i, j); nespárovaná samice věkové třídy i může najít partnera v jiné věkové třídě než j a podobně pro samce. Tedy pro všechny nezáporné vektory f, m, ˜f, ˜m takové, že f ≤ ˜f a m ≤ ˜m a pro všechny indexy i, j = 1, 2, . . . , k platí Mij(f, m) ≥ Mij( ˜f1, . . . , ˜fi−1, fi, ˜fi+1, . . . , ˜fk, ˜m1, . . . , ˜mj−1, mj, ˜mj+1, . . . , ˜mk). 6.2. VĚKOVĚ STRUKTUROVANÁ DVOJPOHLAVNÍ POPULACE 121 Nejjednodušší funkce partnerství je taková, že množství vzniklých párů typu (i, j) závisí pouze na množství nespárovaných samic věkové třídy i a samců věkové třídy j. V takovém případě ovšem v podmínkách P4) a P5) budou rovnosti. Dostatečně obecná funkce tohoto typu je Hadelerova funkce Mij(f, m) =    pij qfr i + (1 − q)mr j 1/r , fimj > 0, 0, fimj = 0, kde r ∈ R, q ∈ [0, 1] a pij jsou nezáporná čísla taková, aby byla splněna podmínka P2). Nechť fimj > 0. Pro q = 1, resp. q = 0, dostaneme Mij(f, m) = pijfi, resp. Mij(f, m) = pijmj; jedná se tedy o dominanci samic (polygynii), resp. dominanci samců (polyandrii). Nechť nyní q ∈ (0, 1). Pro r = 1 dostaneme Mij(f, m) = pij qfi + (1 − q)mj vážený aritmetický průměr, pro r = −1 dostaneme Mij(f, m) = pij fimj qmj + (1 − q)fi vážený harmonický průměr a pro r → 0 dostaneme Mij(f, m) = pijfq i m1−q j vážený geometrický průměr; průměry jsou nevážené (nevychýlené pro některé pohlaví), pokud q = 1 2. Nakonec pro r → −∞, resp. r → ∞, dostaneme Mij(f, m) = pij min{fi, mj}, resp. Mij(f, m) = pij max{fi, mj}. Realističtější funkce partnerství, která závisí na množství nespárovaných samic a samců všech věkových tříd, může být tvaru Mij(f, m) = pij fimj k i=1 fi + k j=1 mj .