Západočeská univerzita v Plzni Fakulta aplikovaných věd

Z´apadoˇcesk´a univerzita v Plzni Fakulta aplikovan´ych vˇed

Anal´ yza populaˇ cn´ıch model˚ u semestr´aln´ı pr´ace KMA/MM

Radim Hoˇsek, A07237 [email protected]

Kapitola

1

Populaˇcn´ı modely Mezi v´ yznamn´ y okruh studovan´ ych a modelovan´ ych dˇej˚ u patˇr´ı tzv. dynamick´e procesy, tedy procesy vyv´ıjej´ıc´ı se v z´avislosti na ˇcasov´e promˇenn´e. Jedn´ım z takov´ ychto (i pro laika nejsrozumitelnˇejˇs´ıch) dynamick´ ych dˇej˚ u je ’ v´ yvoj urˇcit´e populace, at uˇz lidsk´e, nebo obecnˇe libovoln´e biologick´e, coˇz bude i t´ematem n´ asleduj´ıc´ıch str´anek. Populac´ı rozum´ıme skupinu ˇziv´ ych organism˚ u jednoho biologick´eho druhu. V´ıce populac´ı r˚ uzn´ ych druh˚ u pak naz´ yv´ame biologick´e spoleˇcenstv´ı. Prvn´ı model v´ yvoje poˇctu jedinc˚ u v populaci vytvoˇril anglick´ y ekonom Thomas Robert Malthus (1766–1834) a tento model dnes naz´ yv´ame Malthus˚ uv populaˇcn´ı z´ akon. Malthus˚ uv model, jak za okamˇzik uk´aˇzeme, pˇredpokl´ ad´ a exponenci´ aln´ı r˚ ust populace, kter´ y nen´ı (aˇz na nˇekter´e speci´aln´ı a sp´ıˇse umˇel´e pˇr´ıpady) realistick´ y. Nen´ı proto divu, ˇze jiˇz od 19. stolet´ı byla snaha populaˇcn´ı modely zpˇresnit, pˇr´ıkladem m˚ uˇze b´ yt Verhulst˚ uv model z roku 1845.

1.1

Principy populaˇ cn´ıch model˚ u

At’ uˇz je podoba v´ ysledn´ ych model˚ u jak´akoli, to z´asadn´ı z˚ ustav´a nemˇenn´e: Proces: V´ yvoj (r˚ ust) mnoˇzstv´ı biologick´ ych jedinc˚ u jednoho druhu v nˇejak´e oblasti. Z´ akonitost: Populaˇcn´ı z´akon r˚ ustu pro jednu populaci s jist´ ym typem rozmnoˇzov´ an´ı.

2

´ 1.2. OD ZAKONITOSTI K MODELU

Radim Hoˇsek

Obecn´ e matematick´ e vyj´ adˇ ren´ı: Okamˇzit´ y pˇr´ır˚ ustek populace z´avis´ı na mnoˇzstv´ı jedinc˚ u (resp. p´ar˚ u schopn´ ych rozmnoˇzov´an´ı) a na migraci. dN (t) = aN (t) + M (t), dt

t ∈ h0, T i,

(1.1)

kde a > 0 je populaˇcn´ı koeficient, hustota populace N (t) oznaˇcuje pr˚ umˇernou stˇredn´ı koncentraci populace v uvaˇzovan´em prostoru v ˇcasov´em okamˇziku t, M (t) je migraˇcn´ı funkce. Obecn´e bilanˇcn´ı pˇredpoklady o zmˇen´ach hustoty populace se naz´ yvaj´ı z´ akony populaˇcn´ı dynamiky.

1.2

Od z´ akonitosti k modelu

Pˇredpokl´ adejme nyn´ı i prostorov´e rozloˇzen´ı populace, pak oznaˇcujeme N (x, t) poˇcet jedinc˚ u populace v jednotce objemu v okamˇziku t (resp. ”kolem bodu x v okamˇziku t”). Integr´ al pˇres bilanˇcn´ı oblast Z

N (x, t) dx

N (ΩB , t) = N (t) = ΩB

oznaˇcuje poˇcet jedinc˚ u v oblasti v okamˇziku t. Budeme pˇredpokl´ adat standartn´ı sloˇzen´ı populace z hlediska rozmnoˇzov´an´ı (tj. zpravidla rovnomˇern´e rozloˇzen´ı mezi obˇe pohlav´ı) a naopak nebudeme do u ´vah zahrnovat zpoˇzd’ovac´ı efekty (tˇehotenstv´ı, larv´aln´ı stadia, ˇcas od narozen´ı k sexu´ aln´ı vyspˇelosti,. . . ). Jakkoli zn´ı tyto restrikce zvl´aˇstnˇe, napˇr´ıklad u modelov´ an´ı populac´ı jednobunˇeˇcn´ ych organism˚ u jsou pomˇernˇe opodstanˇen´e. Podobnˇe jako u model˚ u fyzik´aln´ıch jev˚ u zmˇena hustoty populace v kaˇzd´em (libovoln´em) ˇcasov´em u ´seku (t1 , t2 ) je rovna celkov´emu pˇr´ır˚ ustku (´ ubytku) hustoty populace na tomto u ´seku, coˇz lze vyj´adˇrit integr´alem Zt2

[M (t) + rN (t)] dt, t1

v nˇemˇz

Rt2

M (t) dt je celkov´ a migrace v ht1 , t2 i (M je tzv. migraˇcn´ı hustota) a

t1

Rt2

rN (t) dt je zmˇena vyvolan´a zrozen´ım a u ´mrt´ım. Koeficient r z´avis´ı obecnˇe

t1

na N , tj. r = r(N (t)) a urˇcuje tzv.relativn´ı vitalitu populace. Glob´aln´ı tvar z´ akona pak m˚ uˇzeme ps´ at ve tvaru 3

