Równanie logistyczne

Trochę historii

W latach siedemdziesiątych XX wieku, na Uniwersytecie w Oxford, australijski uczony Robert May zajmował się teoretycznymi aspektami dynamiki populacyjnej. Swoje prace podsumował w artykule, który ukazał się w Nature pod prowokującym tytułem „Proste modele matematyczne z bardzo skomplikowaną dynamiką” x[may76]_. Artykuł ten po latach stał się jedną z najczęściej cytowanych prac z teoretycznej ekologii. Co wzbudziło tak wielkie zainteresowanie w tej pracy?

May zajmował się zastosowaniem matematyki w ilościowym opisie zjawisk ekologicznych. Klasycznym zadaniem w tej dziedzinie jest obliczenie wielkości populacji pewnego gatunku w czasie, znając jej stan liczebny w chwili początkowej. Najprostszym, z punktu widzenia modelowania matematycznego, rodzajem ekosystemów wydawały się takie w których życie jednego pokolenia populacji trwa jeden sezon. Dobrym przykładem jest populacja owadów, które w ciągu jednego sezonu przechodzą pełną metamorfozę np. motyle. Czas jest w naturalny sposób podzielony na dyskretne okresy, odpowiadające cyklom życia populacji. Równania opisujące taki ekosystem mają więc formę dyskretnych układów iterowanych w których bieżąca liczebność osobników danego gatunku w ekosystemie jest funkcją liczebności w poprzednim okresie.

Robert May zajmował się między innymi właśnie taką dynamiką. Badając układy iteracyjne, uprościł ekosystem do jednego gatunku w którym populacja była funkcją kwadratową populacji w roku poprzednim. Skąd taki model? Najprostszym równaniem dyskretnym opisującym ewolucję populacji jest model liniowy:

(1)\[N_{i+1} = \alpha \; N_{i},\]

gdzie \(N_i\) to liczebność w i-tym sezonie. Łatwo się przekonać, że takie równanie może prowadzić to trzech scenariuszy. Jeżeli \(\alpha>1\) to populacja będzie nieograniczenie rosnąć, jeżeli \(\alpha<1\) to zaniknie oraz dla \(\alpha=1\) ewolucja nie będzie zmieniać stanu liczebnego populacji. Najprostszym rozwinięciem tego modelu jest wprowadzenie zależności stałej \(\alpha\) od wielkości populacji. Wyobraźmy sobie populacje szkodników w zamkniętym ekosystemie. Szkodniki zjadają zboże, którego jest dokładnie taka sama ilość co roku. Jeżeli owadów jest mało w porównaniu do ilości pożywienia, to mogą rozmnażać się z pełną siła rozrodczą - na przykład w następnym sezonie będzie ich cztery razy więcej niż w poprzednim. Jednak w miarę wzrostu liczebności szkodników, pożywienia nie będzie wystarczać i siła rozrodcza będzie maleć. W krytycznym przypadku można sobie wyobrazić ze owady zjadają latem całe zboże po czym wszystkie osobniki umierają z głodu przed osiągnięciem zdolności rozrodczej. Załóżmy więc, że nasza stała rozrodu będzie liniową funkcją populacji:

(2)\[\alpha = \alpha( N_{i} ) = A - B N_{i},\]

gdzie \(A\) to stała wzrostu populacji w warunkach dostatku pożywienia a \(B\) jest stałą, która określa jak szybko brak pożywienia będzie zmniejszał siłę rozrodczą. W szczególności jeśli \(N_i=A/B\) to pożywienia jest na tyle mało, że żaden osobnik nie przeżywa sezonu żerowania.

Równanie (1) ze stałą (2), można przeskalować do postaci matematycznie równoważnej, zależnej tylko od jednego parametru. Równanie takie obecnie jest znane pod nazwą odwzorowania logistycznego:

(3)\[x_{i+1} = a x_{i} (1 - x_{i}),\]

gdzie \(a<=4\) jest pewną dodatnią stałą a \(x_i\in(0,1)\) jest proporcjonalne do liczebności populacji w i-tym sezonie.

Informacja

