Korzystając z doświadczeń związanych z omówionymi już modelami oddziaływań typu drapieżnik – ofiara przejdziemy teraz do opisu oddziaływań innego typu. Jak poprzednio zakładamy, że w środowisku występują dwa gatunki P1 i P2, ale tym razem osobniki tych gatunków konkurują ze sobą o zasoby środowiska. Zwykle budując taki model przyjmujemy, że jest to konkurencja o pożywienie, ale może ona dotyczyć także dostępnej przestrzeni życiowej, a często mamy do czynienia z obydwoma rodzajami konkurencji. Ponieważ gatunki konkurują ze sobą, to nie widać powodu do wyróżniania któregokolwiek z nich, zatem opis powinien być symetryczny — jeśli zamienimy te gatunki i nazwiemy gatunek Pi (i=1, 2) gatunkiem Pj (j=2, 1) i na odwrót, to układ powinien pozostać taki sam.

Skoro mamy skorzystać z nabytych doświadczeń, zbudujemy model w oparciu o wewnętrzną dynamikę logistyczną dla każdego z gatunków Pi. Wobec tego, jeśli w środowisku występuje tylko jeden z tych gatunków, to jego zagęszczenie Ni⁢t w chwili t opisuje równanie

gdzie ri oznacza współczynnik rozrodczości netto dla gatunku Pi, zaś Ki — pojemność środowiska dla tego gatunku. Oczywiście założenie ograniczonej pojemności środowiska implikuje też występowanie konkurencji między osobnikami tego samego gatunku. Tę konkurencję nazywamy w tym kontekście konkurencją wewnątrzgatunkową lub krócej — konkurencją wewnętrzną. Oprócz tego mamy też do czynienia z konkurencją zewnątrzgatunkową, albo konkurencją zewnętrzną, którą opisujemy podobnie do konkurencji wewnętrznej. Zakładamy więc, że osobniki konkurują ze sobą w trakcie spotkań między dwoma osobnikami różnego gatunku i liczba tych spotkań jest proporcjonalna do liczebności każdego z gatunków Pi, i=1, 2. Wobec tego składnik konkurencji zewnętrznej w obu równaniach układu jest proporcjonalny do N1⁢t⁢N2⁢t.

Ostatecznie otrzymujemy dwóch układ równań

gdzie ai⁢j są współczynnikami konkurencji zewnętrznej, przy czym dla uproszczenia obliczeń związanych z układem (8.1) współczynnik ten odnosi się do liczebności populacji Pj w stosunku do jej pojemności środowiska. Zauważmy, że dzięki takiemu zapisowi możemy łatwo przeprowadzić ubezwymiarowienie układu (8.1) wprowadzając zmienne xi=NiKi, i=1,2. Przypomnijmy, że taką zamianę zmiennych wprowadziliśmy w przypadku dyskretnego równania logistycznego, co umożliwiło ograniczenie rozważań do odcinka 0,1. W nowych zmiennych układ (8.1) przyjmuje postać

przy czym xi odzwierciedla procentową eksploatację środowiska przez gatunek i.

Przeprowadzimy teraz analizę dynamiki modelu konkurujących gatunków opisanego przez układ (8.2) podobnie jak dla układu (6.1). Zaczniemy od zbadania podstawowych własności takich jak istnienie, jednoznaczność, nieujemność rozwiązań dla nieujemnych warunków początkowych.

Po zbadaniu podstawowych własności możemy przejść do przeanalizowania portretu fazowego. Izokliny zerowe układu (8.2) wyznaczamy jako

• x˙1=0⟺x1=0 lub x2=1a12⁢1-x1;

• x˙2=0⟺x2=0 lub x2=1-a21⁢x1.

Wzajemne położenie izoklin nietrywialnych zależy od wielkości współczynników ai⁢j. Mamy trzy istotnie różne przypadki generyczne

1. jeden ze współczynników jest większy, a drugi mniejszy niż 1 — założymy, że a12>1 i a21<1, drugi układ parametrów implikuje taką samą dynamikę z dokładnością do zamiany miejscami gatunków P1 i P2;

2. a12⁢,⁢a21>1;

3. a12⁢,⁢a21<1.

