Zagadnienia

2.1 Równanie logistyczne — model Verhulsta
2.2 Dyskretne równanie logistyczne

2. Modelowanie pojedynczej populacji II

2.1. Równanie logistyczne — model Verhulsta

Jak sprawdziliśmy, model Malthusa z emigracją nie rozwiązuje problemu ograniczonej pojemności siedliska. Wszystkie modele rozważane do tej pory bazują na modelu Malthusa, przy czym są liniowe. Teraz omówimy model zaproponowany ok. 1840 r. przez belgijskiego matematyka Pierre'a F. Verhulsta w [16], który jest prawdopodobnie pierwszym znanym w biologii populacyjnej modelem nieliniowym.

Model ten powstał w wyniku dyskusji grupy uczonych nad modelem Malthusa. Większość osób biorąca udział w tej dyskusji zgadzała się co do tego, że istnieją naturalne procesy, które hamują postulowany przez Malthusa geometryczny (nieograniczony) przyrost liczebności populacji. Według Verhulsta przyrost liczebności hamuje np. konkurencja o zasoby siedliska, która wiąże się w oczywisty sposób z ograniczonością tych zasobów. W związku z tym w równaniu (1.2) — które w dalszym ciągu będziemy utożsamiać z równaniem (1.4) dla $r_{n}>0$ — należy uwzględnić składnik opisujący konkurencję. W tym kontekście mówimy o konkurencji wewnątrzgatunkowej w odróżnieniu od konkurencji zewnątrzgatunkowej, występującej pomiędzy osobnikami różnych gatunków. Trzeba się więc zastanowić, jak reprezentować tę konkurencję w języku matematyki. Wszystkie modele biologiczne powstające w XIX (czy nawet w pierwszej połowie XX wieku) miały swoje korzenie w ugruntowanych już modelach fizycznych. Zatem nie należy się dziwić, że konkurencję między osobnikami tego samego gatunku Verhulst opisał w taki sam sposób jak losowe zderzenia cząsteczek gazu elementarnego. Liczba takich zderzeń jest proporcjonalna do kwadratu liczby cząstek $N^{2}$ , a biorąc pod uwagę, że osobniki mogą konkurować ze sobą tylko wtedy, gdy się bezpośrednio zetkną można przyjąć, że konkurencja jest proporcjonalna do kwadratu liczebności, czyli równanie (1.2) zamienia się na

$\dot{N}(t)=rN(t)-bN^{2}(t),$

(2.1)

gdzie $b$ oznacza współczynnik konkurencji wewnątrzgatunkowej i odzwierciedla ,,szkodliwość” tej konkurencji dla populacji $\mathcal{P}$ .

To samo równanie możemy otrzymać na podstawie innego modelu heurystycznego. Wydaje się oczywiste, że przyjęty w modelu Malthusa (1.4) stały współczynnik reprodukcji $r_{n}=\text{const}$ stanowi często zbyt duże uproszczenie. Zwykle współczynnik ten zależy od wielu czynników, w tym także bezpośrednio od liczebności populacji. Jeśli liczebność jest niewielka, to osobniki mają dobry dostęp do zasobów siedliska i mogą się bez przeszkód rozmnażać, ale przy wzroście liczebności ten dostęp zmniejsza się. W takiej sytuacji przyrost per capita (czyli względny przyrost na jednego osobnika w populacji)

$\frac{\dot{N}(t)}{N(t)}=\frac{\frac{dN}{dt}(t)}{N(t)},$

stały w modelu Malthusa (1.4), staje się funkcją zależną od $N(t)$ . Ponieważ według naszych założeń reprodukcja zmniejsza się wraz ze wzrostem liczebności populacji, to najprostsze matematyczne odzwierciedlenie takiej zależności stanowi liniowa funkcja malejąca. Wobec tego

$\frac{\dot{N}(t)}{N(t)}=h(N(t)),\quad\text{gdzie}\quad h(x)=r\left(1-\frac{x}{K}\right)$

i przy $b=\tfrac{r}{K}$ równanie (2.1) i równanie logistyczne

$\dot{N}(t)=rN(t)\left(1-\frac{N(t)}{K}\right)$

(2.2)

pokrywają się. Nazwa ,,równanie logistyczne” wiąże się z typowym logistycznym kształtem części rozwiązań równania (2.2). Równanie to nazywa się często także równaniem Verhulsta lub równaniem Verhulsta – Pearla, przy czym nazwisko biologa R. Pearla pojawia się w tej nazwie w związku z zastosowaniami równania (2.2) do opisu konkretnych populacji.

