Diskrétní modely Diskrétní modely jednodruhových populacíjednodruhových populací

Mathematical Biology - Murray

1989

Jednoduché modely Kontinuální modely - diferenciální rovnice

Diskréní modely - diferenční rovnice

- populace roste v diskrétních krocích (např. dělení buňky)

- jen přibližné řešení

- různě dlouhé kroky mezi jednotlivými generacemi

- časový krok 1

- vztah mezi populací v čase t+1 (Nt+1) a t (Nt)

Diferenční rovnice)()(1 tttt NfNFNN

Hledáme obecně nelineární funkci , popisující systém. Vycházíme z vypozorovaných faktů.

)( tNf

známe funkci Nt rekurzivně vypočítáme Nt+1

Funkce Nt 0

počáteční stav N0

Popis funkce vyjádřené diferenční rovnicí

• Na začátku pro malé hodnoty Nt – nic nebrání k množení růst

relativně největší

• Pro velké Nt – přemnožení nedostatek potravy vymírání.

Samoregulace růstů - funkce má maximum

Graf závislosti následujícího stavu na aktuálním.

Nejjednodušší popis populačního modeluPopulace v následujícím kroku bude násobek aktuálního stavu. Platí vztah

tt NrN .1 0NrN tt 0)( rNF t

Populace roste geometrickou řadou.

• Nevýhoda: model není skutečný pro většinu populací ani pro dlouhé časy.

• Využití: pro počáteční stádia růstu bakterií.

Musíme také uvažovat, že ne všechny bakterie přežijí. Zobecníme vzorec na

r > 1 …zrůst r < 1 …vymírání

St NrN .1 b

tS NN 1tS NN

NS …všechny bakterie co přežijí a můžou se množit,

b = konst…vyjadřuje ty bakterie co nepřežijí

Dva modely popisují praktické situace v biologii

1. model - Verhulstův proces – ale pro velká Nt nabývá záporných

hodnot

, r > 0, K > 0

K

NrNN t

tt 11

2. model - lze použít pro velké Nt.

, r > 0, K > 0

Exponenciála vyjadřuje faktor úmrtnosti v závislosti na čase.

• pro malé hodnoty Nt je význam exponenciály malý (lineární charakter)

• pro velké hodnoty Nt je její význam velký

K

Nr

tt

t

eNN1

1

K

Nr t

e

Cobwebbing „vytváření pavučiny“:Grafický postup řešení informací o dynamickém chování populace

Rovnovážné stavy popisuje forma

nebo

Rov. stavy jako průsečíky křivky a přímky

)()( **** NFNNfN 0* N 1)( * NF

)(1 tt NfN tt NN 1

*NNm Nt spěje monotónně

k N* pro *

0 NN *0 NN

) (0 1N f N

t

Iterace = převod funkční hodnoty zpět na argument

*10 ... NNN

Vlastní hodnota rovnovážného stavu N*

)(

)( ** Nf

NNdN

Ndf

tt

t … důležitý parametr stability pro

lokální chování okolo rovnovážné polohy pro malé perturbace1)(0 * NfRovnovážné stavy N* je stabilní pro

nestabilní pro

Pro

1)( * Nf

mNN * 0)( * Nf

a) stabilní rovnovážný stav N* pro

b) neutrální stav N* pro

0)(1 * Nf

klesající oscilace

1)( * Nf

periodické oscilace

c) nestabilní stav N* pro

rostoucí oscilace