Oprócz wymienionych powyżej przypadków generycznych mamy też przypadki niegeneryczne, gdy co najmniej jeden ze współczynników konkurencji jest równy 1. Analizę tych przypadków pomijamy. Od strony biologicznej nieprawdopodobne wydaje się, żeby wartość jakiegokolwiek parametru utrzymywała się na stałym poziomie — zwykle występują drobne wahania, a co za tym idzie przypadki niegeneryczne są nieistotne z biologicznego punktu widzenia.

Zbadajmy najpierw stany stacjonarne i ich stabilność. Mamy zawsze trzy stany stacjonarne na brzegu przestrzeni fazowej R+2: 0,0, 1,0 i 0,1, które odzwierciedlają wymarcie co najmniej jednego z gatunków. Jeśli nietrywialne izokliny przecinają się, to mamy dodatni stan stacjonarny

Widzimy, że dodatni stan stacjonarny istnieje, jeśli albo oba współczynniki a12⁢,⁢a21 są poniżej 1, albo oba są powyżej. W celu zbadania stabilności wyznaczamy macierz Jacobiego

i obliczamy wartości własne

• dla 0,0 mamy λi=ri>0, i=1,2;

• dla 1,0: λ1=-r1<0, λ2=r2⁢1-a21 i dla 0,1 — symetrycznie;

• dla x¯1,x¯2 wartości własne są rozwiązaniami równania charakterystycznego

  λ2+λ⁢r1⁢x¯1+r2⁢x¯2+r1⁢r2⁢x¯1⁢x¯2⁢1-a12⁢a21.

Wynika stąd, że 0,0 jest zawsze węzłem niestabilnym, stabilność 1,0 zależy od parametru a21 — jeśli a21>1, to mamy węzeł stabilny, a jeśli a21<1, to siodło (analogicznie dla 0,1 w zależności od a12). Natomiast stabilność dodatniego stanu stacjonarnego zależy od iloczynu a12⁢a21 — jeśli jest on większy niż 1, to x¯1,x¯2 jest siodłem, a dla nierówności przeciwnej — stabilnym węzłem. We wszystkich tych przypadkach następuje bifurkacja siodło – węzeł dla krytycznej wartości parametrów (czyli gdy są równe 1).

Analizę zachowania rozwiązań w przestrzeni fazowej R+2 zaczniemy od przypadku 1, tj. a21<1<a12. Przy takim układzie parametrów nietrywialne izokliny dla zmiennych x1 i x2 nie przecinają się w omawianej przestrzeni fazowej — nie istnieje dodatni stan stacjonarny. Wiemy już, że 1,0 jest siodłem, a 0,1 węzłem stabilnym. Izokliny dzielą przestrzeń fazową na 3 obszary: obszar A pod izokliną dla x1, obszar B pomiędzy izoklinami oraz C — powyżej izokliny dla x2, por. rys. 8.1.