Przejdziemy teraz do analizy przebiegu rozwiązań równania (2.2). Jak poprzednio jest to równanie o zmiennych rozdzielonych, które możemy rozwiązać korzystając z rozkładu ułamka $\tfrac{1}{N\left(1-N/K\right)}$ na ułamki proste,

$\frac{1}{N\left(1-\frac{N}{K}\right)}=\frac{1}{N}+\frac{\frac{1}{K}}{1-\frac{N}{K}}.$

Naszym głównym zadaniem nie jest jednak poszukiwanie analitycznych rozwiązań budowanych modeli, tylko analiza własności i przebiegu tych rozwiązań. Zauważmy, że prawa strona równania (2.2) jest funkcją klasy $\mathbf{C}^{1}$ , co implikuje istnienie i jednoznaczność rozwiązań. Łatwo też sprawdzić, że dla $N_{0}\geq 0$ mamy $N(t)\geq 0$ dla $t\geq 0$ . Jedną z metod tej wykazania tej własności stanowi zapisanie równania (2.2) w równoważnej postaci całkowej

$\frac{\dot{N}}{N}=r\left(1-\frac{N}{K}\right)\ \Longrightarrow\ \int\limits _{{N_{0}}}^{{N(t)}}\frac{dN}{N}=\int\limits _{0}^{t}r\left(1-\frac{N(s)}{K}\right)ds\ \Longrightarrow\ N(t)=N_{0}\text{e}^{{r\int\limits _{0}^{t}\left(1-\frac{N(s)}{K}\right)ds}},$

skąd od razu wynika nieujemność. Z kolei nieujemność implikuje oszacowanie

$\dot{N}(t)\leq rN(t)\quad\Longrightarrow\quad N(t)\leq N_{0}\text{e}^{{rt}},$

czyli wzrost jest co najwyżej wykładniczy i rozwiązanie istnieje dla wszystkich $t\geq 0$ .

$Wykres funkcji $f(N)=rN\left(1-\tfrac{N}{K}\right)$ (po lewej) oraz odpowiadający mu portret fazowy równania~\eqref{logistyczne} (po prawej).$

Rys. 2.1. Wykres funkcji $f(N)=rN\left(1-\tfrac{N}{K}\right)$ (po lewej) oraz odpowiadający mu portret fazowy równania (2.2) (po prawej).

Zbadamy teraz zależność $(N,\dot{N})$ , która determinuje monotoniczność rozwiązań. Przyjrzyjmy się wykresowi funkcji $f(N)=rN\left(1-\tfrac{N}{K}\right)$ , por. rys. 2.1. Jest to parabola o miejscach zerowych $N=0$ i $N=K$ skierowana w dół. Wynika stąd, że $N\equiv 0$ i $\bar{N}\equiv K$ są rozwiązaniami stacjonarnymi. Jeśli $N(t)\in(0,K)$ , to $\dot{N}(t)>0$ , czyli rozwiązania dla warunku początkowego $N_{0}<K$ są rosnące i ograniczone, ponieważ z jednoznaczności rozwiązań wynika, że nie mogą przeciąć stacjonarnego rozwiązania $\bar{N}$ . Zatem jeśli $0<N_{0}<K$ , to

$\lim\limits _{{t\to\infty}}N(t)=K.$

Analogicznie jeśli $N_{0}>K$ (rozwijająca się w naturalnych warunkach populacja nie może przekroczyć granicy $K$ , jak pokazaliśmy powyżej, czyli taka sytuacja jest możliwa jedynie w przypadku introdukcji osobników), to $\dot{N}(t)<0$ , więc dla takich wartości początkowych liczebność maleje i ponownie asymptotycznie osiąga wartość $K$ . Jednak aby dokładniej przeanalizować przebieg rozwiązań musimy znać nie tylko wartość pierwszej pochodnej, ale także drugiej, która determinuje wypukłość/wklęsłość funkcji $N(t)$ .

Obliczmy zatem

$\ddot{N}(t)=\frac{d}{dt}\dot{N}(t)=\frac{d}{dt}\left(rN(t)\left(1-\frac{N(t)}{K}\right)\right)=$

$r\dot{N}(t)\left(1-\frac{N(t)}{K}\right)+rN(t)\left(-\frac{\dot{N}(t)}{K}\right)=r\dot{N}(t)\left(1-\frac{2N(t)}{K}\right).$