1)( * Nf

Shrnutí:

rovnováha je stabilní pro

nestabilní pro

kritická bifurkace

11

),1()1,(

1 11

tečná bifurkace vidlová bifurkace

Globální chování Je-li rovnovážný stav nestabilní , grafické řešení nezkonverguje do rovnovážné polohy – LOKÁLNĚ NESTABILNÍ

Je-li populace je ohraničená hodnotou Nmax a Nmin, ať startujeme

z kterýchkoli počátečních podmínek – GLOBÁLNĚ STABILNÍ

Pak řešení náhodně cestuje kolem nedefinovatelné hodnoty chaotické chování

1

Diskrétní logistický model: Chaos Uvazujme nelineární model , r > 0, K > 0

Určíme si , r > 0

Zajímá nás řešení ut 0, to odpovídá intervalu 0 < u0 < 1

Stabilní stav a odpovídající vlastní

hodnoty z

K

NrNN t

tt 11

K

Nu tt )1(1 ttt uruu

0* u rf )0(

r

ru

1* ruf 2)( *

)1( *** uruu

Rovnovážné stavy a bifurkaceZvyšujeme parametr r

0 < r < 1 … rovnováha jen pro u* = 0 0 < < 1 … stabilní

r = 1 … rovnováha jen pro u* = 0 = 1 … nestabilní (1. bifurkace)

1 < r < 3 … pro u* = 0 1 < < 3 … nestabilní

… pro –1 < < 1 … stabilní

r = 3 … pro = – 1 … nestabilní (2. bifurkace)

r

ru

1*

r

ru

1*

Iterativní procesJe zřejmé, že platí:

1. iterace:

2. iterace:

Úpravou a dosazením za ut+2 = ut = u2* dostaneme

Stabilní stavy:

když r > 1 … jeden stabilní stav

když r > 3 … dva stabilní stavy

)( 01 ufu )())(( 0

202 ufuffu

)( 0ufu tt

)1(1 ttt uruu

)1(1)1()(22 tttttt uruururufu

0)1()1()1( *2

*2

2*2

*2 rurrurrruu

0*2 u 0

1*2

r

ru

02

)3)(1()1(*2

r

rrru

A, B, C … rovnovážné body u2*

B > 1 … nestabilní stav pro 2. iteraci

–1 < A < 1, –1 < C < 1 … stabilní stavy

Při každé bifurkaci se předchozí stav stává nestabilním (přerušovaná čára) a stabilní řešení se rozdělí na dvě větve.Úsek stabilního řešení má periodu 2, 22, 23,…

Pro r4 < r < r8

4-cyklé periodické řešení

Kritické body

)( *2*2 AA ufu

***1 )( AAA uufu

**1

**1 , ACCA uuuu

Zvyšováním r:

Přibývá tak pravděpodobných řešení v intervalu (0,1). Vzdálenosti mezi bifurkačními větvemi se zmenšují. Jednotlivá stabilní řešení se pak prolínají.

3-periodické řešení

Existuje limitní hodnota rC, ve kterém jsou nestabilní všechny

periodické řešení periody 2n.

Např. 3. iterace se 3 rovnovážnými stavy, kde = 1 (tečná bifurkace)

Pro případ, kdy se stabilní řešení stavy překrývají, vznikne liše-periodické řešení.

Sarkovsiiho theorém: Jestliže existuje liše-periodické řešení pro parametr rC, pak existuje chaotické

řešení pro r > rC

Řešení ut modelu pro různá r

K

Nr

tt

t

eNN1

1

2-cyklé periodické řešení

8-cyklé periodické řešení

4-cyklé periodické řešení

chaotické chování

chaotické chování 3-cyklé periodické řešení

Časová iterace diskrétného modelu

Iterativní diagram diskrétního modelu pro 1,9 < r < 3,0

Pro r2 < r < r4 … řešení ut osciluje mezi dvěma body A a B

Pro r4 < r < r8 … ut představuje 4-periodické řešení

Pro rC < r < rP … řešení je chaotické (neperiodické)

Pro r > rP … pravidelné periodické řešení, po něm následuje opět

neperiodické, atd.

Stejný sled bifurkace se opakuje ve fraktálním duchu

)1(1 tttt xrxxx

Výzkum chaosu vedl k zajímavým zjištěním:

Jestliže r2, r4, … r2n, … jsou sledem periodického zdvojení

bifurkačních hodnot, pak

univerzální konstanta

Jestliže pro některé ut a libovolné a existuje liché číslo n

takové, že platí , pak existuje liché periodické řešení, které předznamená chaos.

...66920,4)1(2)2(2

2)1(2

lim

nn

nn

n rr

rr

)( tuf

),(),( rufuruf tttn

Stabilita, periodické řešení a bifurkace (analyticky)

Parametr r mění řešení modelu

Zvyšováním r bifurkace vyšší period. řešení chaotické řešení

Rovnovážné body jsou řešením u*(r)

Zkoumáním lineární stability u*

Substitucí tohoto do a převedením do Taylorovy řady:

),(1 rufu tt rC = (– 1, 1)

),( ** rufu

tt vuu *1tv

)()()()( 2***1

*tttt vOufvufvufvu

)(1 tt ufu

)( ** ufu ,

ttt vufvv .)( *1 )( *uf … vlastní hodnota 1.iterace v u*

Řešení: pro když

Čili u* je jestliže

00vv t

t t

1

1

1)(

1)(1*

*

uf

uf

Pro u* STABILNÍ každá malá perturbace z u* slábne k nule

1)(0 * ufa) monotónně když

b) s klesajícími oscilacemi 0)(1 * uf

1)( * uf

1)( * uf

Pro u* NESTABILNÍ každá malá perturbace z u* roste

a) monotónně když

