Zagadnienia

4.1 Ciąg Fibonacciego
4.2 Macierze Lesliego

4. Modele pojedynczej populacji z uwzględnieniem wieku I

Przejdziemy teraz do opisu sytuacji, gdy populacja nie jest jednorodna — rozróżniamy osobniki w różnym wieku.

4.1. Ciąg Fibonacciego

W najprostszym przypadku w populacji wyodrębniamy dwie grupy wiekowe:

osobników niedojrzałych (w wieku przed reprodukcyjnym);
osobników dojrzałych (w wieku reprodukcyjnym).

Co ciekawe, model opisujący taką sytuację jest najstarszym znanym modelem populacyjnym i znamy go pod nazwą ciągu Fibonacciego.

W I połowie XIII wieku Leonardo z Pizy (Fibonacci, czyli ,,filus Bonacci” — syn Bonacciego) zastosował ten ciąg do opisu następującego zagadnienia populacyjnego.

,,Pewien człowiek wziął parę królików i umieścił je w miejscu otoczonym ze wszystkich stron murem. Ile par królików urodzi się z tej pary w ciągu roku, jeśli założymy, że z każdej pary po miesiącu rodzi się nowa para, która staje się płodna po upływie kolejnego miesiąca?”

Liber abaci rozdział III.

Rys. 4.1. Ilustracja zagadnienia populacyjnego opisanego przez Leonarda z Pizy.

Jako rozwiązanie tego zagadnienia Fibonacci zaproponował ciąg:

$1,\; 2,\; 3,\; 5,\; 8,\; 13,\; 21,\; 34,\; 55,\; 89,\; 144\ldots,$

znany dziś właśnie jako ciąg Fibonacciego. Ogólny wyraz tego ciągu opisujemy za pomocą formuły rekurencyjnej

$F_{{n+1}}=F_{n}+F_{{n-1}},\quad F_{1}=1,\ F_{2}=1.$

(4.1)

W Liber abaci Fibonacci pominął pierwszy wyraz ciągu — zauważmy, że powinny być dwa wyrazy równe $1$ , gdyż pierwsza para też staje się płodna dopiero po upływie miesiąca, por. rys. 4.1. Analogiczny schemat opisuje liczbę przodków pojedynczego trutnia w roju pszczół, por. rys. 4.2.

Rys. 4.2. Ilustracja przedstawiająca schemat opisujący liczbę przodków pojedynczego trutnia w roju pszczół.

Nazwijmy kolejne wyrazy ciągu Fibonacciego liczbami Fibonacciego. Okazuje się, że liczby Fibonacciego opisują różne wielkości przyrodnicze, jak np. liczbę płatków kwiatów (stokrotki stanowią sztandarowy przykład, typowo liczba płatków równa jest 34, 55 albo 89!), liczbę pędów roślin w kolejnych fazach wzrostu, czy też liczbę spiral w różnych konstrukcjach spiralnych (spirale prawoskrętne i lewoskrętne), takich jak kwiatostany słonecznika, owoc ananasa czy szyszki, por. rys. 4.3.

Zdjęcia przedstawiające wielkości przyrodnicze, które są opisywane przez liczby Fibonacciego. Od lewej: liczba płatków stokrotki; liczba pędów roślin; liczba spiral w kwiatostanie słonecznika.

Rys. 4.3. Zdjęcia przedstawiające wielkości przyrodnicze, które są opisywane przez liczby Fibonacciego. Od lewej: liczba płatków stokrotki; liczba pędów roślin; liczba spiral w kwiatostanie słonecznika.

Równanie (4.1) jest liniowym równaniem dyskretnym, przy czym możemy je traktować jako równanie z opóźnionym argumentem — typowo dyskretny układ dynamiczny bez opóźnienia zapisujemy w postaci

$X_{{n+1}}=G(X_{n}),$

(4.2)

gdzie $X_{n}$ opisuje liczebność populacji (czy inną wielkość biofizyczną) w chwili $n$ , a $G$ odzwierciedla prawa rządzące dynamiką tej populacji. W ogólnym przypadku $X_{n}\in\mathbb{R}^{n}$ . Oczywiście postać funkcji $G$ zależy od modelu heurystycznego (np. dla modelu Malthusa (1.3) funkcja $G$ jest liniowa) i w istotny sposób wpływa na przebieg rozwiązań, czyli dynamikę ciągu $\left(X_{n}\right)_{{n=1}}^{{\infty}}$ . Dyskretny układ z opóźnieniem zawsze możemy sprowadzić do układu postaci (4.2) wprowadzając pomocnicze zmienne. Przyjrzyjmy się tej procedurze na przykładzie ciągu Fibonacciego (4.1). Niech $Y_{n}$ oznacza wektor kolejnych liczebności $Y_{n}=(X_{n},X_{{n-1}})^{{\text{T}}}$ . Wtedy

