Równanie logistyczne II

Tkamy pajęczynę

Jednym z ciekawych sposobów poznania dynamiki układów dyskretnych są wykresy „pajęczynowe” (z ang. cobweb plot), znane także jako wykresy Verhulsta. Spróbujemy samodzielnie skonstruować taki wykres. Na osiach będą umieszczone wartości populacji w kolejnych iteracjach: \(x_i,x_{i+1}\), zakresy obu osi będą więc \(0..1\). Zaczynamy od narysowania prostej \(x_{i+1}=x_i\), będącej przekątną wykresu, a następnie wykresu zależności \(x_{i+1}=f(x_i)\), dla pewnego ustalonego parametru \(a\). Chcemy przedstawić trajektorię ewolucji pewnego stanu początkowego \(x_0\). Procedura rysowania składa się z czterech etapów:

  1. Znajdujemy punkt przecięcia pionowej prostej przechodzącej
    przez punkt \((x_0,0)\) z wykresem funkcji \(f\), czyli: \(x_0, f(x_0)\)
  2. Łączymy ten punkt poziomą linią z przekątną, tzn. z punktem \(f(x_0), f(x_0)\).
  3. Linią pionową łączymy powyższy punkt z wykresem funkcji \(f\), czyli z punktem \(f(x_0), f(f(x_0))\).
  4. Potwarzamy dowolną ilość razy kroki 2 i 3.

Powyższy algorytm łatwo jest wykonać nawet na kartce papieru, bez użycia komputera. Wystarczy na wykresie zawierającym przekątną oraz krzywą \(f(x)\), łączyć naprzemienne funkcję z przekątną i przekątną i z funkcją, odcinkami, odpowiednio: poziomymi i pionowymi.

Poniższa implementacja, oprócz rysowania wykresu, koloruje pierwszych pięć iteracji na niebiesko a ostatnie pięć na czerwono, co pozwala na lepsze dostrzeżenie pojawiających się cykli. Zachęcamy do eksperymentowania z poniższym kodem i manipulacji sposobem wizualizacji.

Badanie układu można rozpocząć od przyglądania się jak układ dąży do zerowego punktu stałego dla \(a<1\). W tym przypadku widać brak punktu przecięcia się paraboli z przekątną, z wyjątkiem zera. W obszarze parametru \(1<a<3\) parabola ma niezerowy punkt przecięcia się z przekątną. Zwiększanie parametru powyżej \(a=2\) powoduje, że ewolucja coraz wolniej dąży do punktu stałego, a gdy się zbliżymy do trzech np. \(a=2.9\) układ wykonuje wiele oscylacji zanim znajdzie się w otoczeniu atraktora. Wygląda to tak jakby atraktor coraz słabiej przyciągał. Jeżeli zwiększymy parametr niewiele powyżej trójki np.: \(a=3.5\) to otrzymujemy rozwiązanie, które jest zamkniętą krzywą owijającą się jeden raz wokół niestabilnego punktu stałego, co odpowiada rozwiązaniu o okresie 2. Dla \(a=3.5\) krzywa owija się już dwa razy. Odpowiada to rozwiązaniu o okresie 4, co sugeruje, że układ pomiędzy wartościami parametru \(a=3.2\) a \(a=3.5\) przeszedł kolejną bifurkację! Ustalmy teraz parametr na największą wartość \(a=4.0\). Zachowanie się układu jest w pełni chaotyczne i nie wskazuje na obecność cykli. Możemy zwiększyć ilość iteracji lub zmienić punkt początkowy. Za każdym razem otrzymamy niepowtarzającą się trajektorię.

Czułość na warunki początkowe

Rozwiązanie jest czułe na warunki początkowe. Możemy się o tym przekonać, obliczając ciąg \(x_i\) dla dwóch mało róźniących się warunków początkowych:

Można też narysować to samo zjawisko na wykresie pajęczynowym, w tym celu możemy skonstruować nieznacznie zmienioną wersję poprzedniego programu:

Widzimy, że niezależnie od tego jak blisko siebie wystartujemy, zawsze po pewnym - i to do tego niezbyt długim czasie, rozwiązania rozbiegają się. Jest to cecha układów chaotycznych. Cytując Edwarda Lorenza:

Chaos

„… when the present determines the future, but the approximate present does not approximately determine the future …”

Diagram bifurkacyjny

Z poprzednich rysunków, widzieliśmy, że równanie logistyczne ma rosnącą do nieskończoności liczbę orbit coraz wyższego rzędu. Przy próbie ich wyznaczenia, pojawił się problem z zerami wielomianów wysokiego rzędu i musielismy się poddać. Możemy jednak bardzo prosto wyznaczyć diagram na którym będą przedstawione bifurkacje. Będziemy potrzebowali do tego dość dużej mocy obliczeniowej. Postąpimy następująco:

  1. podzielimy przedział zmienności paramteru na \(N_a\) punktów
  2. dla każdej wartości \(N_a\) przeprowadzimy jednoczesną symulację \(N_x\) układów z różnymi warunkami początkowymi, wykonami np. 1000 iteracji.
  3. zapamietamy tylko ostatnie wartości stanów tych 1000.
  4. naniesiemy punkty \((a_i,x_{1000})\) dla każdej z wybranych wartości parametru \(a\).

Jak można się domyslić wymaga to wykonania dużej ilości obliczeń. Spróbujmy „ostrożnie” z nastepującymi wartościami:

Taki wykres nazywa się diagramem bifurkacyjnym. Nazwa pochodzi od faktu, że przedstawia on bifurkacje układu przy zmianie parametru a. Porównajmy otrzymany diagram, z wykresem punktów stałych \(x=f(x)\) oraz \(x=f(f(x))\), który otrzymaliśmy poprzednio. Nalepiej przeprowadzić porównanie, umieszczając oba wykresy jeden na drugim:

Widzimy, że taki prosty algorytm umożliwił poznanie struktury punktów stałych mapy logistycznej. Jedynie stabilne punkty stałe są widoczne na takim wykresie. Wykonując wiele iteracji, zbliżamy się do tzw. atraktora układu. Czasem jest on jednym punktem, a czasem ma bardzo skomplikaowaną budowę.

Zadania:

  1. Zbadaj co się dzieje w przedziale (3,3,78)?
  2. Eksplorując diagram bifurkacyjny, czy zaobserwujesz własności samopodobieństwa atraktorów równania logistycznego?
[may76]May, R. M. „Simple mathematical models with very complicated dynamics”. Nature 261 (5560): 459–467,1976.