Jeśli populacja ma liczebność równą jeden, to nie dożywa do następnego pokoleniu. Tak samo było by w przypadku gdy jest większa od jednego, dlatego wystarczy się ograniczyć do \(x_i\in(0,1)\). Z tego samego powodu nie rozważamy parametru \(a>4\) - bowiem \(a<=4\) odwzorowanie logistyczne przeprowadza zawsze odcinek (0,1) w odcinek (0,1).

Mogło by się wydawać, że tak prosty model będzie dawał proste wyniki. Spróbujmy sami!

Rozważmy model (3) dla parametru \(a=0.5\), startując z liczebności \(x=0.45\). Kolejne wartości populacji można otrzymać stosując przekształcenie kwadratowe (3) do wartości z poprzedniego sezonu, na przykład za pomocą poniższego programu:

Wykonując ten przykład otrzymujemy kolejne wartości populacji, które wraz z upływem czasu dążą do zera. Eksperymentując z powyższym kodem łatwo też jest się przekonać, że niezależnie od wartości z której startujemy, zawsze populacja ginie.

Możemy sobie też ułatwić zadanie, wykorzystując w Sage narzędzie do szybkiego prototypowania elementów interaktywnych - dekorator @interact. Ponadto, zamiast wypisywać wartości liczbowe przedstawmy je graficzne rysując wykres liczebności populacji od czasu.

W powyższym kodzie, elementy slider pozwalają nam na wykonanie funkcji myf dla wybranych interaktywnie wartości \(x\) i \(a\).

Zwiększmy teraz parametr \(a\) do dowolnej wartości z przedziału \(a\in(1,3)\). Okazuje się, że wtedy ciąg \(x_i\) dąży do pewnej wielkości - tym razem jednak nie jest to zero. Interpretując w kategoriach ekologii, możemy powiedzieć, że wielkość populacji ustala się na pewnym poziomie, który nie zmienia się z sezonu na sezon. Podobnie jak poprzednim razem, ta wartość graniczna nie zależy od punktu startowego. Czyli niezależnie od tego czy populacja wystartuje bardzo małą liczebnością czy dużą, po kilku pokoleniach i tak będzie taka sama. W takim przypadku mamy efekt dążenia ekosystemu do stabilizacji - populacja dostosowuje swoją liczebność do możliwości wyżywienia się.

Taki efekt był oczekiwany przez badaczy i równanie logistyczne (3) nie przyciągnęło by szczególnej uwagi gdyby nie pewna niespodzianka. Okazało się bowiem, że dla pewnych wartości parametru \(a\) model nie zachowuje się w przewidywalny sposób. Pojawiają się nie tylko stany okresowe, ale i stany w których populacja z roku na rok zmienia się w chaotyczny sposób i występuje czułość na warunki początkowe - wszystkie cechy, które są charakterystyczne dla chaosu deterministycznego.

Zbadajmy to! Na początek ustalmy wartość parametru na \(a = 3.2\) i przyjrzyjmy się ewolucji. Zaskoczeniem może być fakt, że tym razem populacja nie osiąga jednej wartości, ale dwie, które występują kolejno po sobie co drugi sezon. Przyjrzyjmy się bliżej temu zjawisku. Po pierwsze jeżeli ciąg kolejnych wartości \(x_i\) dąży do pewnej granicy, to możemy napisać dokładny warunek na jej wartość \(x_g\). Musi bowiem zachodzić \(x_g=f(x_g)\). Jeżeli taki punkt istnieje dla pewnej funkcji \(f\), to mówimy, że jest to punkt stały odwzorowania. Możemy więc dokładnie wyznaczyć wartość punktów stałych odwzorowania logistycznego w zależności od parametru \(a\). Prosty rachunek pokazuje, że mamy dwa rozwiązania: \(x_g = 0\) oraz \(x_g=1-\frac{1}{a}\). O ile \(x_g = 0\) jest punktem stałym dla dowolnej wartości parametru, to pamiętając, że sens mają tylko wartości \(x_i\in(0,1)\), drugi punkt stały istnieje dla wartości \(a\in(1,4)\). Możemy narysować więc wykres punktów stałych od parametru:

Jeżeli mamy równanie zależne od parametru i ilość rozwiązań zmienia się wraz z tymże parametrem to mówimy, że następuje bifurkacja. W punkcie \(a=1\) następuje właśnie bifurkacja i układ zamiast jednego rozwiązania ma dwa. Jednak zauważmy jeszcze jedno ciekawe zjawisko. Z dowolnego warunku początkowego dla \(a<1\) zawsze otrzymywaliśmy malejący ciąg populacji, który wydawał się być przyciągany do jedynego w tym obszarze punktu stałego - do zera. Taki punkt, do którego układ jest przyciągany zwany jest też atraktorem układu. Dla \(a>1\) mamy dwa punkty stałe. Okazuje się, że w tym obszarze startując z dowolnego punktu z wyjątkiem \(x=0\) zawsze będziemy dążyć do drugiego rozwiązania, który jest atraktorem! Oznacza to, że jeżeli rozwiązanie \(x=0\) zaburzymy dowolnie małą liczbą np. \(x=0.0001\) to i tak po kilkunastu iteracjach populacja będzie dążyła do \(x_g=1-\frac{1}{a}\) (Poeksperymentujmy!). Stabilny dla \(a<1\) punkt stały \(x=0\) staje się niestabilny dla \(a>1\).

Wróćmy więc do naszej sytuacji, w której mamy \(a = 3.2\). Według poprzednich wyliczeń dalej powinniśmy mieć punkt stały \(x_g=1-\frac{1}{a}\)! I mamy, sprawdźmy:

Dodajmy jednak do wartości początkowej pewną małą liczbę np. niech x=x+1e-6. Zobaczmy co się stanie? Okazuje się, że we wcześniejszym punkcie (jak się okaże \(a=3\)) nastąpiła kolejna bifurkacja w wyniku której rozwiązanie \(x_g=1-\frac{1}{a}\) utraciło stabilność na rzecz oscylacji. Ponieważ oscylacje te są w pomiędzy dwoma wartościami, to mówimy, że dla \(a=3.2\) układ ma punkt okresowy z okresem 2. Właściwie to możemy tylko przypuszczać, że tak jest bo wynika to tylko z zabaw podczas których liczba iteracji była skończona. Możemy jednak w tym przypadku pokazać to dokładnie. Jeżeli populacja co drugi sezon przechodzi w tą samą to możemy rozważyć odwzorowanie \(g(x)=f(f(x))\), które przeprowadza układ o dwa sezony do przodu. W taki przypadku powinniśmy punkt stały dla \(g\) odpowiada punktowi okresowemu o okresie 2 dla \(f\). Zastosujmy tą chytrą sztuczkę, tym razem z pomocą Sage:

Dobrze, że możemy wyręczyć się systemem algebry komputerowej, bo niestety równanie \(f(f(x))=x\) jest równaniem czwartego stopnia! Sage na szczęście „potrafi” rozwiązywać analitycznie równania czwartego stopnia i otrzymujemy rozwiązania. Od razu widzimy wśród pierwiastków punkty stałe odwzorowania \(f\), co jest zrozumiałe, bo jeśli zachodzi \(f(x)=x\) to tym bardziej \(f(f(x))=x\). Narysujmy zatem nasz wynik.

Wykres ten, zwany diagramem bifurkacyjnym, nie jest do końca kompletny - skoro pojawiły się dwie bifurkacje to nie ma powodu, żeby zakładać, że więcej się nie pojawi! W dalszej analizie pojawia się jednak zasadniczy problem. Otóż nie możemy badać analitycznie punktów stałych dalszych złożeń odwzorowania \(f(f(f(x)))=x\), bo w poprzednim przypadku wyczerpaliśmy możliwość dokładnego znajdywania miejsc zerowych wielomianów. Zgodnie z Teoria Galois wzory analityczne na pierwiastki wielomianu kończą się w przypadku ogólnym na stopniu cztery. Oczywiście można zastosować metody przybliżone, lub metodę graficzną. Jednak okazuje się, że całkiem niezłym sposobem na poznanie struktury cykli układu jest po prostu jego symulacja na tyle długa by układ zdążył dojść wystarczająco blisko do atraktora. Zanim użyjemy tego sposobu, zapoznajmy się z metodą graficzną - jak mawiano, ilustracja jest warta tysiąca słów. Zapraszamy do lektury części II!