$Y_{{n+1}}=\left(\begin{array}[]{cc}1&1\\ 1&0\end{array}\right)Y_{n},$

zatem jedno równanie z opóźnieniem rzędu jednej jednostki czasu zamienia się na układ dwóch równań postaci (4.2). Liczba równań w układzie wyjściowym zależy od wielkości opóźnień w układzie początkowym.

Jak już przypominaliśmy przy rozwiązywaniu równania (1.3), rozwiązań układu (4.2) dla liniowej funkcji $G$ szukamy w postaci $X_{n}=\lambda^{n}X_{0}$ , gdzie $\lambda$ jest wartością własną macierzy przekształcenia $G$ . Podstawiając taką postać rozwiązania do zależności rekurencyjnej (4.1) otrzymujemy

$\lambda^{{n+2}}X_{0}=\lambda^{{n+1}}X_{0}+\lambda^{n}X_{0},$

zatem $\lambda$ spełnia równanie

$\lambda^{2}-\lambda-1=0,$

(4.3)

czyli równanie charakterystyczne dla macierzy $\left(\begin{array}[]{cc}1&1\\ 1&0\end{array}\right)$ . Rozwiązując równanie (4.3) dostajemy

$\lambda _{1}=\frac{1+\sqrt{5}}{2},\quad\lambda _{2}=\frac{1-\sqrt{5}}{2}.$

Wobec tego wyraz ogólny ciągu Fibonacciego obliczamy jako kombinację liniową

$F_{n}=a\lambda _{1}^{n}+b\lambda _{2}^{n},$

(4.4)

gdzie $a$ i $b$ musimy wyznaczyć korzystając z warunku początkowego $F_{1}=F_{2}=1$ . Dostajemy zatem układ równań

$\begin{array}[]{lclcc}a(1+\sqrt{5})&+&b(1-\sqrt{5})&=&2,\\ a{(1+\sqrt{5})^{2}}&+&b{(1-\sqrt{5})^{2}}&=&4,\\ \end{array}$

skąd obliczamy $a=\tfrac{\sqrt{5}}{5}=-b$ i ostatecznie ogólny wyraz ciągu Fibonacciego

$F_{n}=\frac{\sqrt{5}}{5}\cdot\frac{(1+\sqrt{5})^{n}-(1-\sqrt{5})^{n}}{2^{n}}.$

(4.5)

Zauważmy, że ciąg $\left(F_{n}\right)_{{n=1}}^{{\infty}}$ jest rosnący i $F_{n}\to\infty$ przy $n\to\infty$ . Co więcej, tę samą własność mają wszystkie ciągi określone formułą rekurencyjną (4.1) bez względu na warunek początkowy, o ile tylko $F_{1}$ , $F_{2}>0$ . Jeśli wzór (4.5) przepiszemy w postaci

$F_{n}=\frac{\sqrt{5}}{5}\left(\frac{1+\sqrt{5}}{2}\right)^{n}\left(1-\left(\frac{1-\sqrt{5}}{1+\sqrt{5}}\right)^{n}\right),$

to możemy zauważyć, że wraz z rosnącym $n$ dynamika ciągu $\left(F_{n}\right)_{{n=1}}^{{\infty}}$ jest coraz bliższa dynamice ciągu $\left(\frac{\sqrt{5}}{5}\left(\frac{1+\sqrt{5}}{2}\right)^{n}\right)_{{n=1}}^{{\infty}}$ determinowanej przez dominującą wartość własną $\lambda _{1}$ . Wnioskujemy więc, że wraz z rosnącym $n$ dynamika ta coraz mniej różni się od dynamiki dyskretnego modelu Malthusa (1.3). Podobne zachowanie zaobserwujemy dla większości wartości początkowych $F_{1}$ , $F_{2}>0$ .

Widzimy więc, że dynamika ciągu Fibonacciego jest raczej ,,uboga”. Spójrzmy jednak na ciąg $\tfrac{F_{n}}{F_{{n-1}}}$ . Zdefiniujmy