1.3. MALTHUS˚ UV MODEL

Radim Hoˇsek

N (t2 ) − N (t1 ) =

Zt2

[M (t) + rN (t)] dt,

∀ht1 , t2 i ⊂ h0, T i.

(1.2)

t1

Je-li M + rN spojit´ a funkce, pak limitn´ım pˇrechodem pro t1 → t, t2 → t, z´ısk´ ame lok´ aln´ı tvar populaˇcn´ıho z´akona dN (t) = r(N (t))N (t) + M (t), t ∈ (0, T ). (1.3) dt V dalˇs´ıch modelech budeme migraci pro jednoduchost zanedb´avat.

1.3

Malthus˚ uv model

V Malthusovˇe modelu pˇredpokl´ad´ame konstantn´ı koeficient r(N (t)) = α, lok´ aln´ı tvar Malthusova populaˇcn´ıho z´akona je proto: dN (t) = αN (t), t ∈ (0, T ). (1.4) dt Pˇredpokl´ adejme, ˇze zn´ ame poˇcet jedinc˚ u populace v ˇcase t = 0: N (0) = N0 ≥ 0. Pak ˇreˇs´ıme rovnici(1.4) jako poˇc´ateˇcn´ı u ´lohu element´arn´ı metodou separace promˇenn´ ych: Z Z dN = αdt, N po integraci: lnN = αt + c, a po u ´pravˇe a dosazen´ı poˇca´teˇcn´ı podm´ınky pak: N (t) = N0 eαt . Snadno ovˇeˇr´ıme, ˇze plat´ı: lim N (t) = +∞

t→+∞

pro α > 0

a lim N (t) = 0

t→+∞

pro α < 0.

Volba α = 0 je z biologick´eho hlediska nesmysln´a a ˇreˇsen´ım je konstantn´ı stav populace N (t) ≡ N0 . D´ıky jednoduchosti Malthusova modelu jsme mohli urˇcit analytick´e ˇreˇsen´ı lok´ aln´ıho tvaru z´ akona. Ne vˇzdy to ale bude moˇzn´e, pro zkoum´an´ı vlastnost´ı sloˇzitˇejˇs´ıch model˚ u budeme muset pouˇz´ıt jin´ y apar´at, kter´ y pˇredstav´ıme pr´ avˇe na tomto jednoduch´em modelu. 4

1.3. MALTHUS˚ UV MODEL

Radim Hoˇsek

ˇ sen´ı Malthusova modelu pro α = 1 a r˚ Obr´ azek 1.1: Reˇ uzn´e hodnoty N0 .

1.3.1

Anal´ yza Malthusova modelu

T´ımto apar´ atem jsou prostˇredky matematick´e anal´ yzy pro vyˇsetˇrov´an´ı stability stacion´ arn´ıch bod˚ u obyˇcejn´ ych diferenci´aln´ıch rovnic. Pˇekn´ y souhrn naleznete napˇr´ıklad v [3]. Mˇejme rovnici tvaru dx(t) = f (x(t)). (1.5) dt Postup je n´ asleduj´ıc´ı: Nejprve urˇc´ıme stacion´arn´ı body rovnice, tzn. takov´e hodnoty xi , pro kter´e je derivace nulov´a, tzn. z (1.5) plat´ı, ˇze f (x) = 0. df (xi ) Z hodnoty f 0 (xi ) = pak urˇc´ıme stabilitu stacion´arn´ıho bodu. Je-li dx f 0 (xi ) kladn´e, je xi nestabiln´ım rovnov´aˇzn´ ym bodem, naopak, je-li z´aporn´e, je asymptoticky stabiln´ı. Pˇr´ıpad, kdy f 0 (xi ) = 0 nelze snadno rozhodnout a vyˇzaduje dalˇs´ı anal´ yzu. V pˇr´ıpadˇe Malthusova modelu m´ame stacion´arn´ı bod pouze jedin´ y, a to N = 0. A protoˇze derivace prav´e strany podle N je rovna α pro vˇsechna N a pˇredpokl´ ad´ ame α > 0, je N = 0 nestabiln´ım stacion´arn´ım bodem, t´eˇz negativn´ım atraktorem. Prakticky to znamen´a, ˇze pro nulov´ y poˇc´ateˇcn´ı stav populace z˚ ust´ av´ a velikost populace nemˇenn´a a pro libovolnˇe mal´ y (kladn´ y) 5

1.4. VERHULST˚ UV MODEL

Radim Hoˇsek

poˇc´ ateˇcn´ı stav populace bude stav populace r˚ ust. Viz obr´azek (1.1). Naopak, je-li α < 0, je N = 0 asymptoticky stabiln´ım stacion´arn´ım bodem, neboli positivn´ım atraktorem. To je koneckonc˚ u v souladu s analytick´ ym ˇreˇsen´ım rovnice modelu.

1.4

Verhulst˚ uv model

Malthus˚ uv model m´ a jednu z´asadn´ı nev´ yhodu - umoˇzn ˇuje neomezen´ y r˚ ust populace, coˇz je ˇcasto v rozporu s naˇs´ı pˇredstavou. Asi se shodneme, ˇze r˚ ust libovoln´e pozemsk´e populace je omezen pˇrinejmenˇs´ım velikost´ı naˇs´ı planety. Proto v roce 1825 do Malthusova modelu zavedl urˇcit´ y jednoduch´ y regulaˇcn´ı mechanismus belgick´ y matematik Pierre Fran¸cois Verhulst. Ten zvolil jakou specifickou m´ıru r˚ ustu line´ arn´ı funkci r(N (t)) = aN + b, a > 0, ˇcastˇeji se ovˇsem pouˇz´ıv´ a tvar r(N (t)) = α(1 − N/K), α > 0, pro jasnˇejˇs´ı interpretaci konstant. Vznikl tak Verhulst˚ uv model, jindy t´eˇz zvan´ y logistick´ y. Lok´aln´ı podoba populaˇcn´ıho z´ akona podle Verhulsta je tedy: dN (t) N (t) = αN (t)(1 − ), t ∈ (0, T ). (1.6) dt K Snadno nahl´edneme, ˇze pro K → +∞ pˇrech´az´ı Verhulst˚ uv model v Malthus˚ uv. Kladnou konstantu K lze interpretovat jako kapacitu prostˇred´ı, α je pak konstanta urˇcuj´ıc´ı tzv. reprodukˇcn´ı r˚ ust. I rovnici popisuj´ıc´ı chov´an´ı Verhulstova modelu lze ˇreˇsit analyticky, pomoc´ı separace promˇenn´ ych: Z

