ś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 ! :)