$x_{n}=\frac{F_{n}}{F_{{n-1}}}=\frac{F_{{n-1}}+F_{{n-2}}}{F_{{n-1}}}=1+\frac{1}{\frac{F_{{n-1}}}{F_{{n-2}}}}=1+\frac{1}{x_{{n-1}}}.$

Możemy sprawdzić, że ciąg $(x_{n})_{{n=1}}^{{\infty}}$ oscyluje i dzieli się na dwa podciągi $(x_{{2k}})_{{k=1}}^{{\infty}}$ , $(x_{{2k+1}})_{{k=1}}^{{\infty}}$ , monotoniczne, zbieżne do tej samej liczby $\tfrac{1+\sqrt{5}}{2}$ , z czego wynika, że cały ciąg zbiega do tej liczby.

$Początkowe wyrazy ciągu $x_{n}=\frac{F_{n}}{F_{{n-1}}}$.$

Rys. 4.4. Początkowe wyrazy ciągu $x_{n}=\frac{F_{n}}{F_{{n-1}}}$ .

Liczba ta znana jest od starożytności, oznaczamy ją $\Phi$ i definiuje ona tzw. złotą proporcję lub złoty podział. Nazwa ta wiąże się z podziałem odcinka na dwie części, z których dłuższa do krótszej ma się tak jak cały odcinek do dłuższej. W starożytności mówiono także o ,,boskiej proporcji” i wykorzystywano ją w wielu różnych konstrukcjach, także przy budowie Partenonu. Również w czasach nowożytnych wielu ludzi doszukuje się liczb Fibonacciego czy złotej proporcji w różnych zjawiskach, najczęściej jednak nie ma to prawie żadnych podstaw naukowych.

4.2. Macierze Lesliego

W ogólnym przypadku możemy wyodrębnić więcej grup wiekowych. Niech populacja $\mathcal{P}$ będzie podzielona na $k$ grup wiekowych. W tej sytuacji liczebność całej populacji jest opisana za pomocą wektora $\mathbf{N}_{t}=\left(N^{0}_{t},N^{2}_{t},\ldots,N^{k}_{t}\right)^{{\text{T}}}$ , gdzie $N^{i}_{t}$ oznacza liczebność grupy wiekowej $i$ w chwili $t$ , dla $i=0$ mamy osobniki nowonarodzone, a dla $i=k$ — osobniki najstarsze, które nie mają szansy dożyć do następnej chwili. Podobnie jak w przypadku modelu Malthusa chcemy na podstawie znajomości wektora liczebności $\mathbf{N}_{t}$ określić $\mathbf{N}_{{t+1}}$ . Jak zwykle musimy zacząć od modelu heurystycznego, czyli omówienia procesów, które są według nas istotne przy opisie dynamiki populacji $\mathcal{P}$ . Ponownie podobnie jak w modelu Malthusa zakładamy, że osobniki są jednorodne w obrębie każdej grupy wiekowej, w każdej grupie wiekowej mamy do czynienia tylko z dwoma procesami — rozrodczości i starzenia się/śmiertelności, przy czym współczynniki rozrodczości i starzenia się są charakterystyczne dla danej grupy wiekowej. Analogicznie jak w przypadku modelu Malthusa w obrębie danej grupy wiekowej oba procesy przebiegają tak samo dla wszystkich osobników i są rozłożone równomiernie w czasie. Zakładamy też, że jednostka czasu jest równa jednostce zmiany wieku (starzenia się), czyli że po upływie np. jednego roku osobnik starzeje się o rok, co więcej po upływie jednej jednostki czasu przechodzi do następnej grupy wiekowej albo umiera. Dostajemy więc:

$N_{{t+1}}^{{i+1}}=\gamma _{i}N_{t}^{i},\ i\in\mathbb{N},\quad N_{{t+1}}^{0}=\sum\limits _{{i=0}}^{k}\alpha _{i}N_{t}^{i},$

gdzie $\gamma _{i}$ odzwierciedla przeżywalność osobników z grupy wiekowej $i$ (zatem $1-\gamma _{i}$ to śmiertelność), a $\alpha _{i}$ — rozrodczość danej grupy wiekowej. W związku z tym formalnie model Lesliego można przedstawić jako trójkę $(\mathbf{N}_{0},\alpha,\eta)$ , gdzie wektory $\mathbf{N}_{0}=(N_{0}^{0},\ldots,N_{0}^{k})^{T}$ , $\alpha=(\alpha _{0},\ldots,\alpha _{k})$ , $\eta=(\eta _{0},\ldots,\eta _{{k-1}})$ oznaczają odpowiednio liczebności początkowe, współczynniki rozrodczości oraz śmiertelności, $\eta _{i}=1-\gamma _{i}$ .