b) s rostoucími oscilacemi

Vezměme si příklad modelu , r > 0 )1(

1tur

tt euu

Rovnovážné stavy u* = 0 nebo u* = 1

Odpovídající vl. hodnoty pro r > 0 u* = 0 nestabilní

u* = 1 stabilní pro 0 < r < 1 s monotónní konvergencí stabilní pro 1 < r < 2 s klesajícími oscilacemi

u* = 1 nestabilní pro r > 2 s rostoucími oscilacemi

)1( *

1 ure

1)0( ref

rf 1)1(

r = 2 … první bifurkační hodnota v u* = 1 se rozvětví

Pro malé dostaneme Taylorovým rozvojem

přepíšeme do , kde

Z toho lze vyvodit: 4-periodické řešení od r4 2,45

6-periodické řešení od r6 2,54

chaotické chování pro r > rC 2,57

tu1 )]1(1[1 ttt uruu

]1[)1(1 ttt UUrU r

urU t

t

1

.

vysoká citlivost řešení na malou odchylku r

Pro t.iteraci u0 platí

Iterací , vytvoříme množinu

Bod je periodický s periodou m čili m-periodický, jestliže platí

a současně pro

Body u0, u1, …, um-1 tvoří formu m-cyklu

Bod u0 nazýváme vztažným bodem zobrazení fm

)( 0ufu tt

,...,, 210 uuu)()( 01

1 ufufu iii

,...2,1,0i

00 ),( uruf m 00 ),( uruf i 1,...,2,1,0 mi

Stabilita vztažného bodu u0

Pro rovnovážný stav u* to bylo jednoduše

Pro m-cyklus bodů u0, u1, …, um-1 uvažujme pro zjednodušení

Vlastní hodnotu m m-cyklu definujeme jako

, a tak

forma nezávislá na i

)( *uf

),(),( rufruF m ),(),( 1 rufruG m

i

m

m uuu

ruf

),( 1...10 mi

),()),((),( ruGruGfruF iii ),(),( 1 ruGruf ii

i

im

i uuu

rufruf

),(

),(1

1

1

0

),(m

iim ruf

i

m

m uuu

ruf

),(

Shrnutí:

Bifurkace nastává v hodnotě parametru r0 tehdy, když se

kvalitativně mění charakter řešení pro parametr r < r0 a r > r0. Z

předchozí analýzy očekáváme, že se změnou parametru r mění určité periodické řešení v řešení s jinou periodou. Jakmile je možné liché periodické řešení, pak dle Sarkovskeho teorému lze připustit řešení s libovolnou periodou, což implikuje chaos.

V hodnotě bifurkačního parametru přechází hodnota vlastního parametru stabilní rovnovážné polohy 1 nebo –1.

Pro hodnotu –1 se bifurkace nazývá vidlová Pro hodnotu 1 se bifurkace nazývá tečná

Za použití výkonných počítačů můžeme vypočítat všechny hodnoty každé rovnovážné polohy a následně určit bifurkační hodnoty r.

Diskrétní modely s prodlevou Doposud rozebírané modely jsou vhodné jen pro druhy se zanedbatelným dospívajícím obdobím (např. hmyz)

Nyní uvažujme časové období T potřebné k sexuální zralosti (např. u velryb). Tato prodleva vede ke studiu obecného diferenčního modelu formy

Jednoduchý model se zahrnutou prodlevou: , r > 0

Rovnovážné stavy u* = 0 nebo u* = 1

Stabilitu určíme pomocí malé perturbace z rovnovážné polohy (nelze použít metodu vl. hodnot - derivaci) Pro u* = 0

, r > 0 vt pro t u* = 0 nestabilní

),(1 Tttt uufu

)1( 11

t

tt

ureuu

)1( *

1 ure

tt vuu *tt vu 0 1tv

revrv

rev

rve

rev

vrevv ttt

tt

ttt .)1.(...

)1(. 1

111

1

re

Zkoumáním lineární stability u* = 1 Opět substitucí tohoto do a převedením do Tayl. řady:

Hledáme řešení diferenční rovnice ve formě vt = zt z2 – z + r = 0,

dá 2 hodnoty z1 a z2, kde pro

pro

Řešení je pak lin. kombinací , kde A a B jsou libovolné konstanty.

tt vu 1 1tv

)1)(1()1(1 11

1

ttt

tt rvvrv

evv 011 ttt rvvv

rz 4112

12,1

4

1r

iez

2,1 4

1r2

1r 14tan 1 r

tt vuu *

),(1 Tttt uufu

011 ttt rvvv ttt BzAzv 21

Jestliže 0 < r < 1/4 z1 a z2 jsou reálná 0 < z1 < 1, 0 < z2 < 1

vt 0 pro t u* = 1 je lineárně stabilní rovnovážná poloha.

Pro malou perturbaci monotónní.

Jestliže r > 1/4 z1 a z2 jsou komplexní,

Rovněž

12 zz

rzzz 22

121

Řešení je pak , aby bylo řešení reálné

Reálné řešení:

pro 1/4 < r < 1 je vt 0 pro t u* = 1 je stabilní

(A ani cos na stabilitu nemají vliv) pro r > 1 je vt pro t u* = 1 je nestabilní

kolem kritického bodu rC = 1 pro r 1 je

Poněvadž / 3 pro r 1

6-cyklé periodické řešení

121 zz

ttt zBAzv 11 AB

)cos(2 tAv tt Aarg14tan 1 r

33tan 1

11 z

)3/cos(2 tAvt

Řešení diferenční rovnice s prodlevou pro 3 hodnoty r > 1.

a) 6-periodické řešení, bifurkuje mimo stacionární stavy v r = rC

b) Známky 6-cyklu stále existují (nepravidelnost)