Zauważmy, że izokliny trywialne x1=0 oraz x2=0 stanowią rozwiązania układu. Rzeczywiście, jeśli x10=0, to x˙1≡0 i stąd x1⁢t=0 dla t≥0. Wtedy x˙2=r2⁢x2⁢1-x2 i x2 spełnia równanie logistyczne, zatem na osi pionowej rozwiązanie zbiega do punktu 0,1. Dokładnie takie samo rozumowanie pokazuje, że na osi poziomej rozwiązanie zbiega do 1,0. Ponieważ ten ostatni punkt jest siodłem, więc wyznaczyliśmy w ten sposób rozmaitość stabilną. Z przebiegu pola wektorowego w otoczeniu punktu 1,0 wnioskujemy, że interesujący nas fragment rozmaitości niestabilnej przebiega w obszarze B. W obszarze A obie zmienne rosną, nie mogą pozostać w tym obszarze, bo musiałyby być zbieżne do punktu o dodatnich współrzędnych na izoklinie dla zmiennej x1 — ale taki punkt musiałby wtedy być stanem stacjonarnym, a nie ma dodatniego stanu stacjonarnego w tym przypadku. Wobec tego rozwiązanie przechodzi do obszaru B. Z kolei jeśli punkt początkowy należy do C, to obie zmienne maleją i analogiczne rozumowanie pokazuje, że albo rozwiązanie zostaje w C i wtedy zbiega do stanu stacjonarnego 0,1, albo przechodzi do obszaru B. Przebieg pola wektorowego pokazuje, że rozwiązanie nie może wyjść z obszaru B, a ponieważ w tym obszarze jest także monotoniczne — x1 maleje a x2 rośnie — to zbiega do stanu stacjonarnego, czyli do 0,1. Rozmaitość niestabilna stanowi separatrysę, która oddziela rozwiązania zaczynające się w A (dokładniej, jeśli przeprowadzimy analizę przebiegu takich rozwiązań dla t<0, to zauważymy, że zbiegają one do 0,0 przy t→-∞) i w związku z tym mające taką własność, że x2 stale rośnie (pod izokliną dla zmiennej x2), natomiast x1 najpierw rośnie, osiąga wartość maksymalną przekraczając izoklinę, a następnie maleje — musi więc maleć do 0, od rozwiązań zaczynających się w C i przechodzących do B — takie rozwiązanie ma malejącą pierwszą współrzędną, natomiast druga współrzędna osiąga minimum na izoklinie, a następnie rośnie do 1. Ostatecznie wszystkie rozwiązania oprócz x20=0 zbiegają asymptotycznie do stanu stacjonarnego 0,1 i są od pewnego miejsca monotoniczne. Wnioskujemy, że w tym przypadku gatunek P1 ginie, natomiast P2 stabilizuje się na poziomie pojemności środowiska.

Wspomnieliśmy już, że jeśli układ parametrów a12 i a21 zmienia się, tj. gdy a12<1<a21 to rozwiązania zachowują się symetrycznie — dostajemy więc zbieżność wszystkich dodatnich rozwiązań do 1,0. Opisany przebieg rozwiązań wiąże się z tym, że konkurencja zewnętrzna okazuje się znacznie bardziej niebezpieczna (a12>1 lub a21>1) dla jednego z gatunków, a mało wpływa na drugi. Wtedy ten, dla którego jest niebezpieczna — ginie.

Jeśli oba parametry a12⁢,⁢a21>1 lub a12⁢,⁢a21<1, to nietrywialne izokliny przecinają się tworząc dodatni stan stacjonarny. Jeśli jednak oba parametry są powyżej 1, stan ten jest siodłem, natomiast w przeciwnym przypadku — węzłem stabilnym. W przestrzeni fazowej R+2 możemy teraz wydzielić 4 obszary i analogicznie jak w przypadku 1., analizując przebieg pola wektorowego wnioskujemy, że w przypadku 3. dodatni stan stacjonarny jest globalnie stabilny we wnętrzu przestrzeni fazowej, natomiast w przypadku 2. rozmaitość stabilna tworzy separatrysę rozdzielającą baseny przyciągania stanów 1,0 i 0,1 — w zależności od warunku początkowego jeden gatunek przeżywa, a drugi ginie.

Podsumowując przeprowadzoną analizę zauważmy, że jeśli tylko a12⁢a21>1, co oznacza, że co najmniej jeden ze współczynników a12, a21 przewyższa 1, to dynamika modelu odzwierciedla znaną zasadę ekologiczną mówiącą o konkurencyjnym wykluczaniu gatunków — jeśli nisze ekologiczne dwóch konkurujących gatunków zbytnio się pokrywają, to jeden z gatunków wypiera drugi. Tylko przy bardzo szczególnych układach parametrów możliwe jest współistnienie takich gatunków w tym samym siedlisku.

Kolejnym typem oddziaływań występujących między dwoma gatunkami P1, P2 jest symbioza. Analogicznie jak w przypadku modelu konkurujących gatunków budując model heurystyczny zakładamy, że wewnętrzną dynamikę każdego z gatunków można opisać za pomocą modelu logistycznego, ze względu na ograniczoną pojemność środowiska, a co za tym idzie — występowanie konkurencji wewnątrzgatunkowej. Oddziaływania międzygatunkowe pojawiają się oczywiście wtedy, gdy spotykają się osobniki różnych gatunków i — ponownie podobnie do poprzednich przypadków — zakłada się, że wpływ tych oddziaływań na dynamikę każdego z gatunków zależy od liczby spotkań między osobnikami gatunków P1 oraz P2.

