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