Jak już wspomnieliśmy, w modelu Lotki – Volterry zostały pominięte różne istotne czynniki, gdyż w zamyśle Volterry taki model powinien być możliwie najprostszy. W szczególności wewnętrzna dynamika populacji ofiar została zbudowana w oparciu o model Malthusa zakładający nieograniczoność pojemności środowiska dla tego gatunku, co przy niewielkiej liczebności drapieżników prowadzi do gwałtownego wzrostu populacji ofiar i przyczynia się do okresowości rozwiązań. Załóżmy więc, że wewnętrzna dynamika gatunku Po rządzi się innymi prawami, w szczególności liczebność tej populacji jest ograniczona przez pojemność środowiska. Najprostszym znanym nam modelem opisującym taką dynamikę jest równanie logistyczne, zbudujmy zatem model drapieżnik – ofiara w oparciu o to równanie. Wobec tego, przy oznaczeniach wprowadzonych w wyjściowym modelu Lotki – Volterry, model z ograniczoną pojemnością środowiska dla ofiar możemy opisać następującym układem równań

gdzie K oznacza pojemność środowiska, analogicznie jak w modelu logistycznym. Dynamika rozwiązań układu (7.1) zależy w istotny sposób od tego parametru. Zauważmy, że postawowe własności takie jak istnienie, jednoznaczność, nieujemność i przedłużalność rozwiązań możemy zbadać w analogiczny sposób jak dla układu (6.1). Wnioskujemy zatem, że dla dowolnego nieujemnego warunku początkowego V0,P0 istnieje jednoznaczne, nieujemne rozwiązanie układu (7.1) określone dla wszystkich t≥0. Dość prosto możemy wykazać też ograniczoność rozwiązań. Sformułujmy odpowiednie stwierdzenie

Zajmiemy się teraz analizą portretu fazowego. Izokliny zerowe są następującymi prostymi

• dla zmiennej V: V=0 lub P=ra⁢1-VK;

• dla zmiennej P: P=0 lub V=sa⁢b.

W zależności od parametru K — nietrywialne izokliny mogą się przeciąć w analizowanej przestrzeni fazowej, gdy K jest relatywnie duże, albo nie, gdy K jest małe. Zatem liczba stanów stacjonarnych zależy od K: zawsze są 2 stany 0,0 i K,0, natomiast dla K>sa⁢b istnieje dodatni stan stacjonarny V¯,P¯, gdzie V¯=sa⁢b i P¯=ra⁢1-sa⁢b⁢K. Wyznaczmy wartości własne dla poszczególnych stanów stacjonarnych. Macierz Jacobiego układu (7.1) ma postać

Stąd wyznaczamy wartości własne

• dla 0,0 mamy λ1=r>0, λ2=-s<0;

• dla K,0: λ1=-r<0, λ2=a⁢b⁢K-s;

• dla V¯,P¯ wartości własne są rozwiązaniem równania charakterystycznego

  λ2+r⁢V¯K⁢λ+a2⁢b⁢V¯⁢P¯=0,

  więc w zależności od parametrów albo λ1 i λ2 są rzeczywiste ujemne, albo są zespolone i mają ujemne części rzeczywiste.

Wobec tego stan stacjonarny 0,0 jest siodłem bez względu na parametry układu, z kolei K,0 jest albo węzłem stabilnym, o ile K<sa⁢b, czyli dodatni stan stacjonarny nie istnieje, albo siodłem, gdy K>sa⁢b, a wtedy istnieje V¯,P¯, który jest albo węzłem, albo ogniskiem stabilnym. Typ punktu krytycznego V¯,P¯ zależy także od K, jeśli K jest duże, tak że 4⁢a2⁢b2⁢K2-4⁢a⁢s⁢b⁢K-r⁢s>0, to obserwujemy ognisko i gasnące oscylacje, w przeciwnym przypadku — węzeł i zbieżność do stanu stacjonarnego jest od pewnego momentu monotoniczna.

Przystępując do narysowania portretu fazowego zauważmy, że trywialne izokliny zerowe stanowią trajektorie układu. Rzeczywiście: jeśli 0,P0, P0>0, jest warunkiem początkowym, to oczywiście V˙=0, czyli zmienna V nie zmienia się, V⁢t=0 dla t>0, natomiast dla zmiennej P mamy równanie P˙=-s⁢P, czyli P⁢t=P0⁢exp⁡-s⁢t dla t>0. Podobnie dla warunku początkowego V0,0, V0>0, zachodzi P˙=0, czyli P⁢t=0 dla t>0 oraz V˙=r⁢V⁢1-VK, zatem V⁢t spełnia równanie logistyczne, co zauważyliśmy już wcześniej w dowodzie stwierdzenia 7.1.