Widzimy, że jeśli rozwiązanie w pewnej chwili $t_{p}$ przyjmuje wartość $N(t_{p})=\tfrac{K}{2}$ , to zmienia się charakter przebiegu — funkcja $N(t)$ ma punkt przegięcia. Analizując znaki czynników $\dot{N}(t)$ i $1-\tfrac{2N(t)}{K}$ możemy wydzielić kilka obszarów, por. rys. 2.2

$N_{0}\in\left(0,\frac{K}{2}\right)$ , to $N(t)$ rośnie do $K$ przy $t\to\infty$ , więc istnieje $t_{p}$ , takie że $\ddot{N}(t_{p})=0$ i $N(t)$ jest wypukła dla $t<t_{p}$ (szybki wzrost, porównywalny z wykładniczym, dla małych zagęszczeń), a wklęsła dla $t>t_{p}$ (następuje spowolnienie wzrostu dla większych zagęszczeń);
$N_{0}\in\left(\frac{K}{2},K\right)$ , wtedy funkcja $N(t)$ rośnie, ale wzrost jest powolny, funkcja $N$ jest wklęsła;
$N_{0}>K$ , co (należy ponownie podkreślić) odpowiada introdukcji osobników, wtedy $N(t)$ maleje w sposób wypukły do $K$ , zatem spadek liczebności jest dość szybki.

$Przykładowe rozwiązania równania logistycznego~\eqref{logistyczne} dla różnych warunków początkowych$

Rys. 2.2. Przykładowe rozwiązania równania logistycznego (2.2) dla różnych warunków początkowych.

Widzimy więc, że dla $N_{0}\in\left(0,\tfrac{K}{2}\right)$ rozwiązanie ma kształt wydłużonej litery S, jak krzywa logistyczna, która dobrze odzwierciedla wiele procesów, przebiegających w początkowej fazie szybko i intensywnie, a po pewnym czasie ulegających wysyceniu i stabilizujących się na maksymalnym dla danego procesu poziomie.

2.2. Dyskretne równanie logistyczne

Często równolegle z modelem ciągłym rozważa się model dyskretny. Przedstawimy teraz dyskretną wersję równania (2.2). W tym celu przybliżymy pochodną $\dot{N}$ ilorazem różnicowym

$\dot{N}(t)\approx\frac{N(t+\Delta t)-N(t)}{\Delta t},$

założymy (analogicznie jak w przypadku modelu Malthusa), że $\Delta t$ jest równa hipotetycznej jednostce czasu $\Delta t=1$ i zastosujemy notację dyskretną $N(t)=N_{t}$ , $t\in\mathbb{N}$ . Otrzymujemy

$N_{{t+1}}=(1+r)N_{t}-\frac{r}{K}N_{t}^{2}=aN_{t}\left(1-\frac{N_{t}}{K_{1}}\right),$

(2.3)

gdzie $a=1+r$ i $K_{1}=\frac{1+r}{r}K$ . Zauważmy, że o ile z biologicznego punktu widzenia mamy w równaniu logistycznym dwa parametry: współczynnik rozrodczości i pojemność siedliska, to z matematycznego punktu widzenia jeden z nich możemy z łatwością wyeliminować. Niech $x_{t}=\tfrac{N_{t}}{K_{1}}$ . Dzieląc równanie (2.3) stronami przez $K_{1}$ dostajemy

$x_{{t+1}}=ax_{t}(1-x_{t}),$

(2.4)

czyli model jednoparametrowy. Zauważmy, że model ten ma sens, o ile kolejne wyrazy ciągu $\left(x_{t}\right)_{{t=0}}^{{\infty}}$ pozostają w przedziale $[0,1]$ , ponieważ reprezentują one zagęszczenia populacji w proporcji do pojemności siedliska $K_{1}$ . Wyrazy te obliczamy jako kolejne iteracje funkcji $F(x)=ax(1-x)$ , która jest funkcją kwadratową z wartością maksymalną $F\left(\tfrac{1}{2}\right)=\tfrac{a}{4}$ . Zatem aby wszystkie wyrazy pozostawały w przedziale $[0,1]$ musi zachodzić $0<a\leq 4$ . Będziemy więc badać dynamikę ciągu $\left(x_{t}\right)_{{t=0}}^{{\infty}}$ dla $a\in(0,4]$ i warunku początkowego $x_{0}\in[0,1]$ . Oczywiście jeśli $x_{0}=0$ albo $x_{0}=1$ , to $x_{t}=0$ dla $t\geq 1$ , więc interesujące są wartości $x_{0}\in(0,1)$ .

