środa, 16 września 2015

Analiza bifurkacji mapy logistycznej

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)

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:



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.

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:


Aby odnaleźć najmniejszą wartość R, podstawiamy x=2/3 do wyrażenia na R (z rozwiązania "a"):

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 :).

środa, 8 lipca 2015

Analiza nieliniowego układu dynamicznego w Pythonie

Analiza nieliniowego układu dynamicznego w Pythonie

z użyciem bibliotek SymPy, NumPy, Matplotlib


W dzisiejszym (pierwszym !) poście chciałbym zaprezentować sposób rozwiązania prostego zagadnienia dynamiki nieliniowej używając Pythona i jego bibliotek. Jak można zauważyć poniżej, Python wraz z dodatkowymi bibliotekami jest często w stanie z powodzeniem zastąpić komercyjne, drogie oprogramowanie jak Matlab, Mathematica czy Maple. Zaczynamy !

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).

Dla osób które Pythona zupełnie nie znają, polecam zacząć od tutoriala (tutaj). Oczywiście poniższe ćwiczenie można wykonać także będąc zupełnie "zielonym" w temacie :).

Na warsztat wezmę niezbyt skomplikowane zadanie z mojej ulubionej książki o dynamice nieliniowej - "Nonlinear Dynamics and Chaos" autorstwa S. H. Strogatz'a - polecam serdecznie, wszystko ładnie wyjaśnione, w Internecie można znaleźć PDF.

Ponieważ w Pythonie potęgę zapisuje się za pomocą operatora **, taką konwencję będę też stosował tutaj. x**2 będzie więc oznaczał x do potęgi 2.

Rozpracujemy zdanie 6.5.1 ze strony 185. Treść jest następująca:

Rozważmy układ x'' = x**3 - x
a) Znaleźć wszystkie położenia równowagi i zakwalifikować je
b) Znaleźć "wielkość zachowaną" (energię)
c) Narysować portret fazowy

Zacznijmy od początku - podpunkt a). W celu znalezienia położeń równowagi przekształcamy równanie różniczkowe II rzędu (x'' = x**3 - x) na 2 równania różniczkowe I rzędu. Stosujemy podstawienie:

y = x'

Otrzymujemy układ dwóch równań:

x' = y
y' = x^3 - x

Przedstawmy prawe strony równań w postaci funkcji argumentów x, y:

f(x, y) = y
g(x, y) = x**3 - x

Oczywiście w tym wypadku obie funkcje zależą od jednej zmiennej, ale ogólnie w układach II rzędu nie jest to regułą, stąd taki ogólniejszy zapis.

Punkt równowagi spełnia warunki:

x' = f(x, y) = 0
oraz
y' = g(x, y) = 0

Czas użyć Pythona - z menu Start, folder Anaconda, wybieram IPython (Py 2.7) QTConsole. Naszym oczom ukazuje się taki widok:


Na początek rozwiążemy układ równań f(x,y)=0 i g(x,y)=0. Zaczniemy od zaimportowania biblioteki do obliczeń symbolicznych SymPy. W tym celu wpisujemy komendę:

from sympy import *

(Wszystkie komendy, które musimy wpisać do konsoli, będę podawał kursywą. To co wyrzuci konsola zapiszę czcionką pogrubioną).

Deklarujemy x, y oraz f, g jako symbole:

f, g, x, y = symbols('f g x y')

Przypisujemy wzory funkcji f i g:

f = y
g = x**3 - x

Oczywiście po każdej komendzie naciskamy Enter, powinno to wyglądać mniej więcej tak:



Następnie rozwiązujemy układ f(x,y)=0 i g(x,y)=0:

solve([Eq(f, 0), Eq(g, 0)], [x, y])

Oznacza to, że mamy listę dwóch równań: f = 0 i g = 0, a układ ten rozwiązujemy ze względu na x i y. Otrzymujemy wynik:

[(-1, 0), (0, 0), (1, 0)]

Nasz układ posiada więc 3 położenia równowagi. Aby je zakwalifikować, dla każdego z nich liczymy jakobian, czyli macierz pochodnych. Jakobian naszego układu obliczamy w następujący sposób:

A = Matrix([[diff(f, x), diff(f, y)], [diff(g, x), diff(g, y)]])

Powyższa komenda przypisuje do macierzy A następujące elementy:
- W pierwszej kolumnie, pierwszym wierszu: pochodna f po x
- W drugiej kolumnie, pierwszym wierszu: pochodna f po y

i tak dalej.

Żeby zobaczyć rezultat, wpisujemy w konsoli:

A

Otrzymujemy następujący wynik:


W tym zadaniu jakobian zależy tylko od zmiennej x, ale ogólnie może on zależeć od x i y.

Podstawiamy teraz do jakobianu kolejno punkty w których f(x,y)=0 i g(x,y)=0 i w każdym z tych punktów obliczymy wartości własne i wektory własne. Dla punktu (-1, 0) mamy:

A.subs({x: -1, y: 0})

Komenda ta podstawia -1 pod x oraz 0 pod y (gdyby y znajdował się w jakobianie). Wynik jest następujący:

Matrix([
[0, 1],
[2, 0]])

Wartości własne macierzy oblicza się następująco:

A.subs({x: -1, y: 0}).eigenvals()

Wynik:

{-sqrt(2): 1, sqrt(2): 1}