Jeśli izokliny nie przecinają się wewnątrz przestrzeni fazowej, to przestrzeń tę możemy podzielić na trzy obszary, por. rys. 7.1 i rys. 7.2: A pod izokliną dla V, B nad tą izokliną i na lewo od izokliny dla P oraz C — na prawo od izokliny dla P. Zauważmy, że rozwiązania dla warunku początkowego z C muszą przejść do B. Gdyby tak się nie stało, to mielibyśmy obie zmienne monotoniczne, przy czym zmienna V pozostawałaby malejąca i ograniczona, a zmienna P rosłaby nieograniczenie. W takim przypadku z pierwszego równania asymptotycznie dostajemy V˙∼c1-c2⁢P, gdzie c1 i c2 są dodatnimi stałymi, czyli V˙→-∞, co przeczy ograniczoności V. Zatem rozwiązanie przechodzi do B — tu obie zmienne są malejące, zatem albo pozostają w B i wobec tego mają granicę, czyli zbiegają do stanu stacjonarnego K,0 na brzegu obszaru, albo przechodzą do A. W A obie zmienne także są monotoniczne, tyle że V rośnie. Przebieg pola wektorowego pokazuje, że rozwiązania nie mogą wyjść z obszaru A, więc także są zbieżne i zbiegają do K,0.

W drugim przypadku, gdy istnieje dodatni stan stacjonarny, przebieg pola wektorowego wskazuje, podobnie jak dla modelu Lotki – Volterry, że rozwiązania okrążają punkt krytyczny V¯,P¯ w przestrzeni fazowej, por. rys. 7.1 i rys. 7.2. Wartości własne dla tego punktu implikują, że jest on lokalnie stabilny, ale na takiej podstawie nie możemy wykluczyć istnienia orbit zamkniętych, czyli cykli granicznych. Można to próbować wykazać na różne sposoby, w tym przypadku dzięki prostocie modelu działa kilka standardowych metod. Najprostszy sposób w przypadku układów dwuwymiarowych polega na zastosowaniu kryterium Dulaca – Bendixsona. Dla danego układu o prawej stronie opisanej funkcją G⁢V,P=G1⁢V,P,G2⁢V,PT rozważamy pewną funkcję B:R+2→R klasy C1, przy czym

nie zmienia znaku w R+2. Wtedy w R+2 nie ma cykli, w szczególności nie ma cykli granicznych. Typowo stosowana funkcja B ma postać B⁢x,y=1x⁢y. W przypadku układu (7.1) dostajemy

z czego wnioskujemy globalną stabilność stanu stacjonarnego V¯,P¯, gdyż jeśli nie ma cyklu granicznego, to z twierdzenia Poincarégo – Bendixsona wszystkie rozwiązania zbiegają do stanu stacjonarnego, w tym przypadku do jedynego stanu stabilnego V¯,P¯.

Bardziej uniwersalna metoda (działająca także w przypadku modeli o większej liczbie zmiennych) polega na znalezieniu funkcjonału Lapunowa. W tym przypadku można zastosować funkcjonał Lapunowa zaproponowany dla modelu Lotki – Volterry. W celu uproszczenia obliczeń przeprowadźmy najpierw zamianę zmiennych x=V-V¯, y=P-P¯, która przesuwa stan stacjonarny do punktu 0,0. Zauważmy, że po zamianie zmiennych powstaje nowa przestrzeń fazowa D⁢F=-V¯,+∞×-P¯,+∞, a układ (7.1) przyjmuje postać

Funkcjonał Lapunowa dla układu (7.2) definiujemy jako

Należy teraz sprawdzić odpowiednie własności funkcjonału L: w szczególności mamy L⁢x,y≥0 w przestrzeni D⁢F oraz L⁢x,y=0⁢⟺⁢x=0⁢⁢i⁢⁢y=0 (czyli V=V¯ i P=P¯). Policzmy pochodną L wzdłuż trajektorii układu (7.2)

dla α1=b i α2=1. Wynika stąd globalna stabilność dodatniego stanu stacjonarnego. Aby wykazać jego globalną asymptotyczną stabilność wystarczy udowodnić, że L⁢V⁢t,P⁢t jest ściśle malejąca, co pozostawiamy jako ćwiczenie.

