tag:blogger.com,1999:blog-14689645565391706732024-03-05T15:11:11.085-08:00Notatki z dynamikiMarek Balcerzakhttp://www.blogger.com/profile/02116594145485113487noreply@blogger.comBlogger2125tag:blogger.com,1999:blog-1468964556539170673.post-43938636368935470092015-09-16T09:09:00.000-07:002015-09-16T09:23:42.622-07:00Analiza bifurkacji mapy logistycznej<h2>
Analiza bifurkacji mapy logistycznej</h2>
<h3>
z użyciem bibliotek SymPy, NumPy, Matplotlib</h3>
<div>
<br /></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<a href="http://www.complexityexplorer.org/">http://www.complexityexplorer.org/</a></div>
<div>
<br /></div>
<div>
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 :).</div>
<div>
<br /></div>
<div>
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:</div>
<div>
<br /></div>
<div>
x(n+1) = R*x(n)*(1-x(n))</div>
<div>
<br /></div>
<div>
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.</div>
<div>
<br /></div>
<div>
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. </div>
<div>
<br /></div>
<div>
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).</div>
<div>
<br /></div>
<div>
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). </div>
<div>
<br /></div>
<div>
Tyle teorii. W praktyce pojawia się kilka pytań:</div>
<div>
<br /></div>
<div>
P: Jak wybrać x(0) ?</div>
<div>
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). </div>
<div>
<br /></div>
<div>
P: Co to znaczy "obliczamy odpowiednio dużo wartości x(n)" ?</div>
<div>
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 :).</div>
<div>
<br /></div>
<div>
P: Co to znaczy "odrzucamy x(n) dla małych n" ?</div>
<div>
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.</div>
<div>
<br /></div>
<div>
Czas przejść do napisania kodu. Gotowe rozwiązanie, wraz z komentarzami, poniżej:<br />
<br />
<i># matplotlib.pyplot - do rysowania wykresow</i><br />
<i># numpy - do generowania rownomiernie rozlozonej listy liczb</i><br />
<i># Set - do przechowywania listy punktow na wykres</i><br />
<i><br /></i>
<i>import matplotlib.pyplot as plt</i><br />
<i>import numpy as np</i><br />
<i>from sets import Set</i><br />
<i><br /></i>
<i><br /></i>
<i># Funkcja do obliczania mapy logistycznej - zwraca liste od x(1) do x(n)</i><br />
<i>def logistic_map(n, x0, r):</i><br />
<i> x = x0</i><br />
<i> values = []</i><br />
<i> for i in range(n):</i><br />
<i> x = r*x*(1-x)</i><br />
<i> values.append(x)</i><br />
<i> return values</i><br />
<i><br /></i>
<i># Zbior punktow wykresu</i><br />
<i>bif_diagram = Set([])</i><br />
<i><br /></i>
<i># Ile pierwszych krokow odrzucamy</i><br />
<i>minsteps = 5000</i><br />
<i># Ile krokow wykonujemy</i><br />
<i>steps = 10000</i><br />
<i># Warunek poczatkowy</i><br />
<i>x0 = 0.5</i><br />
<i># Przestrzen parametru R - 2001 liczb rozlozonych rownomiernie od 2 do 4</i><br />
<i>r_space = np.linspace(2, 4, 2001)</i><br />
<i><br /></i>
<i># Dla kazdego "R":</i><br />
<i>for r in r_space:</i><br />
<i> # Drukujemy "R"</i><br />
<i> print(r)</i><br />
<i> # Obliczamy mape logistyczna</i><br />
<i> points_array = logistic_map(steps, x0, r)</i><br />
<i> # Odrzucamy pierwsze 5000 (minsteps) krokow</i><br />
<i> points_array = points_array[minsteps:]</i><br />
<i> # Dodajemy punkty do wykresu</i><br />
<i> for y in points_array:</i><br />
<i> bif_diagram.add((r, y))</i><br />
<i><br /></i>
<i>r_array = []</i><br />
<i>y_array = []</i><br />
<i><br /></i>
<i># Rozdzielamy zbior punktow na liste wspolrzednych "r" i "y"</i><br />
<i>for point in bif_diagram:</i><br />
<i> (r, y) = point</i><br />
<i> r_array.append(r)</i><br />
<i> y_array.append(y)</i><br />
<i><br /></i>
<i># Rysujemy rezultat</i><br />
<br />
<i>plt.plot(r_array, y_array, '.', markersize=1)</i></div>
<div>
<br />
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 <a href="http://continuum.io/downloads" target="_blank">tutaj</a> - 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).<br />
<br />
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):<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjfRi4LHGhuDPb4XEFU1iKq1JTt2Gh-eRtkCQCKOckErrKZ6erhS8iO82S-_rFVoelDTijAW-BNtV3PssC9Tbx7yV6v9sIKBpjRPuW5Fwxo_8EGE5yEeG8ZwCpc0CJTS9gUXiaRg9I46fYQ/s1600/screen_spyder.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="225" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjfRi4LHGhuDPb4XEFU1iKq1JTt2Gh-eRtkCQCKOckErrKZ6erhS8iO82S-_rFVoelDTijAW-BNtV3PssC9Tbx7yV6v9sIKBpjRPuW5Fwxo_8EGE5yEeG8ZwCpc0CJTS9gUXiaRg9I46fYQ/s400/screen_spyder.png" width="400" /></a></div>
<br />
<br /></div>
<div>
Efekt działania kodu wygląda następująco:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiKuXj5qvCp7978XFbkwKT3gRSDbZ3ohaViYXfUDw4WY7w8LwdgizwGAKYa6mVqq3HArJ8Iama8lhxFFJPznqEpKDdUTN9SlBURHrPGTJd887qAVYmM8gMlUsemNCoT8lTrCdgN_2JoO8d2/s1600/bifurek.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="237" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiKuXj5qvCp7978XFbkwKT3gRSDbZ3ohaViYXfUDw4WY7w8LwdgizwGAKYa6mVqq3HArJ8Iama8lhxFFJPznqEpKDdUTN9SlBURHrPGTJd887qAVYmM8gMlUsemNCoT8lTrCdgN_2JoO8d2/s400/bifurek.png" width="400" /></a></div>
<br />
Wykres warto zobaczyć w powiększeniu. Można na nim zauważyć kilka ciekawych rzeczy:<br />
<br />
- 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")<br />
<br />
- 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").<br />
<br />
- 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.<br />
<br />
- 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").<br />
<br />
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...<br />
<br />
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:<br />
<br />
<i><u>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 ?</u></i><br />
<br />
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").<br />
<br />
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 :)).<br />
<br />
Parafrazując treść zadania można powiedzieć, że szukamy takiej <i>najmniejszej</i> wartości R, dla której istnieje wartość x(n) spełniająca równanie:<br />
<br />
x(n) = x(n+2)<i> (warunek cyklu granicznego o okresie 2)</i><br />
<br />
dla wszystkich odpowiednio dużych "n", a przy tym nie spełnia równania:<br />
<br />
x(n) = x(n+1) <i>(warunek punktu krytycznego)</i><br />
<br />
Czas użyć SymPy - pakietu Pythona do obliczeń symbolicznych (fajny tutorial jest <a href="http://docs.sympy.org/latest/tutorial/index.html" target="_blank">tutaj</a> ). 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".<br />
<br />
Oczywiście każdą komendę potwierdzamy klawiszem Enter :).<br />
<br />
Najpierw importujemy całą zawartość modułu SymPy:<br />
<br />
<i>from sympy import *</i><br />
<br />
Następnie definiujemy potrzebne symbole (ich znaczenie wyjaśni się dalej):<br />
<br />
<i>(a, b, c, d, R, x, x1, x2) = symbols('a b c d R x x1 x2')</i><br />
<br />
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:<br />
<br />
<i>x1 = R*x*(1-x)</i><br />
<br />
oraz<br />
<br />
<i>x2 = R*x1*(1-x1)</i><br />
<i><br /></i>
W sympy jako x1 oznaczyliśmy x(n+1), zaś pod symbolem x2 mamy x(n+2).<br />
<br />
Możemy teraz podejrzeć zawartość symbolu x2 (odpowiadającego x(n+2)) pisząc w konsoli:<br />
<br />
<i>x2</i><br />
<i><br /></i>
lub, jeśli chcemy postać ładnie rozłożoną na czynniki:<br />
<i><br /></i>
<i>factor(x2)</i><br />
<i><br /></i>
Otrzymujemy następujący rezultat:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh0_1DlVhgliArkL-emb-4hnori3PUBr7d5gJFQgjiZiB96woJm1beysJA-M_ogIgGv8EEt2xeqovYT-L0JMRPQv27OA5RjgHN1NTKsK3iJBWzunmUshvY66YTZJNVsm0AFnSh2YoPvmec3/s1600/x_n_2.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh0_1DlVhgliArkL-emb-4hnori3PUBr7d5gJFQgjiZiB96woJm1beysJA-M_ogIgGv8EEt2xeqovYT-L0JMRPQv27OA5RjgHN1NTKsK3iJBWzunmUshvY66YTZJNVsm0AFnSh2YoPvmec3/s1600/x_n_2.PNG" /></a></div>
Rozwiążmy teraz równania:<br />
<br />
x(n+2) - x(n) = 0<br />
<br />
oraz<br />
<br />
x(n+1) - x(n) = 0<br />
<br />
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.<br />
<br />
Piszemy więc w konsoli:<br />
<i><br /></i>
<i>solve(x2 - x, R)</i><br />
<i><br /></i>
oraz<br />
<i><br /></i>
<i>solve(x1 - x, R)</i><br />
<br />
Znalezione wartości parametru R są następujące:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg5ZFzqKCKZdoSrBuZj8KXQgxRkk4AFQSv9rVsipj4t5-gjBF4PddSLIaMxM5oUSEsW3tdLBpSgVAI96nlqBYfLLZ0jrGB2mszghxzRTjIWX99HJTgGip8ZF1myzWGNDpc8uku1KIQG9svz/s1600/solve.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="125" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg5ZFzqKCKZdoSrBuZj8KXQgxRkk4AFQSv9rVsipj4t5-gjBF4PddSLIaMxM5oUSEsW3tdLBpSgVAI96nlqBYfLLZ0jrGB2mszghxzRTjIWX99HJTgGip8ZF1myzWGNDpc8uku1KIQG9svz/s400/solve.PNG" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
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.</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
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.</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
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:</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
<i>solve(-3*x**2 + 2*x + 1, x)</i></div>
<br />
Równanie ma 2 rozwiązania:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj4AFvGuJlFVLf18Vh0TlnIwSlx0efXP2vhdQFSboCK_I4oUFh555PfCPlVsrdZsepXRU2w9sUCKow_zFv8OhmanHWJkL5eUpVGdOwhj9scy7RG-cLs7PhBazd618mVG8Y1Vz80-Y6O7t8f/s1600/solve2.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj4AFvGuJlFVLf18Vh0TlnIwSlx0efXP2vhdQFSboCK_I4oUFh555PfCPlVsrdZsepXRU2w9sUCKow_zFv8OhmanHWJkL5eUpVGdOwhj9scy7RG-cLs7PhBazd618mVG8Y1Vz80-Y6O7t8f/s1600/solve2.PNG" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
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.</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
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.<br />
<br />
Tyle dyskusji o dziedzinie, teraz szukamy <b><i><u>najmniejszej możliwej</u></i></b> wartości R spełniającej 3 warunki:<br />
- Należy do przedziału (0, 4), zgodnie z założeniami mapy logistycznej<br />
- Jest rozwiązaniem równania x(n+2) - x(n) = 0<br />
- Nie jest rozwiązaniem równania x(n+1) - x(n)<br />
<br />
Spójrzmy jeszcze raz na rozwiązania równania x(n+2) - x(n) = 0:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjBLjVWoJH3C_QzxUkUeUetaLzRjhUvDnFf2DIDesLQvF0mPsKJ7n_mLNCLYcYPp21t1vKqApkNzpn8IUIIw4isJK7Cn3Qo8qRzhhBpGZhd-1aT0QbwldPA8NyOGGEqyJ6gX5GfI7gfjExK/s1600/solve.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="125" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjBLjVWoJH3C_QzxUkUeUetaLzRjhUvDnFf2DIDesLQvF0mPsKJ7n_mLNCLYcYPp21t1vKqApkNzpn8IUIIw4isJK7Cn3Qo8qRzhhBpGZhd-1aT0QbwldPA8NyOGGEqyJ6gX5GfI7gfjExK/s400/solve.PNG" width="400" /></a></div>
<br />
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 "-".<br />
<br />
Użyjmy następującej komendy:<br />
<br />
<i>(a, b, c) = solve(x2 - x, R)</i><br />
<i><br /></i>
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":<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhaAYimz4njQ-ugvpjIQbTNkAMzkmHNvt_ahjvDAqk5cID-7c2P4w07t3g4tkf1E-Wvwuq1wKjYFJghyphenhyphenFksM_cz_EaoXz2kgUqUwRbB0q8fm0UoJH2r1nU9PD8XbGeLLmbv7iZVepluQCKf/s1600/solve3.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhaAYimz4njQ-ugvpjIQbTNkAMzkmHNvt_ahjvDAqk5cID-7c2P4w07t3g4tkf1E-Wvwuq1wKjYFJghyphenhyphenFksM_cz_EaoXz2kgUqUwRbB0q8fm0UoJH2r1nU9PD8XbGeLLmbv7iZVepluQCKf/s1600/solve3.PNG" /></a></div>
<br />
Chcemy aby R było jak najmniejsze, liczymy więc pochodną wyrażenia na R (z rozwiązania "a") po x:<br />
<br />
<i>d = diff(a, x)</i><br />
<i><br /></i>
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjVum0_rT18bkqiKiqXeA2tqWzOpaI8yuoW7nDmlRmGJtlppK-HekpNhkKmBUjziL0tkTkd_G8meSOeVf7h2V6HRwRsfVZN-0kV0uya5caJjRg4omuuy6pq9MRSQF7-iLlu_DEoODapquwi/s1600/diff.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="92" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjVum0_rT18bkqiKiqXeA2tqWzOpaI8yuoW7nDmlRmGJtlppK-HekpNhkKmBUjziL0tkTkd_G8meSOeVf7h2V6HRwRsfVZN-0kV0uya5caJjRg4omuuy6pq9MRSQF7-iLlu_DEoODapquwi/s400/diff.PNG" width="400" /></a></div>
Sprawdzamy dla jakiego x pochodna jest równa 0:<br />
<br />
<i>solve(d, x)</i><br />
<br />
i otrzymujemy:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgeI6cJImgDc7kBy1FvpY_xS0pubhDANmoiVg7bYqY2CGxs5dHYm2PWnMonrRPvRzRFLBbPeJ0USW-xYuAD6NK8ul7k_begBrtDrHcZuV3jxP-6I5NECbxl3aWPVxW9QKjeiUKEdiHr2sUl/s1600/solve4.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgeI6cJImgDc7kBy1FvpY_xS0pubhDANmoiVg7bYqY2CGxs5dHYm2PWnMonrRPvRzRFLBbPeJ0USW-xYuAD6NK8ul7k_begBrtDrHcZuV3jxP-6I5NECbxl3aWPVxW9QKjeiUKEdiHr2sUl/s1600/solve4.PNG" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
Aby odnaleźć najmniejszą wartość R, podstawiamy x=2/3 do wyrażenia na R (z rozwiązania "a"):</div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgEL8ZDxR145QMQQasqi9MT6cpInYY43kqFXcE_BL_ni7OkYg0EXpHinR-GVioP_7AjuB0-b8cWDGH15lPsmZb4_qNjN-rXEjW7hO_sAFIGOoABkgw7z2irtQOhpsN1UYxORffsAye0lmoz/s1600/subs.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgEL8ZDxR145QMQQasqi9MT6cpInYY43kqFXcE_BL_ni7OkYg0EXpHinR-GVioP_7AjuB0-b8cWDGH15lPsmZb4_qNjN-rXEjW7hO_sAFIGOoABkgw7z2irtQOhpsN1UYxORffsAye0lmoz/s1600/subs.PNG" /></a></div>
W ten sposób potwierdziliśmy, że najmniejszą wartość R, przy której istnieje cykl graniczny jest 3 :).<br />
<br />
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 :).<br />
<br /></div>
Marek Balcerzakhttp://www.blogger.com/profile/02116594145485113487noreply@blogger.com0tag:blogger.com,1999:blog-1468964556539170673.post-31904508253756877972015-07-08T08:14:00.000-07:002015-07-08T08:19:34.705-07:00Analiza nieliniowego układu dynamicznego w Pythonie<h2>
Analiza nieliniowego układu dynamicznego w Pythonie</h2>
<h3>
z użyciem bibliotek SymPy, NumPy, Matplotlib</h3>
<div>
<br /></div>
<div>
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 !</div>
<div>
<br /></div>
<div>
Jeśli jeszcze nie zainstalowaliście Pythona, to pod Windows najlepiej moim zdaniem skorzystać z pakietu Anaconda (<a href="http://continuum.io/downloads">do pobrania tutaj</a> - 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).</div>
<div>
<br /></div>
<div>
Dla osób które Pythona zupełnie nie znają, polecam zacząć od tutoriala (<a href="https://docs.python.org/2/tutorial/">tutaj</a>). Oczywiście poniższe ćwiczenie można wykonać także będąc zupełnie "zielonym" w temacie :).</div>
<div>
<br /></div>
<div>
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.</div>
<div>
<br /></div>
<div>
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.</div>
<div>
<br /></div>
<div>
Rozpracujemy zdanie 6.5.1 ze strony 185. Treść jest następująca:</div>
<div>
<br /></div>
<div>
Rozważmy układ x'' = x**3 - x</div>
<div>
a) Znaleźć wszystkie położenia równowagi i zakwalifikować je</div>
<div>
b) Znaleźć "wielkość zachowaną" (energię)</div>
<div>
c) Narysować portret fazowy</div>
<div>
<br /></div>
<div>
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:</div>
<div>
<br /></div>
<div>
y = x'</div>
<div>
<br /></div>
<div>
Otrzymujemy układ dwóch równań:</div>
<div>
<br /></div>
<div>
x' = y</div>
<div>
y' = x^3 - x</div>
<div>
<br /></div>
<div>
Przedstawmy prawe strony równań w postaci funkcji argumentów x, y:</div>
<div>
<br /></div>
<div>
f(x, y) = y</div>
<div>
g(x, y) = x**3 - x</div>
<div>
<br /></div>
<div>
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.</div>
<div>
<br /></div>
<div>
Punkt równowagi spełnia warunki:</div>
<div>
<br /></div>
<div>
x' = f(x, y) = 0</div>
<div>
oraz</div>
<div>
y' = g(x, y) = 0</div>
<div>
<br /></div>
<div>
Czas użyć Pythona - z menu Start, folder Anaconda, wybieram IPython (Py 2.7) QTConsole. Naszym oczom ukazuje się taki widok:</div>
<div>
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjiu9nBCoVzDXkmObetp8GylAJuISyHreirlogWV8bloUBBf8rxYR8viJ03KgtuU8dyGXb6cspMQ0Zn7zeHc2Ke9D_Gl7H0QKSnT3NxbT6330DGNNw-Q6n2kpOadmpeQ0uFq6RiUBnmJ9PL/s1600/konsola.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="255" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjiu9nBCoVzDXkmObetp8GylAJuISyHreirlogWV8bloUBBf8rxYR8viJ03KgtuU8dyGXb6cspMQ0Zn7zeHc2Ke9D_Gl7H0QKSnT3NxbT6330DGNNw-Q6n2kpOadmpeQ0uFq6RiUBnmJ9PL/s320/konsola.PNG" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div>
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ę:</div>
<div>
<br /></div>
<div>
<i>from sympy import *</i></div>
<div>
<br /></div>
<div>
(Wszystkie komendy, które musimy wpisać do konsoli, będę podawał kursywą. To co wyrzuci konsola zapiszę czcionką pogrubioną).</div>
<div>
<br /></div>
<div>
Deklarujemy x, y oraz f, g jako symbole:</div>
<div>
<br /></div>
<div>
<i>f, g, x, y = symbols('f g x y')</i></div>
<div>
<br /></div>
<div>
Przypisujemy wzory funkcji f i g:</div>
<div>
<br /></div>
<div>
<i>f = y</i></div>
<div>
<i>g = x**3 - x</i></div>
<div>
<br />
Oczywiście po każdej komendzie naciskamy Enter, powinno to wyglądać mniej więcej tak:</div>
<div>
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhf8StfKgHiS5oeJvcu_Ls7eaahj2l1Zv0RIz-J41W3sBHRW9rIKFHlO8qb9FpGIC_Mx_3nZNdd769ElJNmclsLGrCNlaONGMIpx-7zP_TPofetIxcjuWtlwpEFMPs0mEvz9SoCJ7QF5Oqp/s1600/komendy_def.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhf8StfKgHiS5oeJvcu_Ls7eaahj2l1Zv0RIz-J41W3sBHRW9rIKFHlO8qb9FpGIC_Mx_3nZNdd769ElJNmclsLGrCNlaONGMIpx-7zP_TPofetIxcjuWtlwpEFMPs0mEvz9SoCJ7QF5Oqp/s1600/komendy_def.PNG" /></a></div>
<div>
<br /></div>
<div>
<i><br /></i></div>
<div>
Następnie rozwiązujemy układ f(x,y)=0 i g(x,y)=0:</div>
<div>
<br /></div>
<div>
<i>solve([Eq(f, 0), Eq(g, 0)], [x, y])</i></div>
<div>
<br /></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<b>[(-1, 0), (0, 0), (1, 0)]</b></div>
<div>
<b><br /></b></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<i>A = Matrix([[diff(f, x), diff(f, y)], [diff(g, x), diff(g, y)]])</i></div>
<div>
<i><br /></i></div>
<div>
Powyższa komenda przypisuje do macierzy A następujące elementy:</div>
<div>
- W pierwszej kolumnie, pierwszym wierszu: pochodna f po x</div>
<div>
- W drugiej kolumnie, pierwszym wierszu: pochodna f po y</div>
<div>
<br /></div>
<div>
i tak dalej.</div>
<div>
<i><br /></i></div>
<div>
Żeby zobaczyć rezultat, wpisujemy w konsoli:</div>
<div>
<br /></div>
<div>
<i>A</i></div>
<div>
<br /></div>
<div>
Otrzymujemy następujący wynik:</div>
<div>
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjBT3gMhz6baftSYQNnbOh50k0xEgXgDNUgN9CrnhE3vhSz12m5RGR4tZK2dJA_AHGn3l35ld6vNvYr3r3s_exKdzbfFsCzCYEA18rL4iRejigOCXjgeCNf7TkayAeUDz5DdtkaWbnFtOap/s1600/jakobian.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="83" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjBT3gMhz6baftSYQNnbOh50k0xEgXgDNUgN9CrnhE3vhSz12m5RGR4tZK2dJA_AHGn3l35ld6vNvYr3r3s_exKdzbfFsCzCYEA18rL4iRejigOCXjgeCNf7TkayAeUDz5DdtkaWbnFtOap/s400/jakobian.PNG" width="400" /></a></div>
<div>
<br /></div>
<div>
W tym zadaniu jakobian zależy tylko od zmiennej x, ale ogólnie może on zależeć od x i y.</div>
<div>
<br /></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<i>A.subs({x: -1, y: 0})</i></div>
<div>
<br /></div>
<div>
Komenda ta podstawia -1 pod x oraz 0 pod y (gdyby y znajdował się w jakobianie). Wynik jest następujący:</div>
<div>
<br /></div>
<div>
<div>
<b>Matrix([</b></div>
<div>
<b>[0, 1],</b></div>
<div>
<b>[2, 0]])</b></div>
</div>
<div>
<br /></div>
<div>
Wartości własne macierzy oblicza się następująco:</div>
<div>
<br /></div>
<div>
<i>A.subs({x: -1, y: 0}).eigenvals()</i></div>
<div>
<i><br /></i></div>
<div>
Wynik:</div>
<div>
<br /></div>
<div>
<b>{-sqrt(2): 1, sqrt(2): 1}</b></div>
<div>
<b><br /></b></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<i>A.subs({x: -1, y: 0}).eigenvects()</i></div>
<div>
<i><br /></i></div>
<div>
<div>
<b>[(-sqrt(2), 1, [Matrix([</b></div>
<div>
<b> [-sqrt(2)/2],</b></div>
<div>
<b> [ 1]])]), (sqrt(2), 1, [Matrix([</b></div>
<div>
<b> [sqrt(2)/2],</b></div>
<div>
<b> [ 1]])])]</b></div>
</div>
<div>
<b><br /></b></div>
<div>
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.</div>
<div>
<br /></div>
<div>
Sprawdzamy punkt (0, 0):</div>
<div>
<br /></div>
<div>
<i>A.subs({x: 0, y: 0}).eigenvals()</i></div>
<div>
<i><br /></i></div>
<div>
Otrzymujemy:</div>
<div>
<br /></div>
<div>
<b>{I: 1, -I: 1}</b></div>
<div>
<b><br /></b></div>
<div>
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). </div>
<div>
<br /></div>
<div>
Na koniec punkt (1, 0):</div>
<div>
<br /></div>
<div>
<i>A.subs({x: 1, y: 0}).eigenvals()</i></div>
<div>
<i><br /></i></div>
<div>
Wynik:</div>
<div>
<br /></div>
<div>
<b>{-sqrt(2): 1, sqrt(2): 1}</b></div>
<div>
<b><br /></b></div>
<div>
Tak jak w (-1, 1) mamy dwie rzeczywiste wartości własne o przeciwnych znakach, czyli znów siodło.</div>
<div>
<br /></div>
<div>
Podpunkt a) rozwiązany, czas na podpunkt b). Ponieważ nasz układ</div>
<div>
<br /></div>
<div>
x'' = x**3 - x</div>
<div>
<br /></div>
<div>
jest w postaci</div>
<div>
<br /></div>
<div>
x'' = F(x)</div>
<div>
<br /></div>
<div>
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:</div>
<div>
<br /></div>
<div>
F(x) = x**3 - x = dV/dx</div>
<div>
<br /></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<i>F = symbols('F')</i></div>
<div>
<i>F = x**3-x</i></div>
<div>
<i>Integral(F, x).doit()</i></div>
<div>
<br /></div>
<div>
Otrzymujemy:</div>
<div>
<br /></div>
<div>
<b>x**4/4 - x**2/2</b></div>
<div>
<br /></div>
<div>
Do energii potencjalnej doliczamy jeszcze kinetyczną (x'**2/2), a więc całkowita energia w układzie wynosi:</div>
<div>
<br /></div>
<div>
E = x'**2/2 + x^4/4 - x^2/2</div>
<div>
<br /></div>
<div>
Do zrobienia został podpunkt c) - rysowanie portretu fazowego. Do tego posłuży nam funkcja streamplot() z pakietu Matplotlib (przykłady <a href="http://matplotlib.org/users/screenshots.html#streamplot">tutaj</a>). </div>
<div>
<br /></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<i>import numpy as np</i></div>
<div>
<i><br /></i></div>
<div>
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.</div>
<div>
<br /></div>
<div>
<i>fl = lambdify((x,y), f)</i></div>
<div>
<i>gl = lambdify((x,y), g)</i></div>
<div>
<i><br /></i></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<i>gl(0.5, 0.5)</i></div>
<div>
<br /></div>
<div>
Otrzymujemy</div>
<div>
<br /></div>
<div>
<b>-0.375</b></div>
<div>
<b><br /></b></div>
<div>
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.</div>
<div>
<br /></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<i>import matplotlib.pyplot as plt</i></div>
<div>
<i><br /></i></div>
<div>
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:</div>
<div>
<br /></div>
<div>
<i>Y, X = np.mgrid[-3:3:100j, -3:3:100j]</i></div>
<div>
<br /></div>
<div>
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).</div>
<div>
<br /></div>
<div>
Obliczamy wartości funkcji f i g dla wszystkich punktów z siatki, korzystając z przygotowanych funkcji Pythona fl i gl:</div>
<div>
<br /></div>
<div>
<i>U = fl(X, Y)</i></div>
<div>
<i>V = gl(X, Y)</i></div>
<div>
<br /></div>
<div>
Na koniec rysuję portret fazowy i wyświetlam go:</div>
<div>
<br /></div>
<div>
<i>plt.streamplot(X, Y, U, V)</i></div>
<div>
<i>plt.show()</i></div>
<div>
<i><br /></i></div>
<div>
Naszym oczom ukazuje się piękny wykres:</div>
<div>
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHCIE-C6UTBueeDRyXxQN8VC_tn61-vu8pQZJpyjtzGHQGz4d-RtyezgSMuwKt5E8fn48m-YCUrJU6kIhndxq-s6GIP2tJeLivCBHGbiJ9kA7WPigOoKSvYysAGJFXvkOxo3g5wSdtJ9MB/s1600/figure_1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="297" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHCIE-C6UTBueeDRyXxQN8VC_tn61-vu8pQZJpyjtzGHQGz4d-RtyezgSMuwKt5E8fn48m-YCUrJU6kIhndxq-s6GIP2tJeLivCBHGbiJ9kA7WPigOoKSvYysAGJFXvkOxo3g5wSdtJ9MB/s400/figure_1.png" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div>
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.</div>
<div>
<br /></div>
<div>
Mam nadzieję, że ten wpis zachęci do dalszych poszukiwań :). Zostawiam linki do tutoriali:</div>
<div>
<br /></div>
<div>
<a href="http://docs.sympy.org/latest/tutorial/index.html">SymPy</a></div>
<div>
<a href="http://wiki.scipy.org/Tentative_NumPy_Tutorial">NumPy</a></div>
<div>
<a href="http://matplotlib.org/users/beginner.html">Matplotlib</a></div>
<div>
<br /></div>
<div>
Powodzenia ! :)</div>
<div>
<br /></div>
Marek Balcerzakhttp://www.blogger.com/profile/02116594145485113487noreply@blogger.com0Łódź, Polska51.7592485 19.45598330000007151.602042499999996 19.133259800000072 51.9164545 19.778706800000069