dN N (1 −

Z N K)

=

α dt,

a substituc´ı u = N/K spolu s rozkladem na parci´aln´ı zlomky z´ısk´ame: ln

N = αt + c. |K − N |

Je-li opˇet N0 = N (0) zn´ am´e, dostaneme v´ ysledn´e ˇreˇsen´ı ve tvaru N (t) =

1.4.1

KN0 eαt , K − N0 + N0 eαt

∀t ∈ (0, T ).

(1.7)

Anal´ yza Verhulstova modelu

Opˇet provedeme anal´ yzu modelu pomoc´ı vlastnost´ı stacion´arn´ıch bod˚ u. Z prav´e strany rovnice (1.6) z´ısk´ ame

6

1.5. GOMPERTZ˚ UV MODEL

Radim Hoˇsek



f (N ) = αN 1 −

N . K

Nulov´ ymi body t´eto funkce jsou N1 = 0 a N2 = K.

Obr´ azek 1.2: Z´avislost prav´e strany Verhulstova z´akona na promˇenn´e N. Konstanty maj´ı hodnoty α = 3, K = 4.

Spoˇc´ıt´ ame derivaci f podle N: 

f 0 (N ) = α 1 − 2

N . K

Dosad´ıme-li za N hodnotu N1 a N2 , zjist´ıme, ˇze f 0 (N1 ) = α > 0 ⇒ N1 = 0 je nestabiln´ı stacion´arn´ı bod, zat´ımco f (N2 ) = −α ⇒ N2 = K je asymptoticky stabiln´ı stacion´arn´ı bod. D´ıky t´eto anal´ yze si vˇsimneme dalˇs´ı vlastnosti stacion´arn´ıch bod˚ u: Je-li funkce f(N) spojit´ a a m´ a-li rovnice f(N) = 0 pouze jednoduch´e koˇreny, pak se pravidelnˇe stˇr´ıd´ a jejich stabilita. Toto tvrzen´ı lze dok´ azat prostˇredky elemenent´arn´ı matematick´e anal´ yzy, staˇc´ı si uvˇedomit, ˇze stabilitu zde pˇredstavuje znam´enko prvn´ı derivace.

1.5

Gompertz˚ uv model

Dalˇs´ım modelem je model vytvoˇren´ y anglick´ ym matematick´ ym samoukem Benjaminem Gompertzem (1779–1865). Ten uvaˇzoval relativn´ı vitalitu jako funkci logaritmickou r(N (t)) = −α ln (N/K), α > 0. Kladn´a konstanta K opˇet pˇredstavuje kapacitu prostˇred´ı. Model je tedy reprezentov´an vztahem 7

1.5. GOMPERTZ˚ UV MODEL

Radim Hoˇsek

N (t) dN (t) = −αN (t) ln , t ∈ (0, T ). (1.8) dt K Zde uˇz je analytick´e ˇreˇsen´ı ponˇekud sloˇzitˇejˇs´ı a proto se omez´ıme na rozbor stacion´ arn´ıch bod˚ u. Opˇet m´a rovnice dva stacion´arn´ı body, jsou to N1 = 0 a N2 = K. Opˇet je i zde nula negativn´ım a K positivn´ım atraktorem, nebot’ derivace prav´e strany rovnice (1.8) podle N : N K



f 0 (N ) = −α K + ln

je pro N1 = 0 kladn´ a a pro N2 = K z´aporn´a. Jeˇstˇe dod´ame, ˇze ˇreˇsen´ım (1.8) je tzv. Gompertzova kˇrivka s pˇredpisem αt )

