Analiza bifurkacji mapy logistycznej
z użyciem bibliotek SymPy, NumPy, Matplotlib
Witam serdecznie po dłuższej przerwie spowodowanej wakacjami i nieuleczalnym lenistwem :). Do napisania dzisiejszego posta zainspirował mnie kurs dynamiki nieliniowej dostępny pod adresem:
który bardzo serdecznie polecam :). Tam właśnie natknąłem się na interesujące zadanie dotyczące mapy logistycznej, które można rozwiązywać dobrą godzinę na kartce, lub w 10 minut używając Pythona. Tutaj zaprezentuję drugie podejście :).
Najpierw jednak parę słów o tym co to jest mapa logistyczna (zwana też odwzorowaniem logistycznym). Jest to bardzo prosty układ dynamiczny czasu dyskretnego z jedną zmienną stanu x(n). Układ ten opisany jest równaniem:
x(n+1) = R*x(n)*(1-x(n))
gdzie n - liczba naturalna, R - liczba rzeczywista (parametr) z przedziału (0, 4). Jeśli warunek początkowy x(0) należy do przedziału [0, 1], to x(n) na zawsze (dla dowolnego n) w tym przedziale zostanie.
Zaczniemy od narysowania wykresu bifurkacyjnego dla powyższego układu dynamicznego. Wykres bifurkacyjny pokazuje co dzieje się ze zmienną stanu układu przy czasie dążącym do nieskończoności w zależności od parametru układu.
Dlaczego miałoby nas to interesować ? Okazuje się, że mapa logistyczna w zależności od parametru R może ukazywać kilka zupełnie różnych zachowań w czasie dążącym do nieskończoności. Układ może "utknąć" na pewnej wartości tak, że x(n+1) = x(n) dla odpowiednio dużego n. Mogą pojawić się też rozwiązania okresowe lub aperiodyczne (chaotyczne).
Jak stworzyć wykres bifurkacyjny ? Bardzo prosto :). Jako oś poziomą wybieramy parametr R, jako oś pionową - wartości x(n). Dla każdej interesującej nas wartości R wybieramy warunek początkowy x(0), a następnie obliczamy odpowiednio dużo wartości x(n). Ponieważ interesuje nas tylko zachowanie w czasie dążącym do nieskończoności, odrzucamy wartości x(n) dla małych n, a pozostałe nanosimy na wykres x(R).
Tyle teorii. W praktyce pojawia się kilka pytań:
P: Jak wybrać x(0) ?
O: Dowolnie, byle nie 0 ani 1 (dlaczego ? wystarczy spojrzeć na wzór - dla tych dwóch wartości początkowych x(n) = 0). Układ ten dla dowolnego R z przedziału (0, 4) posiada tylko jeden atraktor, a więc tylko jeden obiekt, który "ściąga" do siebie wartości x(n). Obszar przyciągania atraktora ("basin of attraction" z angielska) obejmuje cały przedział (0, 1).
P: Co to znaczy "obliczamy odpowiednio dużo wartości x(n)" ?
O: Obliczamy tyle wartości, aby dobrze widoczna była struktura topologiczna atraktora (w prostych żołnierskich słowach: musimy dokładnie widzieć, co się dzieje z x(n) dla dużych n - czy jest stałe, czy okresowo zmienne, czy nieokresowo zmienne - chaotyczne). Najprostsza recepta ? Oblicz tyle x(n) ile uważasz. a potem oblicz o połowę więcej i porównaj. Jeśli nic się nie zmieniło na wykresie - dobrze wybrałeś liczbę cykli :).
P: Co to znaczy "odrzucamy x(n) dla małych n" ?
O: Musimy odrzucić te x(n), które zawierają "odpowiedź przejściową" (ang. "transient response"), czyli oscylacje występujące zanim x(n) osiągnie atraktor. Znowu posiłkujemy się doświadczeniem.
Czas przejść do napisania kodu. Gotowe rozwiązanie, wraz z komentarzami, poniżej:
# matplotlib.pyplot - do rysowania wykresow
# numpy - do generowania rownomiernie rozlozonej listy liczb
# Set - do przechowywania listy punktow na wykres
import matplotlib.pyplot as plt
import numpy as np
from sets import Set
# Funkcja do obliczania mapy logistycznej - zwraca liste od x(1) do x(n)
def logistic_map(n, x0, r):
x = x0
values = []
for i in range(n):
x = r*x*(1-x)
values.append(x)
return values
# Zbior punktow wykresu
bif_diagram = Set([])
# Ile pierwszych krokow odrzucamy
minsteps = 5000
# Ile krokow wykonujemy
steps = 10000
# Warunek poczatkowy
x0 = 0.5
# Przestrzen parametru R - 2001 liczb rozlozonych rownomiernie od 2 do 4
r_space = np.linspace(2, 4, 2001)
# Dla kazdego "R":
for r in r_space:
# Drukujemy "R"
print(r)
# Obliczamy mape logistyczna
points_array = logistic_map(steps, x0, r)
# Odrzucamy pierwsze 5000 (minsteps) krokow
points_array = points_array[minsteps:]
# Dodajemy punkty do wykresu
for y in points_array:
bif_diagram.add((r, y))
r_array = []
y_array = []
# Rozdzielamy zbior punktow na liste wspolrzednych "r" i "y"
for point in bif_diagram:
(r, y) = point
r_array.append(r)
y_array.append(y)
# Rysujemy rezultat
plt.plot(r_array, y_array, '.', markersize=1)
# matplotlib.pyplot - do rysowania wykresow
# numpy - do generowania rownomiernie rozlozonej listy liczb
# Set - do przechowywania listy punktow na wykres
import matplotlib.pyplot as plt
import numpy as np
from sets import Set
# Funkcja do obliczania mapy logistycznej - zwraca liste od x(1) do x(n)
def logistic_map(n, x0, r):
x = x0
values = []
for i in range(n):
x = r*x*(1-x)
values.append(x)
return values
# Zbior punktow wykresu
bif_diagram = Set([])
# Ile pierwszych krokow odrzucamy
minsteps = 5000
# Ile krokow wykonujemy
steps = 10000
# Warunek poczatkowy
x0 = 0.5
# Przestrzen parametru R - 2001 liczb rozlozonych rownomiernie od 2 do 4
r_space = np.linspace(2, 4, 2001)
# Dla kazdego "R":
for r in r_space:
# Drukujemy "R"
print(r)
# Obliczamy mape logistyczna
points_array = logistic_map(steps, x0, r)
# Odrzucamy pierwsze 5000 (minsteps) krokow
points_array = points_array[minsteps:]
# Dodajemy punkty do wykresu
for y in points_array:
bif_diagram.add((r, y))
r_array = []
y_array = []
# Rozdzielamy zbior punktow na liste wspolrzednych "r" i "y"
for point in bif_diagram:
(r, y) = point
r_array.append(r)
y_array.append(y)
# Rysujemy rezultat
plt.plot(r_array, y_array, '.', markersize=1)
Jak uruchomić powyższy program ? Bardzo prosto :). Jeśli jeszcze nie zainstalowaliście Pythona, to pod Windows najlepiej moim zdaniem skorzystać z pakietu Anaconda (do pobrania tutaj - korzystam z Pythona w wersji 2.7). Użytkownicy Linuksa z pewnością poradzą sobie z instalacją sami :) (potrzebny będzie Python 2.7 oraz biblioteki: SymPy, NumPy oraz Matplotlib).
Jeśli korzystamy z Anacondy pod Windows, w menu Start wybieramy foler Anaconda, a następnie uruchamiamy środowisko Spyder. Wklejamy cały powyższy kod do pola tekstowego po lewej i klikamy przycisk Run file (lub naciskamy F5 na klawiaturze):
Efekt działania kodu wygląda następująco:
Wykres warto zobaczyć w powiększeniu. Można na nim zauważyć kilka ciekawych rzeczy:
- Dla wartości R mniejszych od około 3 mamy tylko jedną wartość x na wykresie. Oznacza to, że dla R < 3 x(n) dąży pewnej wartości i zostaje w niej na zawsze - jest to stabilny punkt krytyczny ("stable fixed point")
- Dla każdej wartości R większej od około 3, a mniejszej od około 3,45 mamy dwie wartości x na wykresie. Oznacza to, że dla tego przedziału R wartość x(n) dla n dążącego do nieskończoności "skacze" między dwiema wartościami - jest to orbita okresowa albo cykl graniczny ("limit cycle").
- Dla R większych od około 3,45 a mniejszych od 3,54 mamy także cykl graniczny, ale składający się z 4 wartości. Dla R nieco większych od 3.54 cykl graniczny zawiera już 8 wartości.
- Dla większości wartości R większych od około 3,57 zaobserwować można zachowanie chaotyczne, układ nie wykonuje powtarzalnych cykli. Atraktor chaotyczny nazywany jest także "dziwnym" ("chaotic attractor / strange attractor").
Uważny czytelnik zauważył z pewnością, że przedziały wartości R szacowałem na podstawie wykresu. Są jednak inne możliwości...
Jeśli wykres bifurkacyjny jest gotowy, warto przejść do właściwego zadania, o którym wspominałem na początku posta. Brzmi ono następująco:
Dla jakiej dokładnie wartości R zachowanie x(n) w czasie dążącym do nieskończoności zmienia się ze stabilnego punktu krytycznego na cykl graniczny ?
Wartość parametru, dla której zmienia się jakościowo zachowanie x(n) dla n dążącego do nieskończoności, nazywamy "punktem bifurkacji", a samo zjawisko - "bifurkacją". (ang. "bifurcation").
Zgodnie z obserwacją wykresu, punkt ten znajduje się w pobliżu wartości R = 3. Okazuje się, że jest to wartość dokładna, ale dla pewności przeprowadzimy dowód (z pomocą pakietu SymPy :)).
Parafrazując treść zadania można powiedzieć, że szukamy takiej najmniejszej wartości R, dla której istnieje wartość x(n) spełniająca równanie:
x(n) = x(n+2) (warunek cyklu granicznego o okresie 2)
dla wszystkich odpowiednio dużych "n", a przy tym nie spełnia równania:
x(n) = x(n+1) (warunek punktu krytycznego)
Czas użyć SymPy - pakietu Pythona do obliczeń symbolicznych (fajny tutorial jest tutaj ). Jeśli mamy nadal otwarty program Spyder, użyjmy go pisząc poniższe komendy w konsoli znajdującej się w prawym dolnym rogu okna programu. Można też użyć "IPython QT console".
Oczywiście każdą komendę potwierdzamy klawiszem Enter :).
Najpierw importujemy całą zawartość modułu SymPy:
from sympy import *
Następnie definiujemy potrzebne symbole (ich znaczenie wyjaśni się dalej):
(a, b, c, d, R, x, x1, x2) = symbols('a b c d R x x1 x2')
Musimy zdefiniować równanie x(n) = x(n+2) i podstawić do niego równanie mapy logistycznej: x(n+1) = R*x(n)*(1-x(n)). W tym celu użyjmy następujących komend:
x1 = R*x*(1-x)
oraz
x2 = R*x1*(1-x1)
W sympy jako x1 oznaczyliśmy x(n+1), zaś pod symbolem x2 mamy x(n+2).
Możemy teraz podejrzeć zawartość symbolu x2 (odpowiadającego x(n+2)) pisząc w konsoli:
x2
lub, jeśli chcemy postać ładnie rozłożoną na czynniki:
factor(x2)
Otrzymujemy następujący rezultat:
Rozwiążmy teraz równania:
x(n+2) - x(n) = 0
oraz
x(n+1) - x(n) = 0
ze względu na R. Jak już ustaliliśmy, pierwsze z nich musi być spełnione (jest to warunek cyklu granicznego o okresie 2), a drugie nie może być spełnione, bo oznaczałoby to punkt krytyczny, a nie cykl graniczny.
Piszemy więc w konsoli:
solve(x2 - x, R)
oraz
solve(x1 - x, R)
Znalezione wartości parametru R są następujące:
Równanie ma 2 rozwiązania:
Interesują nas tylko wartości x z przedziału od 0 do 1, więc wartość pod pierwiastkiem będzie dla nich na pewno dodatnia.
Tyle dyskusji o dziedzinie, teraz szukamy najmniejszej możliwej wartości R spełniającej 3 warunki:
- Należy do przedziału (0, 4), zgodnie z założeniami mapy logistycznej
- Jest rozwiązaniem równania x(n+2) - x(n) = 0
- Nie jest rozwiązaniem równania x(n+1) - x(n)
Spójrzmy jeszcze raz na rozwiązania równania x(n+2) - x(n) = 0:
Wiemy, że x należy do przedziału (0, 1), więc mianownik pierwszych dwóch rozwiązań jest na pewno ujemny. Chcemy aby R było dodatnie i należało do przedziału (0, 4), więc dla ujemnego mianownika chcemy też ujemnego licznika. Bierzemy więc pierwsze rozwiązanie, które przed pierwiastkiem ma znak "-".
Użyjmy następującej komendy:
(a, b, c) = solve(x2 - x, R)
W ten sposób pod symbolami a, b, c mamy rozwiązania równania x(n+2)-x(n) = 0. Interesuje nas rozwiązanie "a":
Chcemy aby R było jak najmniejsze, liczymy więc pochodną wyrażenia na R (z rozwiązania "a") po x:
d = diff(a, x)
Sprawdzamy dla jakiego x pochodna jest równa 0:
solve(d, x)
i otrzymujemy:
W ten sposób potwierdziliśmy, że najmniejszą wartość R, przy której istnieje cykl graniczny jest 3 :).
Jak widać używane tutaj pakiety Pythona mają całkiem duże możliwości i mogą znacznie ułatwić rozwiązywanie różnego rodzaju problemów. Zachęcam do dalszych, samodzielnych poszukiwań w tym temacie i serdecznie pozdrawiam - do napisania :).
Wykres warto zobaczyć w powiększeniu. Można na nim zauważyć kilka ciekawych rzeczy:
- Dla wartości R mniejszych od około 3 mamy tylko jedną wartość x na wykresie. Oznacza to, że dla R < 3 x(n) dąży pewnej wartości i zostaje w niej na zawsze - jest to stabilny punkt krytyczny ("stable fixed point")
- Dla każdej wartości R większej od około 3, a mniejszej od około 3,45 mamy dwie wartości x na wykresie. Oznacza to, że dla tego przedziału R wartość x(n) dla n dążącego do nieskończoności "skacze" między dwiema wartościami - jest to orbita okresowa albo cykl graniczny ("limit cycle").
- Dla R większych od około 3,45 a mniejszych od 3,54 mamy także cykl graniczny, ale składający się z 4 wartości. Dla R nieco większych od 3.54 cykl graniczny zawiera już 8 wartości.
- Dla większości wartości R większych od około 3,57 zaobserwować można zachowanie chaotyczne, układ nie wykonuje powtarzalnych cykli. Atraktor chaotyczny nazywany jest także "dziwnym" ("chaotic attractor / strange attractor").
Uważny czytelnik zauważył z pewnością, że przedziały wartości R szacowałem na podstawie wykresu. Są jednak inne możliwości...
Jeśli wykres bifurkacyjny jest gotowy, warto przejść do właściwego zadania, o którym wspominałem na początku posta. Brzmi ono następująco:
Dla jakiej dokładnie wartości R zachowanie x(n) w czasie dążącym do nieskończoności zmienia się ze stabilnego punktu krytycznego na cykl graniczny ?
Wartość parametru, dla której zmienia się jakościowo zachowanie x(n) dla n dążącego do nieskończoności, nazywamy "punktem bifurkacji", a samo zjawisko - "bifurkacją". (ang. "bifurcation").
Zgodnie z obserwacją wykresu, punkt ten znajduje się w pobliżu wartości R = 3. Okazuje się, że jest to wartość dokładna, ale dla pewności przeprowadzimy dowód (z pomocą pakietu SymPy :)).
Parafrazując treść zadania można powiedzieć, że szukamy takiej najmniejszej wartości R, dla której istnieje wartość x(n) spełniająca równanie:
x(n) = x(n+2) (warunek cyklu granicznego o okresie 2)
dla wszystkich odpowiednio dużych "n", a przy tym nie spełnia równania:
x(n) = x(n+1) (warunek punktu krytycznego)
Czas użyć SymPy - pakietu Pythona do obliczeń symbolicznych (fajny tutorial jest tutaj ). Jeśli mamy nadal otwarty program Spyder, użyjmy go pisząc poniższe komendy w konsoli znajdującej się w prawym dolnym rogu okna programu. Można też użyć "IPython QT console".
Oczywiście każdą komendę potwierdzamy klawiszem Enter :).
Najpierw importujemy całą zawartość modułu SymPy:
from sympy import *
Następnie definiujemy potrzebne symbole (ich znaczenie wyjaśni się dalej):
(a, b, c, d, R, x, x1, x2) = symbols('a b c d R x x1 x2')
Musimy zdefiniować równanie x(n) = x(n+2) i podstawić do niego równanie mapy logistycznej: x(n+1) = R*x(n)*(1-x(n)). W tym celu użyjmy następujących komend:
x1 = R*x*(1-x)
oraz
x2 = R*x1*(1-x1)
W sympy jako x1 oznaczyliśmy x(n+1), zaś pod symbolem x2 mamy x(n+2).
Możemy teraz podejrzeć zawartość symbolu x2 (odpowiadającego x(n+2)) pisząc w konsoli:
x2
lub, jeśli chcemy postać ładnie rozłożoną na czynniki:
factor(x2)
Otrzymujemy następujący rezultat:
x(n+2) - x(n) = 0
oraz
x(n+1) - x(n) = 0
ze względu na R. Jak już ustaliliśmy, pierwsze z nich musi być spełnione (jest to warunek cyklu granicznego o okresie 2), a drugie nie może być spełnione, bo oznaczałoby to punkt krytyczny, a nie cykl graniczny.
Piszemy więc w konsoli:
solve(x2 - x, R)
oraz
solve(x1 - x, R)
Znalezione wartości parametru R są następujące:
Spójrzmy na rozwiązania pierwszego równania: x(n+2) - x(n) = 0. Trzecie z nich jest identyczne jak rozwiązanie równania x(n+1) - x(n) = 0. Widać więc, że trzecie rozwiązanie oznacza wartość parametru R zapewniającą stabilny punkt równowagi, a nie cykl graniczny. Pozostają więc dwa pierwsze rozwiązania.
Określmy teraz dziedzinę pierwszych dwóch rozwiązań. W mianowniku mamy 2x^2 - 2x = 2x(x-1). Stąd wiadomo, że x nie może być równy 0 ani 1.
Pod pierwiastkiem mamy wyrażenie -3x^2 + 2x + 1, które musi być dodatnie. Będąc skrajnie leniwym, do rozwiązania równania -3x^2 + 2x + 1 = 0 także użyję SymPy:
solve(-3*x**2 + 2*x + 1, x)
Równanie ma 2 rozwiązania:
Teraz, zgodnie z licealną matematyką (serdeczne pozdrowienia dla prof. Mańko, gdyby to czytała :)), zauważamy że parabola -3x^2 + 2x + 1 ma ujemny współczynnik przy x^2, a więc ramiona paraboli są skierowane w dół, stąd dla argumentów x pomiędzy miejscami zerowymi (z przedziału (-1/3, 1)) wyrażenie -3x^2 + 2x + 1 przyjmuje wartości dodatnie i może być pierwiastkowane.
Tyle dyskusji o dziedzinie, teraz szukamy najmniejszej możliwej wartości R spełniającej 3 warunki:
- Należy do przedziału (0, 4), zgodnie z założeniami mapy logistycznej
- Jest rozwiązaniem równania x(n+2) - x(n) = 0
- Nie jest rozwiązaniem równania x(n+1) - x(n)
Spójrzmy jeszcze raz na rozwiązania równania x(n+2) - x(n) = 0:
Wiemy, że x należy do przedziału (0, 1), więc mianownik pierwszych dwóch rozwiązań jest na pewno ujemny. Chcemy aby R było dodatnie i należało do przedziału (0, 4), więc dla ujemnego mianownika chcemy też ujemnego licznika. Bierzemy więc pierwsze rozwiązanie, które przed pierwiastkiem ma znak "-".
Użyjmy następującej komendy:
(a, b, c) = solve(x2 - x, R)
W ten sposób pod symbolami a, b, c mamy rozwiązania równania x(n+2)-x(n) = 0. Interesuje nas rozwiązanie "a":
Chcemy aby R było jak najmniejsze, liczymy więc pochodną wyrażenia na R (z rozwiązania "a") po x:
d = diff(a, x)
solve(d, x)
i otrzymujemy:
Aby odnaleźć najmniejszą wartość R, podstawiamy x=2/3 do wyrażenia na R (z rozwiązania "a"):
Jak widać używane tutaj pakiety Pythona mają całkiem duże możliwości i mogą znacznie ułatwić rozwiązywanie różnego rodzaju problemów. Zachęcam do dalszych, samodzielnych poszukiwań w tym temacie i serdecznie pozdrawiam - do napisania :).
Brak komentarzy:
Prześlij komentarz