Zauważmy, że tak zdefiniowany układ ma stany stacjonarne asymptotycznie stabilne, przy czym dla K<sa⁢b wszystkie rozwiązania (oprócz rozwiązań na brzegu przestrzeni fazowej) zbiegają do stanu stacjonarnego K,0, natomiast dla K>sa⁢b mamy zbieżność do dodatniego stanu stacjonarnego V¯,P¯, a przy K=sa⁢b następuje bifurkacja siodło – węzeł. Taki typ bifurkacji oznacza, że zmienia się charakter punktu krytycznego z siodła na węzeł, przy czym dla bifurkacyjnej wartości K=sa⁢b parametru K ze stanu stacjonarnego K,0 bifurkuje dodatni stan stacjonarny i dla takiego K rozwiązania te sklejają się w jedno i choć nie przeprowadziliśmy analizy portretu fazowego dla tego przypadku, to doświadczenia związane z badaniem przebiegu orbit w pozostałych przypadkach pozwalają nam wnioskować, że także i dla bifurkacyjnej wartości K mamy globalną stabilność stanu stacjonarnego K,0. Zauważmy dalej, że niewielka zmiana prawej strony układu nie prowadzi w przypadku tego modelu do drastycznych zmian zachowania rozwiązań, zatem taki układ nie ma już niepożądanej własności niestabilności strukturalnej charakteryzującej układ (6.1).

Inny sposób na uzyskanie stabilności strukturalnej stanowi przyjęcie założenia, że część ofiar jest stale niedostępna dla drapieżników — w środowisku występuje pewna liczba kryjówek, w których może się schować określona liczba ofiar. Niech liczba ofiar, która może się schować przed drapieżnikiem wynosi K. Jeśli zatem w chwili początkowej t=0 liczebność całkowita ofiar V⁢t nie przekracza K, to wszystkie ofiary chowają się przed drapieżnikami i drapieżniki nie są w stanie nic upolować. W takiej sytuacji liczebność ofiar rośnie i w pewnym momencie przekracza K — dopiero wtedy możemy mówić o właściwym układzie drapieżnik – ofiara i zastosować odpowiedni model. W modelu tym zakładamy, że liczba ofiar dostępnych dla drapieżników wynosi V-K, V>K, zatem tylko takie ofiary mogą zostać upolowane i wpłynąć na rozwój populacji drapieżników. Wyjściowy model Lotki – Volterry przyjmuje przy takim założeniu postać

przy czym parametry r, a, s, b mają takie same znaczenie jak dla układu (6.1), natomiast K oznacza tu liczbę kryjówek (w odróżnieniu od typowego oznaczenia pojemności środowiska).

Lokalne istnienie i jednoznaczność rozwiązań układu (7.3) ponownie wynikają z wielomianowej postaci prawej strony układu. Ponieważ model ma sens dla V≥K, to ograniczymy nasze rozważania do przestrzeni fazowej D⁢F=K,+∞×R+. Zauważmy, że

Z niezmienniczości zbioru D⁢F wynika od razu globalne istnienie rozwiązań: dla zmiennej V dostajemy oszacowanie liniowe V˙≤r⁢V, czyli V⁢t≤V0⁢er⁢t i dalej na każdym ograniczonym przedziale 0,t* szacujemy P˙≤A⁢P, gdzie A=a⁢b⁢V0⁢er⁢t*-s, podobnie jak w przypadku układu (6.1).

Asymptotykę układu (7.3) badamy znanymi już metodami. Przeprowadzenie analizy portretu fazowego z zastosowaniem np. kryterium Dulaca – Bedixsona pozwala wysnuć następujący wniosek

Zauważmy też, że podobnie jak dla układu (7.1), dodatni stan stacjonarny w zależności od parametrów (np. liczby ofiar mogących schować się w kryjówkach) może być albo ogniskiem, albo węzłem stabilnym, por. rys. 7.4.

Podsumowując — zarówno w przypadku modelu z kryjówkami dla ofiar jak i z pojemnością środowiska dostaliśmy globalną stabilność jednego ze stanów stacjonarnych, przy czym w przypadku modelu z kryjówkami jest to dodatni stan stacjonarny, natomiast w przypadku modelu z ograniczoną pojemnością środowiska może to być albo dodatni stan stacjonarny, albo stan brzegowy opisujący wyginięcie (ekstynkcję) populacji drapieżników. Jeśli stabilny stan stacjonarny jest węzłem, to rozwiązania od pewnego momentu t¯ pozostają monotoniczne, natomiast jeśli jest ogniskiem, to obserwujemy oscylacje i — ze względu na globalną stabilność — są to oscylacje gasnące.