N (t) = N (0)eK(1−e

.

(1.9)

Ve srovn´ an´ı s logistickou kˇrivkou jako ˇreˇsen´ı Verhulstova modelu roste Gompertzova kˇrivka pro niˇzˇs´ı N rychleji. (Viz obr´azek (1.3).) Gompertz˚ uv model dobˇre modeluje napˇr´ıklad chov´an´ı rakovinn´ ych bunˇek, kter´e se pˇri n´ızk´em poˇctu zaˇcnou br´ anit sv´emu z´aniku zv´ yˇsenou reprodukc´ı.

Obr´ azek 1.3: Srovn´an´ı v´ysledk˚ u Gompertzova a Verhulstova modelu.

8

Radim Hoˇsek

1.6

1.6. HARVESTING

Harvesting

Harvesting je anglick´e oznaˇcen´ı pro sklizeˇ n. Modely harvestingu jsou tedy populaˇcn´ı modely, kde je vnˇejˇs´ım z´asahem zmenˇsov´ana velikost populace, at’ uˇz je to sklizn´ı nebo tˇreba predac´ı. Sklizeˇ n m˚ uˇze m´ıt nejr˚ uznˇejˇs´ı strategie, pro jednoduchost se zamˇeˇr´ıme na dvˇe nejjednoduˇsˇs´ı - konstantn´ı a line´arn´ı.

1.6.1

Harvesting - konstantn´ı sklizeˇ n

Pro u ´ˇcely harvestingu n´ am poslouˇz´ı Verhulst˚ uv model obohacen´ y o funkci g(N ) ≡ h, tedy:  dN (t) N (t)  = αN (t) 1 − − h, (1.10) dt K kde h je konstatn´ı rychlost skl´ızen´ı. Nen´ı asi ˇz´adn´ ym pˇrekvapen´ım, ˇze nelze skl´ızet v´ıce, neˇz kolik m´ame. Pokud tedy zn´azorn´ıme nulov´e body prav´e strany rovnice (1.10) jako pr˚ useˇc´ıky dvou funkc´ı, z jejichˇz rozd´ılu se skl´ad´a, bude konstanta vˇzdy prot´ınat parabolu v hodnot´ach

N1,2

K = α± 2

s

α2 − 4

h K

a nav´ıc omez´ıme platnost naˇseho z´akona na interval N ∈ (N1 , K).

Obr´ azek 1.4: Hled´an´ı stacion´arn´ıch bod˚ u modelu harvestingu s konstantn´ı sklizn´ı. Uvˇedomme si, ˇze hodnota h neovlivn´ı hodnotu derivace prav´e strany, bilance dvou stacion´ arn´ıch bod˚ u tedy z˚ ust´av´a stejn´a jako u logistick´eho modelu. Prvn´ı je nestabiln´ı, druh´ y je stabiln´ı. Kdybychom neuˇcinili ono omezen´ı, skl´ızeli bychom, i kdyby nebylo co skl´ızet. Avˇsak i tento model 9

Radim Hoˇsek

1.6. HARVESTING

m˚ uˇze m´ıt sv´e opodstatnˇen´ı — napˇr´ıklad ve finanˇcn´ım sektor jako model neomezen´eho kontokorentn´ıho u ´ˇctu.

1.6.2

Harvesting - progresivn´ı sklizeˇ n

Pr´ avˇe nepˇr´ıjemnost s dod´ av´an´ım doplˇ nuj´ıc´ıch podm´ınek do pˇredchoz´ıho modelu n´ as m˚ uˇze v´est k jin´e strategii skliznˇe. Tou je napˇr´ıklad jej´ı line´arn´ı z´ avislost na velikosti populace. Koeficient line´arn´ı z´avislosti pak znaˇc´ıme e (z anglick´eho effort = u ´sil´ı). Lok´aln´ı tvar populaˇcn´ıho z´akona s progresivn´ı sklizn´ı je tedy  dN (t) N (t)  = αN (t) 1 − − eN (t), t ∈ (0, T ), (1.11) dt K v´ yznam konstant je zn´ am´ y z pˇredchoz´ıch pˇr´ıpad˚ u: α > 0 je koeficient reprodukˇcn´ıho r˚ ustu, K > 0 je kapacita prostˇred´ı, e > 0 je u ´sil´ı, se kter´ ym je skl´ızeno. Vyˇsetˇr´ıme opˇet stacion´arn´ı body, tedy pˇr´ıpady, kdy je prav´a strana vztahu (1.11) nulov´ a. Odtud

N = eN. (1.12) K Jedn´ım ˇreˇsen´ım je N1 = 0. Pokud existuje druh´ y stacion´arn´ı bod, mus´ı platit e < α, nebot’ hodnota v´ yrazu v z´avorce v (1.12) je menˇs´ı neˇz 1. Zjistili jsme, ˇze hodnota e = α je jakousi mezn´ı hodnotou pro charakter chov´an´ı modelu, vyˇsetˇrov´ an´ı stac.bod˚ u tedy rozdˇel´ıme na dva pˇr´ıpady. 

αN 1 −

1. e < α Jak jsme uk´ azali, za tohoto pˇredpokladu bude m´ıt (1.12) dvˇe ˇreˇsen´ı: N1 = 0 a N2 = K(1 − αe ). Derivac´ı prav´e strany rovnice (1.12) podle N z´ısk´ ame pˇredpis N − e, (1.13) K a tedy f 0 (0) = α − e > 0, nula je nestabiln´ım bodem, zat´ımco N2 je stabiln´ı, nebot’ f 0 (N2 ) = e − α < 0, coˇz lze z´ıskat pouh´ ym dosazen´ım pˇredpisu pro N2 do (1.13). 

f 0 (N ) = α 1 − 2

2. e > α Rovnice (1.12) m´ a jedin´e ˇreˇsen´ı, a to N1 = 0, opˇet plat´ı f 0 (0) = α − e, coˇz je ovˇsem narozd´ıl od pˇredchoz´ıho pˇr´ıpadu z´aporn´e a nula je tedy stabiln´ım stacion´ arn´ım bodem. 10

Radim Hoˇsek

1.6. HARVESTING

V bodˇe e = α tedy doch´ az´ı ke kvalitativn´ı zmˇenˇe v prostoru ˇreˇsen´ı rovnice (1.11), tedy bifurkaci. Odpov´ıd´a to i ”selsk´emu rozumu” - nadmˇern´e zvyˇsov´ an´ı skliznˇe/lovu m´ a za n´asledek vyhuben´ı populace. Viz napˇr´ıklad regulace lovu velryb, coˇz nen´ı nic jin´eho neˇz horn´ı omezen´ı konstanty e v tomto modelu. Bifurkace nen´ı pˇr´ıliˇs pevnˇe vymezen´ y pojem a proto se lze setkat s r˚ uznou podobou bifurkaˇcn´ıch sch´emat a diagram˚ u. V naˇsem pˇr´ıpadˇe jej lze zn´azornit napˇr´ıklad v rovinˇe eN . Pˇreruˇsovan´a ˇc´ara znamen´a nestabiln´ı stacion´arn´ı bod, zat´ımco pln´ a stabiln´ı. Viz obr´azek (1.5).

Obr´ azek 1.5: Bifurkaˇcn´ı diagram - progresivn´ı harvesting: z´avislost a stabilita stacion´ arn´ıch ˇreˇsen´ı na parametru e. Parametry: α = 2, K = 4.

1.6.3

Harvesting pro modelov´ an´ı predace

Modely harvestingu se ovˇsem nemus´ı t´ ykat nutnˇe jen skliznˇe. Obecnˇe ˇreˇceno, jde o vnˇejˇs´ı z´ asah do rozvoje populace, kter´ ym m˚ uˇze b´ yt tˇreba regulace populace nˇejak´ ym pred´ atorem. Jsou zde ovˇsem jin´e poˇzadavky na funkci h. Lze pˇredpokl´ adat, ˇze pˇri velk´em rozsahu populace jiˇz pred´ator nelov´ı, nebot’ se nasyt´ı, a pˇri velmi mal´e populaci lov´ı tak´e m´alo - je pro nˇej koˇrist z d˚ uvodu obt´ıˇzn´eho nalezen´ı nezaj´ımav´a. Takov´eto podm´ınky splˇ nuje funkce BN 2 , (1.14) A2 + N 2 kde A a B jsou kladn´e konstanty. Pro anal´ yzu tohoto modelu pouˇzijeme substituˇcn´ı vztahy h(N ) =

u=

N K αA Bt ,κ = ,% = ,τ = . A A B A 11

Radim Hoˇsek

1.6. HARVESTING

Rovnice  dN (t) N (t)  BN (t)2 , = αN (t) 1 − − 2 dt K A + N (t)2

∀t ∈ (0, T )

(1.15)

B  T . A

(1.16)

pak pˇrech´ az´ı do tvaru  u(τ )2 du(τ ) u(τ )  − = %u(τ ) 1 − , dτ K 1 + u(τ )2



∀τ ∈ 0,

Nyn´ı model opˇet podrob´ıme n´am jiˇz dobˇre zn´am´e anal´ yze. Hled´ame stacion´ arn´ı body, tedy takov´ a u, pro kter´a je prav´a strana vztahu (1.16) nulov´a. Snadno kaˇzd´ y nahl´edne, ˇze takov´ ym bodem je u =0 a protoˇ e ze derivace prav´ u 2u strany rovnice (1.16) podle u je rovna v´ yrazu % 1 − 2 K − (1+u , kter´ y 2 )2 je pro u = 0 roven % > 0, je nula negativn´ım atraktorem. Existuj´ı-li dalˇs´ı stacion´ arn´ı body, mus´ı b´ yt stacion´arn´ımi body rovnice (1.16) i po vydˇelen´ı prav´e strany faktorem u, tedy: 

0=% 1−

u u − . K 1 + u2

(1.17)

Obr´ azek 1.6: Nenulov´e stacion´arn´ı body rovnice (1.16) v z´avislosti na parametrech % a κ.

Naˇs´ım c´ılem je tedy naj´ıt ˇreˇsen´ı t´eto rovnice. Zn´azorn´ıme graficky (viz obr´ azek (1.6)). Rovnice (1.17) m´a bud’ jeden nebo tˇri jednoduch´e koˇreny a nebo dva koˇreny, z toho jeden dvojn´asobn´ y. Pr´avˇe ty budeme vyˇsetˇrovat jako pˇrechodov´ y stav mezi prvn´ımi dvˇema variantami. Je-li u ˜ dvojn´asobn´ ym koˇrenem rovnice (1.17), pak mus´ı b´ yt koˇrenem i jeho rovnice 12

Radim Hoˇsek

1.6. HARVESTING

% 1 − u2 , (1.18) − K 1 + u2 kter´ a vznikla derivov´ an´ım (1.17) podle u. Rovnice (1.17) a (1.18) n´am tak tvoˇr´ı soustavu rovnic promˇenn´ ych u, % a κ, jejichˇz postupnou u ´pravou z´ısk´ ame vztah 0=−

0 = u3 − κu2 + κ.

(1.19)

Obr´ azek 1.7: Grafick´e zn´azornˇen´ı rovnice (1.19). Je to polynom tˇret´ıho ˇra´du, kter´ y m´a v z´avislosti na κ bud’ jeden nebo tˇri re´ aln´e koˇreny. To je v souladu s grafickou pˇredstavou - jak je vidˇet na obr´ azku (1.7), jeden koˇren bude pro κ > 0 (jin´e neuvaˇzujeme) z´aporn´ y, a ten nem´ a v naˇsem modelu v´ yznam. Tuto ˇc´ ast jiˇz nelze ˇreˇsit analyticky, v literatuˇre obvykle najdeme na tomto m´ıstˇe u ´koly pro ˇcten´ aˇre typu ”ovˇeˇrte, ˇze hraniˇcn´ı kˇrivkou (v z´avislosti na κ) je...”. Vyˇreˇs´ıme ji tedy pomoc´ı numerick´ ych ˇreˇsiˇc˚ u SW Matlab. V´ ysledkem je bifurkaˇcn´ı diagram v rovinˇe κ% (1.8), kter´ y lze interpretovat takto: Svˇetl´ a oblast odpov´ıd´ a takov´ ym dvojic´ım koeficient˚ u κ, %, pro nˇeˇz m´a rovnice (1.18) 3 stacion´ arn´ı body, zat´ımco tmav´a tˇem dvojic´ım, kde je stacion´ arn´ı bod pouze jeden. U tohoto modelu tak´e doch´az´ı k zaj´ımav´emu jevu, kter´ y lze pozorovat pˇri n´ asleduj´ıc´ı manipulaci s parametry: Mˇejme m´ alo poˇcetnou populaci. Zafixujeme-li parametr κ a budeme mˇenit pouze parametr % (populace zvyˇsuje svou rozmnoˇzovac´ı aktivitu). V mal´ ych hodnot´ ach m´ a model (kromˇe st´ale nestabiln´ı nuly) dalˇs´ı stacion´arn´ı bod (u1 ). Ten mus´ı b´ yt podle tvrzen´ı vyˇrˇcen´eho v minul´e kapitole stabiln´ı a bude se u nˇej ustalovat i naˇse populace. Pokud budeme pro dan´e κ d´ale 13

Radim Hoˇsek

1.6. HARVESTING

Obr´ azek 1.8: Bifurkaˇcn´ı diagram modelu predace. zvyˇsovat hodnotu %, pˇribydou v modelu dalˇs´ı dva stacion´arn´ı body - u2 (nestabiln´ı) a u3 stabiln´ı. Ten je ovˇsem pro naˇsi populaci pˇres nestabiln´ı bod u2 nedosaˇziteln´ y. Teprve kdyˇz hodnota % pˇrekroˇc´ı urˇcitou mez, kdy zmiz´ı 1. a 2. stacion´arn´ı bod, pˇriklon´ı se stav populace ke stabiln´ımu bodu u3 . Kdyby nyn´ı zaˇcala populace sniˇzovat svou rozmnoˇzovac´ı aktivitu, a t´ım i koeficient %, bude se stav populace ustalovat kolem stabiln´ıho stacion´arn´ıho bodu u3 , a to i po vzniku dalˇs´ıch dvou stac. bod˚ u u1 a u2 . N´aˇs model tedy nen´ı reversibiln´ı (nelze se vr´ atit do stejn´eho stavu pˇri stejn´em nastaven´ı parametr˚ u) - projevuje se jev zvan´ y hystereze.

14

Kapitola

2

Lotk˚ uv-Volterr˚ uv model v´ıce populac´ı V pˇredchoz´ı ˇc´ asti jsme doˇsli aˇz k modelu, kde doch´azelo k ovlivnˇen´ı jedn´e populace populac´ı jin´eho druhu — zde se jednalo o vztah dravec-koˇrist. Pˇresto jsme modelovali pouze onu prvn´ı populaci. Nyn´ı alespoˇ n kr´atce nahl´edneme do modelov´ an´ı v´ yvoje obou interaguj´ıc´ıch populac´ı. Nejjednoduˇsˇs´ım modelem, kter´ y vysvˇetl´ıme pr´ avˇe na pˇr´ıkladu koˇristi a dravce, je Lotk˚ uv-Volterr˚ uv model, kter´ y byl publikov´ an v polovinˇe dvac´at´ ych let 20. stolet´ı nez´avisle na sobˇe americk´ ym matematikem a chemikem Alfredem Lotkou (1880–1949) a italsk´ ym matematikem Vito Volterrou (1860–1940). Je vystavˇen na nˇekolika pil´ıˇr´ıch: • Kaˇzd´ a z obou populac´ı by se bez pˇr´ıtomnosti druh´e chovala podle Malthusova modelu (tj. exponenci´aln´ı r˚ ust nebo pokles). • Populaˇcn´ı koeficient je sniˇzov´an (zvyˇsov´an) vlivem vz´ajemn´e interakce populac´ı. • Nen´ı uvaˇzov´ ana ˇz´ adn´ a saturace, tj. jeden dravec m˚ uˇze seˇzrat libovoln´e dostupn´e mnoˇzstv´ı koˇristi — interakce tak z´avis´ı jen na pravdˇepodobnosti setk´ an´ı. Ze zm´ınˇen´ ych pˇredpoklad˚ u pak plynou n´asleduj´ıc´ı vztahy: M˙ (t) = M (t)(a − bN (t)), N˙ (t) = N (t)(−c + dM (t)), ∀t ∈ (0, T ) 15

´ ´ ´ICH ROVNIC Radim Hoˇsek 2.1. ANALYZA SOUSTAV DIFERENCIALN

pro situaci, kde M je (koncentrace) populace koˇristi, a N (koncentrace) populace dravce, a, b, c, d jsou kladn´e konstanty.

2.1

Anal´ yza soustav diferenci´ aln´ıch rovnic

V pˇr´ıpadˇe jedn´e diferenci´ aln´ı rovnice bylo vyˇsetˇrov´an´ı stability stacion´arn´ıch bod˚ u pomˇernˇe jednoduchou z´aleˇzitost´ı. Ve dvou dimenz´ıch, tedy u soustavy dvou diferenci´ aln´ıch rovnic se postupuje obdobnˇe. Cel´a problematika je samozˇrejmˇe sloˇzitˇejˇs´ı, pokus´ım se zde pouze o kr´atk´ y (a ne´ upln´ y) n´ahled. Mˇejme soustavu x(t) ˙ = f (x(t), y(t)), y(t) ˙ = g(x(t), y(t)).

(2.1)

Nejprve hled´ ame stacion´arn´ı body, tj. dvojice (˜ x, y˜), pro kter´e plat´ı f (˜ x, y˜) = 0, g(˜ x, y˜) = 0. Analogi´ı derivace ve stacion´arn´ım bodˇe je zde Jacobiho matice prvn´ıch parci´ aln´ıch derivac´ı: J(x, y) =

∂f ∂x (x, y) ∂g ∂x (x, y)

∂f ∂y (x, y) ∂g ∂y (x, y)

!

Z Jacobiho matice pro dan´ y stac.bod pak spoˇc´ıt´ame vlastn´ı ˇc´ısla (λ1 , λ2 ) a jim pˇr´ısluˇsn´e vlastn´ı vektory (v1 , v2 ). Charakter stacion´arn´ıho bodu je urˇcen znam´enkem vlastn´ıch ˇc´ısel. • λ1 < 0, λ2 < 0: stabiln´ı uzel (positivn´ı atraktor) — pˇritahuje v obou smˇerech v1 , v2 , • λ1 < 0, λ2 > 0 (nebo naopak): nestabiln´ı stac. bod (sedlo) — pˇritahuje ve smˇeru v1 , odpuzuje ve smˇeru v2 (nebo naopak), • λ1 > 0, λ2 > 0: nestabiln´ı uzel (negativn´ı atraktor) — odpuzuje v obou smˇerech v1 , v2 , • λ1 , λ2 ∈ / R, Reλ1 > 0 : nestabiln´ı ohnisko (negativn´ı atraktor) — odpuzuje ve spir´ale, • λ1 , λ2 ∈ / R, Reλ1 < 0 : stabiln´ı ohnisko (positivn´ı atraktor) — pˇritahuje ve spir´ale, 16

ˇ 2.2. MODEL DRAVEC-KORIST

Radim Hoˇsek

• λ1 , λ2 ∈ / R, Reλ1 = 0 : bod rotace - ˇreˇsen´ı ob´ıhaj´ı periodicky okolo stac.bodu. Vˇse ozˇrejm´ıme na n´ asleduj´ıc´ıch pˇr´ıkladech.

2.2

Model dravec-koˇ rist

Jak jsme jiˇz uvedli, Lotk˚ uv-Volterr˚ uv model interakce dvou populac´ı v rol´ıch dravce a koˇristi vypad´ a n´ asledovnˇe: M˙ (t) = M (t)(a − bN (t)), N˙ (t) = N (t)(−c + dM (t)),

(2.2)

∀t ∈ (0, T ) : M, N ≥ 0, a, b, c, d > 0. Provedeme nyn´ı anal´ yzu modelu. Pro sn´ıˇzen´ı poˇctu parametr˚ u pˇrevedeme rovnici pomoc´ı substituce x=

d M, c

y=

b N a

(2.3)

do tvaru

x(t) ˙ = a(1 − y(t))x(t), y(t) ˙ = −c(1 − x(t))y(t).

(2.4)

Stacion´ arn´ı body t´eto rovnice jsou [0, 0] a [1, 1]. Spoˇcteme Jacobiho matici 

J(x, y) = 

a(1 − y) cy 

−ax −c(1 − x)



.

a 0 Protoˇze J(0, 0) = je diagon´aln´ı, jsou vlastn´ı ˇc´ısla λ1 = a 0 −c a λ2 = −c, pˇr´ısluˇsn´e vlastn´ı vektory pak v1 = (1, 0) a v2 = (0, 1). Bod [0, 0] je tedy sedlem, kter´e pˇritahuje populaci N populaci M .  a odpuzuje  0 −a Zat´ımco bod [1, 1] s Jacobiho matic´ı J(1, 1) = m´a vlastn´ı ˇc´ısla c 0 √ ryze imagin´ arn´ı λ1,2 = ±i ac, a je tedy bodem rotace. Znamen´a to, ˇze stavy populac´ı se periodicky opakuj´ı, coˇz potvrzuje i v´ ystup z numerick´eho experimentu na obr´ azku (2.1). 17

Radim Hoˇsek

ˇ ´I DVOU POPULAC´I 2.3. MODEL SOUPEREN

Obr´ azek 2.1: Grafick´e zn´azornˇen´ı periodick´ych ˇreˇsen´ı Lotkova-Volterrova modelu dravec-koˇrist pro r˚ uzn´e poˇc´ ateˇcn´ı podm´ınky.

2.3

Model soupeˇ ren´ı dvou populac´ı

Ze stejn´ ych pˇredpoklad˚ u jako u modelu dravec-koˇrist se vych´az´ı i pˇri tvorbˇe Lotkova-Volterrova modelu soupeˇren´ı dvou populac´ı. M˙ (t) = M (t)(a − bN (t)), N˙ (t) = N (t)(c − dM (t)).

(2.5)

Opˇet plat´ı, ˇze uvaˇzujeme pouze nez´aporn´e hodnoty M, N a a, b, c, d jsou kladn´e konstanty. Za pouˇzit´ı stejn´e substituce jako v pˇredchoz´ım pˇr´ıpadˇe z´ısk´ ame dvouparametrickou soustavu rovnic ve tvaru x(t) ˙ = a(1 − y(t))x(t), y(t) ˙ = c(1 − x(t))y(t).

(2.6)

Stacion´ arn´ı body jsou stejn´e jako u pˇredchoz´ıho modelu — [0, 0] a [1, 1], jejich charakter bude ovˇsem odliˇsn´ y. Jacobiho matice je ve tvaru 

J(x, y) =

a(1 − y) −cy

−ax c(1 − x)

a tedy: • (x, y) = (0, 0) 

J(0, 0) = 18

a 0 0 c



,



,

Radim Hoˇsek

ˇ ´I DVOU POPULAC´I 2.3. MODEL SOUPEREN

λ1 = a > 0, jemu pˇr´ısluˇsn´ y vlastn´ı vektor je v1 = (1, 0), λ2 = c > 0, jemu odpov´ıd´a vlastn´ı vektor v2 = (0, 1) — nestabiln´ı uzel, • (x, y) = (1, 1) 

J(1, 1) =

0 −a −c 0



,

√ λ1,2 = ± ac, — sedlo, kladn´emu vlastn´ımu ˇc´ıslu (nestabiln´ı smˇer) √ odpov´ıd´ a vlastn´ı vektor v1 = ( ac, −c), z´aporn´emu (stabiln´ı smˇer) √ pak v2 = ( ac, c). Model nem´ a stabiln´ı stacion´arn´ı bod - vede k vyhuben´ı jedn´e a exponenci´ aln´ımu r˚ ustu druh´e populace – viz obr´azek (2.2).

Obr´ azek 2.2: Grafick´e zn´azornˇen´ı ˇreˇsen´ı Lotkova-Volterrova modelu soutˇeˇzen´ı. R˚ uzn´e poˇc´ ateˇcn´ı podm´ınky.

2.3.1

Aplikace modelu soupeˇ ren´ı jako modelu bitvy

Pˇ r´ıklad V bitvˇe u Iwo Jima (1945) v Tich´em oce´ anu proti sobˇe st´ aly americk´ a a japonsk´ a arm´ ada. Uvaˇzujeme pouze stˇrelbu na c´ıl, kdy poˇcet zabit´ych voj´ ak˚ u z´ avis´ı pouze na poˇctu u ´toˇcn´ık˚ uau ´spˇeˇsnosti stˇrelby. Ty bereme konstatn´ı s hodnotou αj = 0, 0544, αa = 0, 0106. Poˇc´ ateˇcn´ı stavy vojsk jsou A(0) = 66454, J(0) = 18274. Naˇs´ım c´ılem je na z´ akladˇe tˇechto informac´ı rozhodnout, jak dopadla bitva. Pˇredpoklady odpov´ıdaj´ı tˇem, za kter´ ych byl sestaven Lotk˚ uv-Volterr˚ uv model soupeˇr´ıc´ıch populac´ı, pouze koeficienty rozmnoˇzov´an´ı z˚ ustanou nulov´e. Lze tedy odvodit lok´ aln´ı vztah, platn´ y pro t ∈ (0, T ): 19

´ 2.4. MODEL SYMBIOZY

Radim Hoˇsek

˙ A(t) = −αj J(t), ˙ = −αa A(t). J(t)

(2.7)

Pokud nav´ıc bilanˇcn´ı rovnosti obohat´ıme o zn´am´e poˇc´ateˇcn´ı podm´ınky, totiˇz A(0) a J(0), lze soustavu rovnic snadno numericky ˇreˇsit. Anal´ yzou jedin´eho stacion´ arn´ıho bodu [0, 0] pak m˚ uˇzeme zjist´ıme, ˇze se jedn´a o sedlo se stabiln´ım smˇerem m´ıˇr´ıc´ım do prvn´ıho (a tˇret´ıho) kvadrantu. To odpov´ıd´a pˇredstavˇe vyhuben´ı jedn´e populace (vojsko). O kterou populaci (=arm´adu) p˚ ujde, z´ aleˇz´ı nejen na parametrech, ale i na poˇc´ateˇcn´ıch podm´ınk´ach. Jak vidno z grafu (2.3), bitvu u Iwo Jima prohr´ali Japonci.

Obr´ azek 2.3: Grafick´e zn´azornˇen´ı ˇreˇsen´ı modelu bitvy u Iwo Jima.

2.4

Model symbi´ ozy

Vˇse zavrˇs´ıme modelem symbi´ozy. Plat´ı pˇredpoklady pro Lotk˚ uv-Volterr˚ uv model, tedy exponenci´ aln´ı r˚ ust jednotliv´ ych populac´ı, kdyby nebyla interakce, a vliv interakce na mnoˇzen´ı populac´ı bez saturace. Za tˇechto okolnost´ı lze podobnˇe jako v pˇredchoz´ıch pˇr´ıpadech odvodit model, jehoˇz lok´aln´ı vyj´ adˇren´ı je M˙ (t) = M (t)(a + bN (t)), N˙ (t) = N (t)(c + dM (t)), ∀t ∈ (0, T ) : M, N ≥ 0, a, b, c, d > 0.

20

(2.8)

´ 2.4. MODEL SYMBIOZY

Radim Hoˇsek

Postupujeme jako obvykle — za pouˇzit´ı substituce (2.3) z´ısk´ame soustavu

x(t) ˙ = a(1 + y(t))x(t), y(t) ˙ = c(1 + x(t))y(t),

(2.9)

se stacion´ arn´ımi body [0, 0] a [−1, −1]. Protoˇze ovˇsem uvaˇzujeme pouze kladn´e hodnoty M, N i vˇsech konstant, nen´ı stacion´arn´ı bod [−1, −1] v naˇsem modelu dosaˇziteln´ y a je pro n´as tud´ıˇz nezaj´ımav´ y. Vyˇsetˇr´ıme proto pouze bod [0, 0]. Jacobiho matice je ve tvaru 

J(x, y) = 

a(1 + y) cy



ax c(1 + x)



,

a 0 , λ1 = a > 0, λ2 = c > 0 — poˇc´atek je tedy 0 c nestabiln´ı uzel, coˇz je koneckonc˚ u ve shodˇe s intuitivn´ı pˇredstavou. Maj´ıli populace exponenci´ aln´ı r˚ ust jiˇz bez interakce, s positivn´ı interakc´ı bude jejich r˚ ust tak´e neomezen´ y. Anal´ yza obyˇcejn´ ych diferenci´aln´ıch rovnic a jejich soustav, jimiˇz jsou ˇcasto modely reprezentov´ any, pomoc´ı vyˇsetˇrov´an´ı stacion´arn´ıch bod˚ u nen´ı rozhodnˇe samosp´ asn´ y n´ astroj. Je to vˇsak n´astroj pomˇernˇe jednoduch´ y a poskytuje vˇsak cenn´ y n´ ahled do chov´an´ı pˇr´ısluˇsn´ ych model˚ u, a proto je ˇcasto velmi uˇziteˇcn´e jej vyuˇz´ıvat. a tedy J(0, 0) =

21

Literatura

[1] M´ıka S.: Element´ arn´ı dynamick´e procesy. text k pˇredn´aˇsk´am KMA/MM, ˇ 2006. Plzeˇ n, ZCU [2] Kalas J., Posp´ıˇsil Z.: Spojit´e modely v biologii. Brno, Masarykova univerzita 2001. [3] Istas J.: Mathematical Modelling for the Life Sciences. Grenoble, Universit´e Pierre Mendes, 2005. [4] Murray J.D.: Mathematical Biology. Second edition, Springer 1993. [5] M´ıka S.: Modely ˇziv´e pˇr´ırody. Pozn´amky z pˇredn´aˇsky J. Miloty.

22