Niech Ni⁢t oznaczają liczbę osobników (zagęszczenia) gatunków Pi, i=1,2. Układ symbiotyczny możemy opisać za pomocą następujących równań

gdzie tradycyjnie ri oznacza współczynnik rozrodczości gatunku i, pojemności środowiska dla danego gatunku oznaczamy Ki, natomiast bi⁢j opisuje siłę wpływu oddziaływań symbiotycznych dla gatunku i w symbiozie z gatunkiem j.

Zauważmy, że lokalne istnienie, jednoznaczność i nieujemność rozwiązań układu (8.3) dla nieujemnego warunku początkowego wykazujemy dość łatwo — stosując podobne argumenty jak w przypadku układu (8.1). Natomiast określoność rozwiązań dla wszystkich t>0 staje się tu własnością pożądaną, ale nie zawsze występującą. Rozpatrzmy następującą sytuację. Niech oba gatunki charakteryzują się tymi samymi parametrami, r1=r2, K1=K2 i b12=b21. Załóżmy też, że na początku liczebność obu gatunków jest jednakowa. Wobec tego, ze względu na całkowitą symetrię — liczebność tych gatunków jest jednakowa dla dowolnej chwili t>0, dla której istnieje rozwiązanie układu (8.3). Skoro N1⁢t=N2⁢t, to

i jeśli b12>1K1, to N1 stale rośnie. Co więcej, N˙1>γ⁢N12, gdzie γ=b12-1K1>0. Wiadomo natomiast, że rozwiązania równania x˙=γ⁢x2 dla dodatniego warunku początkowego x0>0 są określone w skończonym przedziale czasu 0,t¯, ponieważ dla t¯=1γ⁢x0 zachodzi limt→t¯-⁡x⁢t=+∞. Rzeczywiście, rozwiązując to równanie metodą rozdzielenia zmiennych dostajemy x⁢t=x01-γ⁢x0⁢t i dla t=t¯ mianownik rozwiązania zeruje się. Tego typu zachowania możemy się spodziewać za każdym razem, gdy w modelu występuje składnik kwadratowy z dodatnim współczynnikiem.

Przeprowadzimy teraz analizę portretu fazowego układu (8.3) w zależności od parametrów. Izokliny zerowe zadane są za pomocą następujących prostych

• dla zmiennej N1: N1=0 lub N2=1b12⁢N1K1-1;

• dla zmiennej N2: N2=0 lub N2=K2⁢1+b21⁢N1.

Widzimy więc, że niezerowe izokliny są prostymi o dodatnim współczynniku kierunkowym, zatem w zależności od parametrów mogą się przecinać lub nie w pierwszej ćwiartce R+2, czyli naszej przestrzeni fazowej. Wobec tego układ ma albo 3 stany stacjonarne 0,0, K1,0 i 0,K2, jeśli izokliny się nie przecinają, albo 4 takie stany, jeśli się przecinają i wtedy mamy dodatni stan stacjonarny

Zauważmy, że stan ten istnieje przy założeniu b12⁢b21⁢K1⁢K2<1, co w praktyce oznacza, że przy zadanych pojemnościach środowiska dla rozważanych gatunków zysk płynący z oddziaływań symbiotycznych nie może być zbyt duży, aby gatunki mogły pozostawać w pewnej równowadze, opisanej przez stan stacjonarny. W takim układzie mamy N¯i>Ki, i=1,2, czyli współrzędne dodatniego stanu stacjonarnego przewyższają pojemności środowiska.

Wyznaczmy teraz wartości własne odpowiadające poszczególnym stanom stacjonarnym. Macierz Jacobiego układu (8.3) jest równa

Podstawiając do wzoru (8.4) współrzędne poszczególnych stanów stacjonarnych otrzymujemy wartości własne

• dla 0,0 mamy λi=ri>0, i=1,2;

• dla K1,0: λ1=-r1<0, λ2=r2⁢1+b21⁢K1>0 i symetrycznie dla 0,K2: λ1=r1⁢1+b12⁢K2>0, λ2=-r2<0;