Jeśli ciąg $\left(x_{t}\right)_{{t=0}}^{{\infty}}$ ma granicę, to równa się ona któremuś ze stanów stacjonarnych. Wyznaczmy więc te stany. Niech $\bar{x}$ oznacza stan stacjonarny. Wtedy

$\bar{x}=a\bar{x}(1-\bar{x})\quad\Longrightarrow\quad\bar{x}=0\quad\text{albo}\quad\bar{x}_{1}=\frac{a-1}{a}.$

Widzimy więc, że dodatni stan stacjonarny istnieje, o ile $a>1$ . Lokalną stabilność tych stanów badamy za pomocą wartości własnych przekształcenia $dF(x)|_{{\bar{x}}}$ . W naszym przypadku $dF(x)=a(1-2x)$ . Jeśli $a\in(0,1)$ , to stan stacjonarny jest wyznaczony jednoznacznie $\bar{x}=0$ i $dF(0)=a$ , przy czym nierówność $0<a<1$ implikuje lokalną stabilność. Co więcej widzimy, że $x_{{t+1}}<x_{t}$ , ponieważ $x_{{t+1}}=F(x_{t})<x_{t}$ dla dowolnego $x_{t}\in(0,1)$ przy założeniu $a\leq 1$ . Mamy więc ciąg malejący i ograniczony, czyli zbieżny, zatem stan $\bar{x}=0$ jest globalnie stabilny w zbiorze $[0,1]$ .

Jeśli $a>1$ , to istnieją dwa stany stacjonarne, przy czym $\bar{x}=0$ staje się niestabilny, gdyż teraz $dF(0)=a>1$ . Dla dodatniego stanu stacjonarnego $\bar{x}_{1}=\tfrac{a-1}{a}$ mamy $dF(\bar{x}_{1})=2-a$ oraz

$|2-a|<1\Leftrightarrow a\in(1,3).$

Prostymi metodami, jak w przypadku zerowego stanu stacjonarnego pokazujemy, że dla $a\in(1,2]$ ciąg $\left(x_{t}\right)_{{t=0}}^{{\infty}}$ jest monotoniczny od pewnego miejsca, natomiast dla $a\in(2,3)$ możemy wyróżnić dwa podciągi monotoniczne od pewnego miejsca. Faktycznie, jeśli $a\in(1,2]$ , to dla $x_{t}\in(0,\bar{x}_{1})$ zachodzi $F(x_{t})>x_{t}$ , natomiast dla $x_{t}\in(\bar{x}_{1},1)$ mamy nierówność przeciwną. Oczywiście wystarczy rozpatrywać $x_{t}\in\left(0,\tfrac{1}{2}\right]$ , gdyż wykres iterowanej funkcji $F(x)$ jest symetryczny względem prostej $x=\tfrac{1}{2}$ oraz $F(x)\leq\tfrac{a}{4}\leq\tfrac{1}{2}$ , więc po pierwszej iteracji wyrazu $x_{0}\in\left(\tfrac{1}{2},1\right)$ dostajemy $x_{1}\in\left(0,\tfrac{1}{2}\right)$ . Wystarczy teraz wykazać, że jeśli $x_{t}<\bar{x}_{1}$ , to $x_{{t+1}}<\bar{x}_{1}$ . Nierówność $ax(1-x)<\bar{x}_{1}$ jest równoważna nierówności kwadratowej

$ax^{2}-ax+\frac{a-1}{a}>0,$

której wyróżnik $\Delta=(a-2)^{2}$ i dla $a\in(1,2]$ mamy dwa pierwiastki $\bar{x}_{1}$ oraz $x=\tfrac{1}{a}\geq\tfrac{1}{2}$ . Wobec tego nierówność ta jest spełniona dla $x\in(0,\bar{x}_{1})$ . Analogicznie wykazujemy, że $x_{t}\in\left(\bar{x}_{1},\tfrac{1}{2}\right]$ implikuje $x_{{t+1}}\in\left(\bar{x}_{1},\tfrac{1}{2}\right]$ . Mamy więc ciągi monotoniczne i ograniczone. Stosując tę samą technikę dowodzimy monotoniczności podciągów w przypadku $a\in\left(2,3\right)$ .

$Przykładowe rozwiązania równania~\eqref{log-dys}.$

Rys. 2.3. Przykładowe rozwiązania równania (2.4) dla $0<a\leq 1$ , czyli gdy $\bar{x}=0$ jest globalnie asymptotycznie stabilny (lewa strona) oraz dla $a\in(1,3)$ , czyli w przypadku stabilności stanu stacjonarnego $\bar{x}_{1}=\tfrac{a-1}{a}$ (strona prawa).

