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 rn>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 N2, 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

gdzie b oznacza współczynnik konkurencji wewnątrzgatunkowej i odzwierciedla ,,szkodliwość” tej konkurencji dla populacji 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 rn=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)

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

i przy b=rK równanie (2.1) i równanie logistyczne

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 1N⁢1-N/K na ułamki proste,

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 C1, co implikuje istnienie i jednoznaczność rozwiązań. Łatwo też sprawdzić, że dla N0≥0 mamy N⁢t≥0 dla t≥0. Jedną z metod tej wykazania tej własności stanowi zapisanie równania (2.2) w równoważnej postaci całkowej

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

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

Zbadamy teraz zależność N,N˙, która determinuje monotoniczność rozwiązań. Przyjrzyjmy się wykresowi funkcji f⁢N=r⁢N⁢1-NK, por. rys. 2.1. Jest to parabola o miejscach zerowych N=0 i N=K skierowana w dół. Wynika stąd, że N≡0 i N¯≡K są rozwiązaniami stacjonarnymi. Jeśli N⁢t∈0,K, to N˙⁢t>0, czyli rozwiązania dla warunku początkowego N0<K są rosnące i ograniczone, ponieważ z jednoznaczności rozwiązań wynika, że nie mogą przeciąć stacjonarnego rozwiązania N¯. Zatem jeśli 0<N0<K, to

Analogicznie jeśli N0>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 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

Widzimy, że jeśli rozwiązanie w pewnej chwili tp przyjmuje wartość N⁢tp=K2, to zmienia się charakter przebiegu — funkcja N⁢t ma punkt przegięcia. Analizując znaki czynników N˙⁢t i 1-2⁢N⁢tK możemy wydzielić kilka obszarów, por. rys. 2.2

• N0∈0,K2, to N⁢t rośnie do K przy t→∞, więc istnieje tp, takie że N¨⁢tp=0 i N⁢t jest wypukła dla t<tp (szybki wzrost, porównywalny z wykładniczym, dla małych zagęszczeń), a wklęsła dla t>tp (następuje spowolnienie wzrostu dla większych zagęszczeń);

• N0∈K2,K, wtedy funkcja N⁢t rośnie, ale wzrost jest powolny, funkcja N jest wklęsła;

• N0>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.

Widzimy więc, że dla N0∈0,K2 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.

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ą N˙ ilorazem różnicowym

założymy (analogicznie jak w przypadku modelu Malthusa), że Δ⁢t jest równa hipotetycznej jednostce czasu Δ⁢t=1 i zastosujemy notację dyskretną N⁢t=Nt, t∈N. Otrzymujemy

gdzie a=1+r i K1=1+rr⁢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 xt=NtK1. Dzieląc równanie (2.3) stronami przez K1 dostajemy

czyli model jednoparametrowy. Zauważmy, że model ten ma sens, o ile kolejne wyrazy ciągu xtt=0∞ pozostają w przedziale 0,1, ponieważ reprezentują one zagęszczenia populacji w proporcji do pojemności siedliska K1. Wyrazy te obliczamy jako kolejne iteracje funkcji F⁢x=a⁢x⁢1-x, która jest funkcją kwadratową z wartością maksymalną F⁢12=a4. Zatem aby wszystkie wyrazy pozostawały w przedziale 0,1 musi zachodzić 0<a≤4. Będziemy więc badać dynamikę ciągu xtt=0∞ dla a∈0,4 i warunku początkowego x0∈0,1. Oczywiście jeśli x0=0 albo x0=1, to xt=0 dla t≥1, więc interesujące są wartości x0∈0,1.

Jeśli ciąg xtt=0∞ ma granicę, to równa się ona któremuś ze stanów stacjonarnych. Wyznaczmy więc te stany. Niech x¯ oznacza stan stacjonarny. Wtedy

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 d⁢F⁢xx¯. W naszym przypadku d⁢F⁢x=a⁢1-2⁢x. Jeśli a∈0,1, to stan stacjonarny jest wyznaczony jednoznacznie x¯=0 i d⁢F⁢0=a, przy czym nierówność 0<a<1 implikuje lokalną stabilność. Co więcej widzimy, że xt+1<xt, ponieważ xt+1=F⁢xt<xt dla dowolnego xt∈0,1 przy założeniu a≤1. Mamy więc ciąg malejący i ograniczony, czyli zbieżny, zatem stan x¯=0 jest globalnie stabilny w zbiorze 0,1.

Jeśli a>1, to istnieją dwa stany stacjonarne, przy czym x¯=0 staje się niestabilny, gdyż teraz d⁢F⁢0=a>1. Dla dodatniego stanu stacjonarnego x¯1=a-1a mamy d⁢F⁢x¯1=2-a oraz

Prostymi metodami, jak w przypadku zerowego stanu stacjonarnego pokazujemy, że dla a∈1,2 ciąg xtt=0∞ jest monotoniczny od pewnego miejsca, natomiast dla a∈2,3 możemy wyróżnić dwa podciągi monotoniczne od pewnego miejsca. Faktycznie, jeśli a∈1,2, to dla xt∈0,x¯1 zachodzi F⁢xt>xt, natomiast dla xt∈x¯1,1 mamy nierówność przeciwną. Oczywiście wystarczy rozpatrywać xt∈0,12, gdyż wykres iterowanej funkcji F⁢x jest symetryczny względem prostej x=12 oraz F⁢x≤a4≤12, więc po pierwszej iteracji wyrazu x0∈12,1 dostajemy x1∈0,12. Wystarczy teraz wykazać, że jeśli xt<x¯1, to xt+1<x¯1. Nierówność a⁢x⁢1-x<x¯1 jest równoważna nierówności kwadratowej

której wyróżnik Δ=a-22 i dla a∈1,2 mamy dwa pierwiastki x¯1 oraz x=1a≥12. Wobec tego nierówność ta jest spełniona dla x∈0,x¯1. Analogicznie wykazujemy, że xt∈x¯1,12 implikuje xt+1∈x¯1,12. Mamy więc ciągi monotoniczne i ograniczone. Stosując tę samą technikę dowodzimy monotoniczności podciągów w przypadku a∈2,3.

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≤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∈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

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

Dzielimy W przez x-x¯1 i obliczamy

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

Wraz z rosnącym a wyczerpują się okresy postaci 2n. Przy przekroczeniu wartości krytycznej a*≈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.

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→0,1 unimodalne.