c) Pravidelnost 6-periodického řešení zcela vymizela (náznak chaosu)

Prodleva zapříčiňuje destabilizační efekt

V předešlém případě pro r = 2 řešení bifurkuje do 2-periodického řešení V pozdějším úloze s prodlevou pro r = 1 bifurkace do 6-periodického řešení

Čím větší prodleva, tím větší destabilizační efekt velké populační výkyvy možná cesta k zániku

Aplikace poznatků na lov velrybVýznam: řízení velrybí populace vůči přemnožení nebo vyhynutí Vyžaduje znát dynamiku velrybí populace a její ekologii

Populační model s prodlevou pro velrybí populaci

0 < < 1 …faktor úmrtnosti velryb R… faktor porodnosti T…doba dospívání ke schopnosti reprodukce mláďat (5 – 10 let)

Předpokládáme: poměr samců a samic je 1, stejnou úmrtnost obou pohlaví

Závislost členu R(Nt-T):

K …počet neulovených jedinců z N možných pro K = N (žádný ulovený) je hranatá závorka je nulová

)()1(1 Tttt NRNN = co přežili + noví jedinci

zT

K

NQPNNR 1)1(

2

1)(

pro K < N je hranatá závorka je záporná z …míra obtížnosti, se kterou jde hustota lovu registrovat Q …maximální porodnost pro nízký počet jedinců P …plodnost na jednu samici (1-)T.N … část přeživších do dospělosti po potřebných T letech 1/2 … polovina velryb jsou samice

Parametry , T a P jsou závislé!

Rovnovážný stav:

Definované h závisí na plodnosti P do úmrtnosti a prodlevě T

Model přeškálujeme

kde ,

KNNNN Tttt 1* hPT )1(

2

1

K

Nu tt )1(1)1(1

zTtTttt uqhuuu

P

Qq Ph T)1(

2

1

Rovnovážné stavyu* = 0 …nestabilní stav u* = 1 …zjistíme perturbací kolem rovnovážné polohy a linearizací

Řešení vt uvažujeme ve tvaru st získáme charakteristickou rovnici

Pro přejde rovnovážný stav do nestabilního

Pro přejde rovnovážný stav do stabilního

tt vu 1 Tttt vqzhvv )1()1(1

0)1()1(1 qzhss TT

1s

1s

Stabilita či nestabilita je závislá už na 4 parametrech (, T, h a qz)

Parametry jsou závislé samý na sobě. Určení zdali je rovnovážná poloho stabilní, je velice složité.

Hospodářský model rybářství Užitečný při vyhodnocení různých strategiích s ohledem na optimalizaci ekonomického přínosu a k jeho udržení. Je vhodný na jakékoli obnovitelné zdroje, které se loví

Hustota populace se řídí vztahem

ht… úlovek populace v čase t

Jaké je maximum trvalé biologické produkce? Jaké je maximum ekonomických výnosů?

Z rovnováhy: Nt = Nt+1 = N*, ht = h*

YM …maximum trvalých rovnovážných výnosů, když N* = NM

ttt hNfN )(1

*** )( NNfh

0*

*

N

h 1)( * Nf MMM NNfY )(0MY

… bez uvažování lovu)(1 tt NfN

… s uvažováním lovu

Cíl hospodářské strategie: Udržování populace tak, aby se získalo maximum výnosu YM (pro udržení rovnovážného stavu)

Ale je těžké znát přesný stav rybí populace. Známý je současný výnos a jak velké úsilí to dalo.

EM …celkové úsilí odpovídající Ym

c…parametr míry úsilí při lovu

Aplikace na model: , 0 < a < b

Substitucí:

Vyloučením NM:

M

MM N

NfcE

)(ln1

t

ttt Na

bNNfN

)(1

2)()(1

MM Na

abNf

)( 2/12/12/1 abaNM

MM

MM N

Na

bNY

MM Na

bcE ln1

1

)()( MMM

cEea

cEbeY

a) Vztah YM a EM = klíčové hledisko hospodářské strategie.