• dla N¯1,N¯2: λi, i=1,2 są rozwiązaniami równania charakterystycznego

  λ2+r1⁢N¯1K1+r2⁢N¯2K2⁢λ+r1⁢r2⁢N¯1⁢N¯2K1⁢K2⁢1-b12⁢b21⁢K1⁢K2⟹λ1,⁢λ2<0.

Wobec tego 0,0 jest węzłem niestabilnym, K1,0 i 0,K2 są siodłami, natomiast dodatni stan stacjonarny, jeśli istnieje, to jest węzłem stabilnym.

Na rysunku 8.3 widzimy przebieg pola wektorowego N˙1,N˙2 w obu przypadkach. Jeśli izokliny nie przecinają się w analizowanej przestrzeni fazowej, to łatwo możemy wywnioskować, że obie współrzędne rozwiązania są nieograniczone i albo są stale rosnące, albo na początku maleją osiągając w pewnej chwili tm minimum, a potem rosną dla t>tm. Rzeczywiście, przestrzeń fazową możemy podzielić na trzy obszary, oznaczone na portrecie fazowym A1 (ponad izokliną zerową dla N2), A2 (pod izokliną zerową dla N1) i B (pomiędzy izoklinami). W obszarze A1 zmienna N2 maleje, a N1 rośnie. Wewnątrz tego obszaru nie ma stanów stacjonarnych, a stan stacjonarny na brzegu jest siodłem, przy czym jego rozmaitość stabilna pokrywa się z osią N2, zatem rozwiązanie nie może zbiegać do tego stanu. Wnioskujemy więc, że rozwiązanie opuszcza ten obszar i przechodzi do obszaru B. Symetrycznie przebiegają orbity w obszarze A2, zatem także wchodzą do B. Z kolei w B obie zmienne są rosnące, przy czym rosną nieograniczenie — gdyby wzrost którejś z nich był ograniczony, to musiałaby być zbieżna, co oznacza, że pochodna też byłaby ograniczona, a nawet dążyłaby do 0, co jest niemożliwe ze względu na przebieg pola wektorowego.

Wróćmy teraz do problemu określoności rozwiązań dla wszystkich t>0. Zauważmy, że gdyby np. zmienna N1⁢t rosła do +∞ w skończonym czasie, to oznaczałoby, że pewnej chwili musiałaby przeciąć izoklinę zerową zmiennej N2, co znów nie zgadza się z przebiegiem pola. Zatem obie zmienne są określone dla wszystkich t>0 i limt→∞⁡Ni⁢t=+∞, i=1,2, oczywiście poza przypadkiem omówionym powyżej, gdy poszczególne współczynniki dla obu gatunków pokrywają się.

Jeśli izokliny zerowe przecinają się w pierwszej ćwiartce układu współrzędnych, to przestrzeń fazowa dzieli się na 4 obszary, A1 — ponad obydwoma izoklinami, A2 — pod nimi, B1 i B2 — pomiędzy nimi. W obszarze B1 obie zmienne rosną, natomiast w B2 — obie maleją. Zauważmy, że każdy prostokąt o jednym wierzchołku w początku układu współrzędnych i dwóch bokach wzdłuż osi oraz naprzeciwległym wierzchołku w obszarze B2 jest zbiorem niezmienniczym, gdyż pole wektorowe wchodzi do wewnątrz. Wobec tego każde rozwiązanie pozostaje ograniczone i analizując pole wektorowe dochodzimy do wniosku, że każde rozwiązanie układu (8.3) z dodatnim warunkiem początkowym zbiega do punktu krytycznego N¯1,N¯2.

Zauważmy, że symbioza zawsze działa korzystnie na liczebność gatunków, które w taki sposób na siebie oddziałują, ale tylko w przypadku drugim, gdy efekt tych oddziaływań nie jest nadmierny, gatunki mogą współistnieć w środowisku. Natomiast jeśli oddziaływania symbiotyczne przekraczają pewien próg, to prowadzi do nieograniczonego wzrostu obu populacji, co oczywiście w warunkach rzeczywistych nie może mieć miejsca i musi zakończyć się katastrofą, por. rys. 8.4.