Mamy dwie rzeczywiste wartości własne o przeciwnych znakach, a więc punkt (-1, 0) jest tzw. "siodłem". Jedynki w wyniku oznaczają, że obie wartości własne są pojedyncze. Możemy też obliczyć wektory własne:

A.subs({x: -1, y: 0}).eigenvects()

[(-sqrt(2), 1, [Matrix([
   [-sqrt(2)/2],
   [         1]])]), (sqrt(2), 1, [Matrix([
   [sqrt(2)/2],
   [        1]])])]

Wynik mówi nam, że pierwsza wartość własna to -2^(1/2), występuje ona raz, zaś wektor własny odpowiadający jej to [-0.5*2^(1/2), 1]. Podobnie z drugą wartością własną i drugim wektorem własnym.

Sprawdzamy punkt (0, 0):

A.subs({x: 0, y: 0}).eigenvals()

Otrzymujemy:

{I: 1, -I: 1}

Mamy dwie urojone wartości własne o przeciwnych znakach (I - jednostka urojona), stąd dookoła tego punktu powinny wystąpić zamknięte orbity (zgodnie z tw. 6.6.1 z książki Strogatz'a). 

Na koniec punkt (1, 0):

A.subs({x: 1, y: 0}).eigenvals()

Wynik:

{-sqrt(2): 1, sqrt(2): 1}

Tak jak w (-1, 1) mamy dwie rzeczywiste wartości własne o przeciwnych znakach, czyli znów siodło.

Podpunkt a) rozwiązany, czas na podpunkt b). Ponieważ nasz układ

x'' = x**3 - x

jest w postaci

x'' = F(x)

to w układzie nie ma tłumienia, a więc zgodnie z podrozdziałem 6,5 z książki Strogatza energia jest zachowana. Siła F(x) w naszym układzie jest więc pochodną potencjału po zmiennej x:

F(x) = x**3 - x = dV/dx

Aby obliczyć energię potencjalną całkujemy więc F(x) = x**3 - x po zmiennej x. Oczywiście można to łatwo zrobić "ręcznie", ale lepiej posłużyć się Pythonem. Wpisujemy więc w konsoli:

F = symbols('F')
F = x**3-x
Integral(F, x).doit()

Otrzymujemy:

x**4/4 - x**2/2

Do energii potencjalnej doliczamy jeszcze kinetyczną (x'**2/2), a więc całkowita energia w układzie wynosi:

E = x'**2/2 + x^4/4 - x^2/2

Do zrobienia został podpunkt c) - rysowanie portretu fazowego. Do tego posłuży nam funkcja streamplot() z pakietu Matplotlib (przykłady tutaj). 

Najpierw przekształcamy funkcje f i g na funkcje Pythona używające pakietu do obliczeń numerycznych NumPy - ułatwi to nam ich wywoływanie. Importujemy pakiet numpy:

import numpy as np

Zamieniamy wyrażenia f i g na funkcje Pythona korzystające z pakietu NumPy - dzięki temu zabiegowi będziemy mogli wywoływać te funkcje tak, jakby były zdefiniowane w kodzie Pythona. W tym celu używamy funkcji lambdify.

fl = lambdify((x,y), f)
gl = lambdify((x,y), g)

Teraz fl i gl stały się funkcjami Pythona o argumentach x, y, które zachowują się tak, jak nasze funkcje f i g. Można to sprawdzić wpisując w konsoli:

gl(0.5, 0.5)

Otrzymujemy

-0.375

a więc tyle, ile oczekiwaliśmy ( g(0.5, 0.5) = 0.5^3 - 0.5 = -0.375 ). Podobnie funkcje fl(1, 0) i gl(1, 0) zwracają 0, bo (1, 0) to położenie równowagi.

Mamy gotowe funkcje obliczające wartości f(x,y) i g(x,y), można więc przystąpić do rysowania portretu fazowego. Importujemy bibliotekę pyplot z pakietu Matplotlib:

import matplotlib.pyplot as plt

Czas przygotować siatkę punktów pod wykres. Zakładam, że wystarczy nam siatka 100x100 punktów. Musi ona obejmować wszystkie punkty równowagi wraz z otoczeniem, aby można było zobaczyć jak zachowuje się nasz wspaniały układ :). Tworzę więc siatkę 100x100 punktów rozciągniętą na przedziałach: -3 <= x <= 3, -3 <= y <= 3:

Y, X = np.mgrid[-3:3:100j, -3:3:100j]

Litera "j" oznacza, że na oba końce przedziałów wchodzą do siatki (tak więc największa wartość z X będzie równa 3, najmniejsza -3).

Obliczamy wartości funkcji f i g dla wszystkich punktów z siatki, korzystając z przygotowanych funkcji Pythona fl i gl:

U = fl(X, Y)
V = gl(X, Y)

Na koniec rysuję portret fazowy i wyświetlam go:

plt.streamplot(X, Y, U, V)
plt.show()

Naszym oczom ukazuje się piękny wykres:


Wykres potwierdza to, co przewidzieliśmy rozwiązując podpunkt a) - w punktach (-1,0) oraz (1,0) mamy siodła, zaś dookoła punktu (0, 0) mamy orbity zamknięte.

Mam nadzieję, że ten wpis zachęci do dalszych poszukiwań :). Zostawiam linki do tutoriali:


Powodzenia ! :)