Když se přírůstkem úsilí sníží výnosy maximum trvalých výnosů překročeno, úsilí má být omezeno, aby se populace mohla obnovit.

b) Uvažujeme maximální ekonomické přínosy včetně ceny za lov a náklady na úsilí. Model s výrazy pro ekonomický návrat p … cena za jednotku výnosu k … cena za jednotku úsilí dostaneme křivku pro maximum výnosu R jako funkci úsilí E

MM kEpYR

Ekologický význam a varování Důležité: porozumět významným řídícím vlastnostem a schopnost předpovídat možné vývoje vyplývajících ze změn parametrů

Každý model závisí na spoustě parametrů. Některé z nich jsou relevantní. Z grafického řešení jsme se dozvěděli, jak se mění, řešení, když se měníme parametry. Pokud bychom vytvořili dokonalý model, mohli bychom předpovědět přesný vývoj populace. To se však nepodaří. Hustota populace je vždy omezená.

Z obrázku lze vyčíst, že vývoj populace je ohraničen mezi dvěmi hranicemi, N populace je ohraničená hodnotou Nmax a Nmin, ať startujeme

z kterýchkoli počátečních podmínek.

Když se velikost populace N dostane do nízkých hodnot hrozí vyhynutí. Nesmí být < 1. Populace může vymřít, kvůli velkému stavu populace (přemnožení). Ten zjistíme z maxima křivky

, který najdeme pomocí Nm

… populace vymře

To aplikujeme na model:

)(1 tt NfN 0tdN

df)(max mNfN

)())(()( 2maxmin mm NfNffNfN

Parametry K a r rozhodují o globální nestabilitě.

Řekněme, když r = 3,5 a K < 1600, populace časem zanikne.

K

Nr

eNNfNt

ttt

1)(1

0)( tNfr

KNm

1max )( r

m er

KNfN

1

min

12))((

r

m

ere

r

KNffN

1)(2min mNfN

0 mNNdN

df

1

12 1

rer

er

K

- další skupinou modelů

vymření

Rovnovážné stavy … Nt = 0, NC a N*

Nt = 0 … stabilní, NC …nestabilním, N* … podle

Oblast Nt < NC … nazývána predation pit (propast)

)( *Nf

Modely s Aleeo efektem

CNNf )( 02 0tN Ct

m NNf )( Nm

Bohatství chování modelů z výsledné nelinearity těchto modelů.