Zapisując ten model w postaci (4.2) otrzymujemy

$\mathbf{N}_{{t+1}}=\mathbf{M}\mathbf{N}_{t},$

(4.6)

gdzie

$\mathbf{M}=\left(\begin{array}[]{ccccc}\alpha _{0}&\alpha _{1}&\alpha _{2}&\ldots&\alpha _{k}\\ \gamma _{0}&0&0&\ldots&0\\ 0&\gamma _{1}&0&\ldots&0\\ \vdots&0&\ddots&\ddots&\vdots\\ 0&\ldots&\ldots&\gamma _{{k-1}}&0\end{array}\right),$

przy czym macierz $\mathbf{M}$ nazywamy macierzą Lesliego, a model (4.6) modelem (macierzowym) Lesliego (patrz [8, 9, 2]). Jeśli mamy zadany początkowy rozkład wieku $\mathbf{N}_{0}$ , to rozkład w dowolnej chwili $t\in\mathbb{N}$ możemy obliczyć ze wzoru

$\mathbf{N}_{t}=\mathbf{M}^{t}\mathbf{N}_{0}.$

Dla dowolnych $t$ , $p\in\mathbb{N}$ zachodzi też $\mathbf{N}_{{t+p}}=\mathbf{M}^{p}\mathbf{N}_{t}$ .

Definicja 4.1

Stanem równowagi modelu Lesliego nazwiemy parę $(\bar{\mathbf{N}},r)$ spełniającą zależność

$r\bar{\mathbf{N}}=\mathbf{M}\bar{\mathbf{N}},$

gdzie $\bar{\mathbf{N}}$ jest wektorem liczebności, który nazywamy ustaloną strukturą wieku, a $r\in\mathbb{R}$ oznacza stały współczynnik wzrostu.

W przypadku macierzy Lesliego liniowa niezależność wierszy i kolumn pozwala zapisać $\mathbf{M}$ w postaci

$\mathbf{M}=\mathbf{V}\Lambda\mathbf{V}^{{-1}},$

gdzie $\Lambda$ to diagonalna macierz wartości własnych, a $\mathbf{V}$ macierz kolumnowych wektorów własnych $V_{0},\ldots,V_{{k}}$ . Wynika stąd, że

$\mathbf{M}^{t}=(\mathbf{V}\Lambda\mathbf{V}^{{-1}})^{t}=\mathbf{V}\Lambda^{t}\mathbf{V}^{{-1}}\ \Longrightarrow\ \mathbf{N}_{t}=\mathbf{V}\Lambda^{t}\mathbf{V}^{{-1}}\mathbf{N}_{0}.$

Oznaczmy $W=\mathbf{V}^{{-1}}\mathbf{N}_{0}$ , $W=(w_{0},\ldots,w_{{k}})$ . Wtedy

$\mathbf{N}_{t}=\sum\limits _{{i=0}}^{{k}}V_{i}\lambda _{i}^{t}w_{i}.$

(4.7)

Własności rozwiązań zależą zarówno od postaci macierzy $\mathbf{M}$ jak i początkowego rozkładu wieku. Twierdzenie Frobeniusa – Perrona gwarantuje istnienie dominującej rzeczywistej wartości własnej $\lambda _{0}\geq|\lambda _{i}|$ , $i=1\ldots,k$ . Dla tej wartości własnej zarówno prawo jak i lewostronny wektor własny jest rzeczywisty i ma nieujemne współrzędne. Jeśli macierz $\mathbf{M}$ spełnia pewne dodatkowe założenia, np. gdy wskaźniki $i$ , dla których $\alpha _{i}>0$ nie mają większego wspólnego dzielnika niż $1$ , to $\lambda _{0}>|\lambda _{i}|$ , $i=1,\ldots,k$ i odpowiadający jej wektor własny $\mathbf{M}V_{0}=\lambda _{0}V_{0}$ ma współrzędne nieujemne i dla większości warunków początkowych asymptotyczne zachowanie rozwiązań jest determinowane przez $aV\lambda _{0}^{t}$ , gdzie $a=\text{const}$ jest stałą zależną od warunku początkowego, ponieważ ze wzoru (4.7) dla dominującej wartości własnej $\lambda _{0}$ wynika rozkład asymptotyczny

