Modelowanie komputerowe w ochronie środowiska Modele ekologiczne
Wybrane modele ekologiczne
oraz metody rozwiązywania równań różniczkowych zwyczajnych
Na podstawie:
�� Urszula Foryś, Matematyka w biologii, Wydawnictwa Naukowo-Techniczne, Warszawa
2005.
�� Janusz Uchmański, Klasyczna ekologia matematyczna, Wydawnictwo Naukowe PWN,
Warszawa 1992.
�� wykład internetowy "Quantitative Population Ecology"
http://www.ento.vt.edu/~sharov/PopEcol/popecol.html
�� Numerical Recipies, roz. 16.1  Runge-Kutta Method
http://www.nrbook.com/a/bookfpdf/f16-1.pdf
Równanie Malthusa  równanie wzrostu wykładniczego
Zakłada się, że w danym środowisku występuje tylko jeden gatunek, oraz że zasoby środowiska
są nieograniczone. Populacja jest jednorodna, osobniki nie umierają.
Równanie Malthusa wyraża się jako:
dN(t)
=� r N(t)
(1) ,
dt
N(t)
gdzie: - zagęszczenie osobników w chwili t;
r >0  współczynnik rozrodczości gatunku.
Aby równanie miało jednoznaczne rozwiązanie niezbędne jest podanie warunku początkowego:
N(0) =� N0
.
Równanie (1) jest równaniem różniczkowym zwyczajnym. Jego rozwiązaniem jest funkcja
wykładnicza
N(t) =� N0 ert
Uwaga:
W modelu ciągłym nie rozpatruje się liczebności populacji, lecz jej zags zczenie, czyli liczbę
N(t)
osobników przypadających na jednostkę powierzchni. Niekiedy można interpretować
jako mas danej populacji.
Proces urodzin i śmierci
Wprowadzając współczynnik śmiertelności populacji s>0 można zbudować model procesu
urodzin i śmierci:
dN(t)
=� (r -� s) N(t)
(2)
dt
Rozwiązaniem tego równania jest funkcja wykładnicza:
N(t) =� N0 e(r-�s) t
r -� s
Zmiany liczebności populacji zależą od znaku współczynnika .
Model wzrostu wykładniczego Proces urodzin i śmierci
30
250
25
200
20
150
15
100
10
50
5
0 0
0 1 2 3 4 5 6 0 1 2 3 4 5 6
czas czas
r=0,1 r=0,4 r=0,6 r-s>0 r=s r-s<0
2007-02-13 1
Modelowanie komputerowe w ochronie środowiska Modele ekologiczne
Model logistyczny
W modelu logistycznym wprowadza się ograniczoność zasobów środowiska (pojemność
środowiska), co powoduje wystąpienie konkurencji między osobnikami (np. zdobywanie
pożywienia, terytorium, itp.).
Równanie logistyczne:
dN(t) N
=� r N(t) (1-� )
(3) ,
dt K
N(t)
gdzie: - zagęszczenie osobników w chwili t;
r >0  współczynnik rozrodczości gatunku;
K>0  pojemność środowiska.
r
a =�
Pojemność środowiska K definiuje się w oparciu o współczynnik konkurencji a, gdzie .
K
Zatem współczynnik konkurencji jest wprost proporcjonalny do współczynnika rozrodczości r i
odwrotnie proporcjonalny do pojemności środowiska K.
Rozwiązanie równania logistycznego:
K
N(t) =�
K
1+� ( -�1) e-�rt .
N0
t �� Ą� N(t) �� K
Przy co oznacza, że liczebność populacji stabilizuje się na poziomie
odpowiadającym pojemności środowiska.
Model logistyczny, K=150
300
250
200
150
100
50
0
0 1 2 3 4 5
N0=10 N0=50 N0=300 N0=150
Dwuwymiarowe modele ekologiczne  model Lotki-Volterry
Model dotyczy zmian liczebności populacji dwóch gatunków: ofiar (ang. prey) oraz
drapieżników (ang. predator) żywiących się osobnikami pierwszego gatunku.
Niech: H(t)  zagęszczenie ofiar w chwili t,
P(t)  zagęszczenie drapieżników w chwili t.
Wtedy populacja dwóch gatunków jest opisana układem równań różniczkowych zwyczajnych:
dH (t)
��
=� r H (t) -� a H (t) P(t)
��
dt
(4)
��
dP(t)
��
=� b H (t) P(t) -� m P(t)
�� dt
gdzie: r  współczynnik rozrodczości gatunku ofiar
a  współczynnik skuteczności polowań
b  współczynnik rozrodczości drapieżników (na jednostkę upolowanej ofiary)
m  współczynnik śmiertelności drapieżników
Rozwiązania H(t) i P(t) są funkcjami okresowymi przesuniętymi w fazie.
2007-02-13 2
Modelowanie komputerowe w ochronie środowiska Modele ekologiczne
Portret fazowy układu drapieżca-ofiara
140
120
25
100
20
80
60
15
40
10
20
0
5
0 50 100 150 200 250 300 350
H P 0
0 20 40 60 80 100 120 140
Wyniki uzyskano dla następujących wartości
zagęszczenie ofiar H
parametrów:
r=0,1 a=0,01, b=0,001, m=0,05
H =50, P =20
0 0
Modyfikacje modelu Lotki-Volterry
�� Model drapieżnik  ofiara z kryjówkami dla ofiar
Przyjmuje się dodatkowe założenie, że pewna cześć ofiar jest niedostępna dla drapieżników
poprzez wykorzystanie kryjówek. Zakłada się, że liczba ukrywających się ofiar jest stała i
oznaczana symbolem z.
Populacja obu gatunków jest opisana za pomocą następująco zmodyfikowanego równania
Lotki-Volterry:
dH
��
=� r H -� a (H -� z) P
��
dt
(5)
��
dP
��
=� b (H -� z) P -� m P
�� dt
gdzie: r  współczynnik rozrodczości gatunku ofiar
a  współczynnik skuteczności polowań
b  współczynnik rozrodczości drapieżników (na jednostkę upolowanej ofiary)
m  współczynnik śmiertelności drapieżników
90 160
80
140
70
120
60
100
50
80
40
60
30
40
20
20
10
0 0
0 50 100 150 200 250 300 0 50 100 150 200 250 300
H P H P
Wyniki uzyskano dla następujących wartości Wyniki uzyskano dla następujących wartości
parametrów: parametrów:
r=0,1 a=0,01, b=0,001, m=0,05, z=20 r=0,1 a=0,01, b=0,001, m=0,05, z=20
H =50, P =20 H =150, P =40
0 0 0 0
2007-02-13 3
zagęszczenie drapieżców
Modelowanie komputerowe w ochronie środowiska Modele ekologiczne
Portret fazowy - różne wartości początkowe Portret fazowy - różne wartości początkowe (powiększenie)
50 20
45
40 18
35
30 16
25
20 14
15
10 12
5
0 10
0 20 40 60 80 100 120 140 160 40 50 60 70 80 90 100
H - ofiary H - ofiary
H0=150, P0=40 H0=50, P0=20 H0=150, P0=40 H0=50, P0=20
�� Model drapieżnik  ofiara z ograniczoną pojemnością środowiska dla ofiar
Przyjmuje się dodatkowe założenie, że wśród ofiar występuje konkurencja spowodowana
ograniczonymi zasobami środowiska (na przykład walką o pożywienie) i w związku z tym w
równaniu Lotki-Volterry zastępuje się wykładniczy model wzrostu populacji ofiar modelem
logistycznym. Pojemność środowiska jest oznaczana symbolem K.
Populacja obu gatunków jest opisana za pomocą następująco zmodyfikowanego układu równań:
dH H
��
=� r H (1-� ) -� a H P
��
dt K
(6)
��
dP
��
=� b H P -� m P
�� dt
gdzie: r  współczynnik rozrodczości gatunku ofiar
a  współczynnik skuteczności polowań
b  współczynnik rozrodczości drapieżników (na jednostkę upolowanej ofiary)
m  współczynnik śmiertelności drapieżników
m m
K Ł� K >�
Warto porównać rozwiązania uzyskane dla oraz dla .
b b
Model drapieżca-ofiara; wzrost ofiar według
Portret fazow y układu drapieżca-ofiara
modelu logistycznego; liniowy wzrost drapieżników.
25
60
50
20
40
30
15
20
10
10
0
5
0 50 100 150 200 250 300
0
H P
0 20 40 60
Wyniki uzyskano dla następujących wartości
zagęszczenie ofiar H
parametrów:
r=0,1 a=0,01, b=0,001, m=0,05, z=20
H =50, P =20, K=30, a zatem K0 0
2007-02-13 4
P - drapieżcy
P - drapieżcy
zagęszczenie drapieżców
Modelowanie komputerowe w ochronie środowiska Modele ekologiczne
Model drapieżca-ofiara; wzrost ofiar według
Portret fazow y układu drapieżca-ofiara
modelu logistycznego; liniowy wzrost drapieżników.
25
80
70
20
60
50
40
15
30
20
10
10
0
5
0 50 100 150 200 250 300
0
H P
0 20 40 60 80
Wyniki uzyskano dla następujących wartości
zagęszczenie ofiar H
parametrów:
r=0,1 a=0,01, b=0,001, m=0,05, z=20
H =50, P =20, K=100, a zatem K>m/b
0 0
Mutualizm
Układy mutualistyczne między gatunkami powstają gdy każdy z gatunków czerpie korzyści z
aktywności partnera (owady  rośliny owadopylne). W klasycznej ekologii matematycznej
korzyści te są wyrażane przyrostami zagęszczeń gatunków.
�� Mutualizm: Model 1
Populacja obu gatunków jest opisana za pomocą następujących równań:
dH H
��
=� r1 H (1-� ) +� d�1 H P
��
�� dt K1
(7)
��
dP P
��
=� r2 P(1-� ) +� d�2 H P
�� dt K2
��
gdzie: r , r  współczynniki rozrodczości gatunków
1 2
K , K  pojemności środowiska dla obu gatunków
1 2
d�1, d�2  intensywność dobroczynnych oddziaływań między gatunkami.
Jeśli oddziaływania między gatunkami są silne w porównaniu z oddziaływaniami
wewnątrzgatunkowymi, to wtedy mutualistyczne oddziaływanie powoduje nieograniczony
wzrost zagęszczeń obu gatunków.
Gdy oddziaływania między gatunkami są słabsze, to wtedy zagęszczenia obu gatunków
zbiegają asymptotycznie do wartości równowagowych, których wartości, wskutek
dobroczynnego oddziaływania, są większe od pojemności środowiska dla poszczególnych
gatunków.
�� Mutualizm: Model 2
Bardziej realistyczny jest model uwzględniający ograniczenia zysków wynikających z
mutualizmu, w którym w wyniku pozytywnego oddziaływania drugiego gatunku zwiększa się
pojemność środowiska gatunku pierwszego. Można to opisać w postaci następującego układu
równań (May, 1981):
dH H
��
=� r1 H (1-� )
��
�� dt K1 +� a�P
(8)
��
dP P
��
=� r2 P(1-� )
�� dt K2 +� b�H
��
gdzie: r , r  współczynniki rozrodczości gatunków
1 2
K , K  pojemności środowiska dla obu gatunków
1 2
a�, b�  intensywność dobroczynnych oddziaływań między gatunkami.
2007-02-13 5
zagęszczenie drapieżców
Modelowanie komputerowe w ochronie środowiska Modele ekologiczne
W modelu tym zakłada się, że każdy z gatunków, jeśli jest izolowany, rozwija się zgodnie z
modelem wzrostu logistycznego. Natomiast oddziaływania drugiego gatunku przejawia się
modyfikacją pojemności środowiska.
a�, b�
Dla dużych wartości parametrów następuje nieograniczony wzrost zagęszczeń obu
a� b� <� 1
gatunków. Natomiast jeśli oddziaływania międzygatunkowe są takie, że zachodzi , to
wtedy istnieje położenie równowagi. Wartości równowagowe H* i P* wyrażają się wzorami:
K1 +� a�K2 K2 +� b�K1
*
H =� P* =�
1-� a�b� 1-� a�b�
Mutualizm - portret fazow y
Mutualizm - zwiększenie pojemności środowiska
35 35
30 30
25 25
20 20
15
15
10
10
5
5
0
0
0 20 40 60 80 100 120
0 10 20 30 40
zagęszczenie gatunku H
H P
Wyniki uzyskano dla następujących wartości
parametrów:
a� b�
r =0,1 r =0,1, K =10, K =20, =0,5, =0,4,
1 2 1 2
H =30, P =5.
0 0
Portret fazowy: mutualizm: zwiększenie pojemności
Mutualizm - zwiększenie pojemności środowiska
środowiska
nieograniczony wzrost
70
1400
60
1200
50
1000
40
800
30 600
20 400
200
10
0
0
0 20 40 60 80 100 120
0 5 10 15 20 25 30 35
zagęszczenie gatunku H
H P
Portret fazowy dla następujących wartości
a� b�
r =0,1 r =0,1, K =10, K =20, =2, =1,2,
1 2 1 2
a�
parametrów: r =0,1 r =0,1, K =10, K =20, =0,5,
1 2 1 2
H =30, P =5.
0 0
b� =0,4, oraz różnych warunków początkowych H i
0
P
0.
�� Mutualizm: Model 3
Można wprowadzić jeszcze inne założenie, opisujące mutualizm, zakładając, że każdy z
gatunków bardziej odczuwa dobrodziejski wpływ partnera gdy jego własne zagęszczenie jest
małe. Założenie to można opisać na przykład za pomocą następujących równań (modyfikacja
układu (7)):
dH H
��
1
=� r1 H (1-� ) +� d�1 H P e-�a� H
��
�� dt K1
(9)
��
dP P
2
��
=� r2 P(1-� ) +� d�2 H P e-�a� P
�� dt K2
��
gdzie: r , r  współczynniki rozrodczości gatunków
1 2
K , K  pojemności środowiska dla obu gatunków
1 2
d�1, d�2  intensywność dobroczynnych oddziaływań między gatunkami.
a�1, a�2  parametry (dodatnie).
2007-02-13 6
zagęszczenie gatunku P
zagęszczenie gatunku P
Modelowanie komputerowe w ochronie środowiska Modele ekologiczne
Człon opisujące związki mutualistyczne dla pierwszego równania jest mały jeśli H jest duże,
natomiast dla małych H jego wartość wzrasta. Analogicznie zachowuje się człon
mutualistyczny w drugim równaniu, tym razem w zależności od wartości P. Tak opisany układ
jest zawsze stabilny, niezależnie od wartości parametrów.
Mutualizm - oddziaływanie wykładnicze Portret fazowy - mutualizm - oddziaływanie wykładnicze
14 14
12
12
10
10
8
8
H
6
6 P
4
4
2
2
0
0 2 4 6 8 10 12
0
0 10 20 30 40
zagęszczenie gatunku H
a� b�
r =0,1 r =0,1, K =10, K =5, =2, =3, H =5,
1 2 1 2 0
P =12.
0
Konkurencja
Rozważa się wzajemne oddziaływanie dwóch gatunków, z których każdy przebywając w
izolacji podlegałby wzrostowi logistycznemu. Jeśli obydwa gatunki występują równocześnie to
pojawia się konkurencja, wpływająca hamująco na wzrost zagęszczenia każdego z gatunków.
Przyjmuje się, że ograniczenie szybkości przyrostu w wyniku konkurencji jest proporcjonalne
do zagęszczenia każdego gatunku.
Powyższe założenia w postaci matematycznej wyglądają następująco:
dH H
��
=� r1 H (1-� ) -� d�1 H P
��
�� dt K1
(10)
��
dP P
��
=� r2 P(1-� ) -� d�2 H P
�� dt K2
��
gdzie: r , r  współczynniki rozrodczości gatunków;
1 2
K , K  pojemności środowiska dla obu gatunków;
1 2
d�1  miara podatności gatunku H na konkurencję ze strony gatunku P;
d�2  miara podatności gatunku P na konkurencję ze strony gatunku H;
Również ten model został zaproponowany przez Lotkę i Volterrę.
W większości przypadków w wyniku konkurencji następuje wyginięcie jednego gatunku.
Natomiast trwałe współwystępowanie obu gatunków jest możliwe tylko wtedy, gdy
oddziaływania między gatunkami są słabe.
Konkurencja
Portret fazowy - konkurencja
45 45
40
40
35
35
30
30
25
25 H
20
20
P
15
15
10
10
5
5
0
0 5 10 15 20 25 30 35
0
0 10 20 30 40 zagęszczenie gatunku H
Wyginięcie gatunku H
r =0,1 r =0,1, K =50, K =40, d�1 =0,03, d�2 =0,06,
1 2 1 2
H =30, P =20.
0 0
2007-02-13 7
zagęszczenie gatunku P
zagęszczenie gatunku P
Modelowanie komputerowe w ochronie środowiska Modele ekologiczne
Konkurencja Konkurencja - współistnienie gatunków
60 40
35
50
30
40
25
H H
30 20
P P
15
20
10
10
5
0 0
0 20 40 60 80 0 50 100 150
Wyginięcie gatunku P Współwystępowanie gatunków
r =0,1 r =0,1, K =50, K =40, d�1 =0,02, d�2 =0,06, r =0,1 r =0,1, K =50, K =40, d�1 =0,002, d�2 =0,006,
1 2 1 2 1 2 1 2
H =30, P =20. H =30, P =20.
0 0 0 0
2007-02-13 8
Modelowanie komputerowe w ochronie środowiska Modele ekologiczne
Metody numeryczne rozwiązywania równań różniczkowych zwyczajnych  podstawowe
wzory
Numeryczne rozwiązanie równania różniczkowego zwyczajnego:
dy(t)
=� f (t, y(t))
dt
z warunkiem początkowym
y(t0 ) =� y0
yn =� y(tn ) tn-�1 tn
polega na wyznaczeniu dyskretnych wartości , n=1,2,.... Kolejne punkty i
znajdują się w odległości h; h jest nazywane krokiem całkowania.
Metoda Eulera
Najprostszą (ale najmniej dokładną) metodą jest metoda Eulera.
yn =� yn-�1 +� h f (tn-�1, yn-�1)
W szczególności:
y1 =� y0 +� h f (t0 , y0 )
,
y2 =� y1 +� h f (t1, y1)
, itd.
Metody Runge-Kutty
Znacznie lepszą jakość rozwiązania przybliżonego uzyskuje się poprzez:
k1 =� h f (tn-�1, yn-�1)
k1
k2 =� h f (tn-�1 +� 0,5h, yn-�1 +� )
2
yn =� yn-�1 +� k2
Jest to metoda Runge-Kutty drugiego rzędu.
Najpopularniejszą i bardzo często wykorzystywaną w praktyce metodą jest metoda Runge-
Kutty czwartego rzędu.
k1 =� h f (tn-�1, yn-�1)
k1
k2 =� h f (tn-�1 +� 0,5h, yn-�1 +� )
2
k2
k3 =� h f (tn-�1 +� 0,5h, yn-�1 +� )
2
k4 =� h f (tn-�1 +� h, yn-�1 +� k3)
k1 k2 k3 k4
yn =� yn-�1 +� +� +� +�
6 3 3 6
Dla naszych potrzeb wystarczająca będzie metoda rzędu drugiego.
Podane metody można stosować również do układów równań różniczkowych zwyczajnych.
Wszystkie wyniki przedstawione w tym opracowaniu zostały wyznaczone w sposób
przybliżony z wykorzystaniem metody Rungego-Kutty drugiego rzędu. Obliczenia i rysunki
wykonano w arkuszu kalkulacyjnym Excel.
2007-02-13 9