Dla omówionych wartości parametrów rozwiązania modelu dyskretnego są zbieżne, podobnie jak w przypadku modelu ciągłego, już dla tych wartości widzimy jednak zasadnicze różnice. Po pierwsze — dla małych $a\leq 1$ rozwiązania dążą do $0$ , czyli populacja wymiera. Nie należy się jednak temu dziwić, gdyż w modelu dyskretnym wartości współczynnika rozrodczości poniżej $1$ odpowiadają procesowi śmiertelności, więc w tym przypadku model dyskretny opisuje sytuację, w której populacja wymiera nawet bez konkurencji, zatem jeśli konkurencja występuje, to sytuacja może się tylko pogorszyć. Przy $a\in(1,2]$ mamy monotoniczną zbieżność do dodatniego stanu stacjonarnego, zatem zachowanie jest analogiczne jak w modelu ciągłym, ale gdy $a$ przekracza $2$ , pojawiają się gasnące oscylacje — takie zachowanie nie jest możliwe w modelu ciągłym. Co się dzieje, gdy $a$ przekracza $3$ ?

Dla $a>3$ oba rozwiązania stacjonarne stają się niestabilne. Możemy jednak sprawdzić, że po przekroczeniu tej wartości rozwiązania zyskują charakter stabilnego cyklu granicznego, najpierw o okresie $2$ , potem o okresie $4$ itd. Cykl o okresie $2$ znajdujemy stosunkowo łatwo, jako rozwiązanie równania

$F^{2}(x)=x,\quad\text{gdzie}\quad F^{2}(x)=F(F(x))=a^{2}x(1-x)(1-ax+ax^{2})$

różne od znanych już stanów stacjonarnych. Szukamy więc dodatnich miejsc zerowych wielomianu

$W(x)=a^{3}x^{3}-2a^{3}x^{2}+a^{2}(a+1)x+1-a^{2}.$

Dzielimy $W$ przez $x-\bar{x}_{1}$ i obliczamy

$c_{{1,2}}=\frac{a+1\pm\sqrt{a^{2}-2a-3}}{2a},\quad C=\{ c_{1},c_{2}\}.$

Zbiór $C$ jest poszukiwanym cyklem o okresie $2$ . Jego lokalną stabilność badamy znanymi już metodami. Dla $a\in(3,1+\sqrt{6})$ cykl jest stabilny. Gdy $a$ przekracza $1+\sqrt{6}$ następuje destabilizacja i pojawia się stabilny cykl o okresie $4$ . Taką zmianę dynamiki rozwiązań nazywamy bifurkacją podwojenia okresu (albo bifurkacją rozwidleniową).

$Przykład cyklicznych rozwiązań równania~\eqref{log-dys} o różnych długościach cyklu.$

Rys. 2.4. Przykład cyklicznych rozwiązań równania (2.4) o różnych długościach cyklu.

Wraz z rosnącym $a$ wyczerpują się okresy postaci $2^{n}$ . Przy przekroczeniu wartości krytycznej $a^{*}\approx 3,569$ pojawiają się rozwiązania o innych okresach. Okresy te są związane z porządkiem Szarkowskiego. Na końcu bifurkuje orbita okresowa o okresie $3$ , która powszechnie kojarzona jest z chaosem. Omawianie pojęcia chaosu wykracza poza ramy tego wykładu. Można jednak udowodnić, że dla $a=4$ rozwiązania równania (2.4) są chaotyczne. Cykl bifurkacji występujący w tym modelu najlepiej obrazuje diagram bifurkacyjny, zwany w tym przypadku drzewem Feigenbauma, por. rys. 2.5.

$Diagram bifurkacyjny modelu~\eqref{log-dys}, zwany drzewem Feigenbauma ({\em Źródło}: wikipedia).$

Rys. 2.5. Diagram bifurkacyjny modelu (2.4), zwany drzewem Feigenbauma (Źródło: wikipedia).

Drzewo Feigenbauma ma pewne własności obiektów zwanych fraktalami, w szczególności charakteryzuje się samopowtarzalnością — jeśli przybliżymy wycinek wykresu, to zobaczymy fragment łudząco podobny do całości. Podkreślić należy, że tego typu dynamikę generują wszystkie funkcje o własnościach podobnych do funkcji logistycznej, czyli $F:[0,1]\to[0,1]$ unimodalne.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Modelowanie matematyczne w biologii i medycynie

Zagadnienia

2. Modelowanie pojedynczej populacji II

2.1. Równanie logistyczne — model Verhulsta

2.2. Dyskretne równanie logistyczne