$\frac{1}{\lambda _{0}^{t}}\mathbf{N}_{t}=V_{0}w_{0}+\sum\limits _{{i=1}}^{{k}}V_{i}\left(\frac{\lambda _{i}}{\lambda _{0}}\right)^{t}w_{i}\to V_{0}w_{0}\ \ \text{przy}\ \ t\to\infty.$

Wobec tego dla większości warunków początkowych populacja asymptotycznie osiąga rozkład wieku $V_{0}$ przy współczynniku rozrodczości $\lambda _{0}$ , gdyż

$\mathbf{N}_{t}\approx w_{0}V_{0}\lambda _{0}^{t}\quad\text{przy}\quad t\to\infty$

i jeśli

$\lambda _{0}>1$ , to liczebności wszystkich grup wiekowych rosną nieograniczenie, mówimy że populacja jest rozwojowa;
$\lambda _{0}<1$ , to populacja jest wymierająca;
$\lambda _{0}=1$ , to liczebność całej populacji i poszczególnych grup wiekowych stabilizuje się — populacja jest stacjonarna, por rys. 4.5.

$Przykładowe rozwiązania modelu~\eqref{Leslie} w przypadku trzech grup wiekowych w zależności od dominującej wartości własnej $\lambda _{0}$. Wstawki pokazują ewolucję procentowego udziału poszczególnych grup wiekowych w całości populacji.$

Rys. 4.5. Przykładowe rozwiązania modelu (4.6) w przypadku trzech grup wiekowych w zależności od dominującej wartości własnej $\lambda _{0}$ . Wstawki pokazują ewolucję procentowego udziału poszczególnych grup wiekowych w całości populacji.

Wektor $V_{0}$ odzwierciedla rozkład liczebności populacji w stanie równowagi, przy czym jeśli $||V_{0}||=1$ , to współrzędne tego wektora odzwierciedlają procentowy wkład poszczególnych grup wiekowych w populację. W takiej sytuacji dynamika modelu Lesliego jest podobna do dynamiki modelu Malthusa, co więcej, jeśli rozpatrzymy całkowitą liczebność populacji $N_{t}=\sum\limits _{{i=0}}^{k}N_{t}^{i}$ , to $N_{t}\sim\lambda _{0}^{t}N_{0}$ .

Możliwe jest jednak także inne zachowanie rozwiązań, w szczególności cykliczne zmiany struktury wieku. Wiąże się to między innymi z dynamiką determinowaną przez początkowe rozkłady wieku należące do podprzestrzeni generowanych przez wektory własne odpowiadające wartościom własnym $\lambda _{1},\ldots,\lambda _{k}$ w sytuacji, gdy wśród tych wartości własnych nie ma już kolejnej dominującej. Oczywiście zależy to w istotny sposób od postaci macierzy Lesliego. Rozpatrzmy przykład macierzy Lesliego odzwierciedlającej sytuację, gdy rozmnażają się tylko osobniki najstarsze. Rozpatrzmy przykład populacji z wydzielonymi dwoma grupami wiekowymi — rozróżniamy osobniki dojrzałe, zdolne do reprodukcji i osobniki niedojrzałe. Mamy więc

$\mathbf{M}=\left(\begin{array}[]{cc}0&m\\ s&0\end{array}\right).$

(4.8)

Obliczając kolejne potęgi macierzy (4.6) dostajemy

$\mathbf{M}^{{2t-1}}=\left(\begin{array}[]{cc}0&m^{t}s^{{t-1}}\\ s^{t}m^{{t-1}}&0\end{array}\right)\quad\text{oraz}\quad\mathbf{M}^{{2t}}=\left(\begin{array}[]{cc}m^{t}s^{{t}}&0\\ 0&s^{t}m^{{t}}\end{array}\right),$

co łatwo wykazać indukcyjnie. Zauważmy, że jeśli $ms=1$ i $s<1$ , to $\mathbf{N}_{{2t}}=\mathbf{N}_{0}$ oraz $\mathbf{N}_{{2t-1}}=\left(\begin{array}[]{c}mN^{2}_{0}\\ sN^{1}_{0}\end{array}\right)$ , czyli mamy zachowanie cykliczne, por. rys. 4.6.