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.

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

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

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.

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.

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

gdzie Xn opisuje liczebność populacji (czy inną wielkość biofizyczną) w chwili n, a G odzwierciedla prawa rządzące dynamiką tej populacji. W ogólnym przypadku Xn∈Rn. 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 Xnn=1∞. 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 Yn oznacza wektor kolejnych liczebności Yn=Xn,Xn-1T. Wtedy

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 Xn=λn⁢X0, gdzie λ jest wartością własną macierzy przekształcenia G. Podstawiając taką postać rozwiązania do zależności rekurencyjnej (4.1) otrzymujemy

zatem λ spełnia równanie

czyli równanie charakterystyczne dla macierzy 1110. Rozwiązując równanie (4.3) dostajemy

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

gdzie a i b musimy wyznaczyć korzystając z warunku początkowego F1=F2=1. Dostajemy zatem układ równań

skąd obliczamy a=55=-b i ostatecznie ogólny wyraz ciągu Fibonacciego

Zauważmy, że ciąg Fnn=1∞ jest rosnący i Fn→∞ przy n→∞. 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 F1, F2>0. Jeśli wzór (4.5) przepiszemy w postaci

to możemy zauważyć, że wraz z rosnącym n dynamika ciągu Fnn=1∞ jest coraz bliższa dynamice ciągu 55⁢1+52nn=1∞ determinowanej przez dominującą wartość własną λ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 F1, F2>0.

Widzimy więc, że dynamika ciągu Fibonacciego jest raczej ,,uboga”. Spójrzmy jednak na ciąg FnFn-1. Zdefiniujmy

Możemy sprawdzić, że ciąg xnn=1∞ oscyluje i dzieli się na dwa podciągi x2⁢kk=1∞, x2⁢k+1k=1∞, monotoniczne, zbieżne do tej samej liczby 1+52, z czego wynika, że cały ciąg zbiega do tej liczby.

Liczba ta znana jest od starożytności, oznaczamy ją Φ 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.

W ogólnym przypadku możemy wyodrębnić więcej grup wiekowych. Niech populacja P będzie podzielona na k grup wiekowych. W tej sytuacji liczebność całej populacji jest opisana za pomocą wektora Nt=Nt0,Nt2,…,NtkT, gdzie Nti 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 Nt określić Nt+1. Jak zwykle musimy zacząć od modelu heurystycznego, czyli omówienia procesów, które są według nas istotne przy opisie dynamiki populacji 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:

gdzie γi odzwierciedla przeżywalność osobników z grupy wiekowej i (zatem 1-γi to śmiertelność), a αi — rozrodczość danej grupy wiekowej. W związku z tym formalnie model Lesliego można przedstawić jako trójkę N0,α,η, gdzie wektory N0=N00,…,N0kT, α=α0,…,αk, η=η0,…,ηk-1 oznaczają odpowiednio liczebności początkowe, współczynniki rozrodczości oraz śmiertelności, ηi=1-γi.

Zapisując ten model w postaci (4.2) otrzymujemy

gdzie

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

Dla dowolnych t, p∈N zachodzi też Nt+p=Mp⁢Nt.

W przypadku macierzy Lesliego liniowa niezależność wierszy i kolumn pozwala zapisać M w postaci

gdzie Λ to diagonalna macierz wartości własnych, a V macierz kolumnowych wektorów własnych V0,…,Vk. Wynika stąd, że

Oznaczmy W=V-1⁢N0, W=w0,…,wk. Wtedy

Własności rozwiązań zależą zarówno od postaci macierzy M jak i początkowego rozkładu wieku. Twierdzenie Frobeniusa – Perrona gwarantuje istnienie dominującej rzeczywistej wartości własnej λ0≥λi, i=1⁢…,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 M spełnia pewne dodatkowe założenia, np. gdy wskaźniki i, dla których αi>0 nie mają większego wspólnego dzielnika niż 1, to λ0>λi, i=1,…,k i odpowiadający jej wektor własny M⁢V0=λ0⁢V0 ma współrzędne nieujemne i dla większości warunków początkowych asymptotyczne zachowanie rozwiązań jest determinowane przez a⁢V⁢λ0t, gdzie a=const jest stałą zależną od warunku początkowego, ponieważ ze wzoru (4.7) dla dominującej wartości własnej λ0 wynika rozkład asymptotyczny

Wobec tego dla większości warunków początkowych populacja asymptotycznie osiąga rozkład wieku V0 przy współczynniku rozrodczości λ0, gdyż

i jeśli

1. λ0>1, to liczebności wszystkich grup wiekowych rosną nieograniczenie, mówimy że populacja jest rozwojowa;

2. λ0<1, to populacja jest wymierająca;

3. λ0=1, to liczebność całej populacji i poszczególnych grup wiekowych stabilizuje się — populacja jest stacjonarna, por rys. 4.5.

Wektor V0 odzwierciedla rozkład liczebności populacji w stanie równowagi, przy czym jeśli V0=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 Nt=∑i=0kNti, to Nt∼λ0t⁢N0.

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 λ1,…,λ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

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

co łatwo wykazać indukcyjnie. Zauważmy, że jeśli m⁢s=1 i s<1, to N2⁢t=N0 oraz N2⁢t-1=m⁢N02s⁢N01, czyli mamy zachowanie cykliczne, por. rys. 4.6.