Przypomnienie wiadomości o symulowaniu ruchów Browna metodą Monte Carlo

W e‑materiale Modelowanie ruchów BrownaP14JSHgfzModelowanie ruchów Browna przeprowadziliśmy proste modelowanie, w którym zakładaliśmy, że w każdej jednostce czasu cząsteczka przemieszcza się o określoną z góry odległość. Zaimplementujemy ten algorytm w języku Python.

Aby zobrazować proces losowego poruszania się cząsteczki w płynie, zakładamy, że każda kolizja powoduje jej przemieszczenie się o określoną długość w dowolnym kierunku. Wybieramy losowy kąt z przedziału  0 , 2 π ) , a następnie przesuwamy cząsteczkę o wybraną długość pod wylosowanym kątem.

Zaimplementujmy algorytm, który pozwoli utworzyć regularne trasy ruchu cząsteczki. Przyjmujemy następujące założenia:

  • cząsteczka na początku znajduje się w środku układu współrzędnych start_xy = (0, 0),

  • każdy ruch cząsteczki ma taką samą długość o wartości wektor = 1,

  • kierunek ruchu wyznaczony jest kątem fi losowanym tylko spośród wielokrotności kąta podstawowego,

  • wielokrotność kąta podstawowego będzie wyliczana na podstawie przyjętej wartości liczby całkowitej k,

  • kolejne współrzędne (x, y) cząstki będą wyliczane według wzorów:

X n = X n 1 + ( w e k t o r cos ( φ ) )
Y n = Y n 1 + ( w e k t o r sin ( φ ) )

Wielokrotność kąta podstawowego 2 π k uzyskamy w następujący sposób:

  • wylosujemy liczbę naturalną z przedziału <0; k) (nie włączając k),

  • podzielimy kąt przez k i pomnożymy przez wylosowaną liczbę.

Przykład 1

Napiszemy program, który wylicza losowe współrzędne poruszającej się cząsteczki i obrazuje jej ruch.

Współrzędne kolejnych punktów wygeneruje funkcja ruchy_browna(). Aby uzyskać takie same zestawy liczb pseudolosowych, użyjemy funkcji seed()seed(x)seed() z modułu random, która zainicjuje wartość początkową generatora liczb argumentem start_seed. Funkcja randint() pozwoli losować wartości całkowite z zadanego przedziału. Losowy kąt ruchu o wielkości 360 / k, dostosowany do uzyskania regularnych tras, będzie wyliczany w instrukcji: fi:float = randint(0, k‑1) * 2 * pi / k.

Wyliczone współrzędne kolejnych punktów, w których znajdzie się cząsteczka, zapisywać będziemy w dwóch listach (lista_x, lista_y) o tej samej długości. Funkcja ruchy_browna zwróci te listy, abyśmy mogli wykreślić ruch cząsteczki.

Do pokazania ruchu cząsteczki posłużymy się funkcją wizualizuj_ruchy_browna() i biblioteką matplotlibmatplotlibmatplotlib, która pozwala m.in. na rysowanie punktów na płaszczyźnie przy użyciu list zawierających ich współrzędne (x, y). Dzięki zastosowaniu odpowiedniej opcji punkty te zostaną automatycznie połączone liniami.

Linia 1. start podkreślnik xy znak równości otwórz nawias okrągły 0 przecinek 0 zamknij nawias okrągły kratka początkowe położenie cząsteczki. Linia 2. n znak równości 150 kratka liczba punktów. Linia 3. k znak równości 3 kratka obrót nastąpi o wielokrotność kąta 360 prawy ukośnik k stopni. Linia 4. wektor znak równości 1 kropka 0 kratka długość przemieszczenia. Linia 5. start podkreślnik seed znak równości 234567 kratka wartość inicjująca generator liczb. Linia 7. def ruchy podkreślnik browna otwórz nawias okrągły start podkreślnik xy dwukropek tuple przecinek n dwukropek int przecinek wektor dwukropek float przecinek start podkreślnik seed dwukropek int znak równości 234567 zamknij nawias okrągły minus zamknij nawias ostrokątny tuple dwukropek. Linia 8. from math import sin przecinek cos przecinek pi. Linia 9. from random import randint przecinek seed. Linia 11. seed otwórz nawias okrągły start podkreślnik seed zamknij nawias okrągły. Linia 12. lista podkreślnik x dwukropek list znak równości otwórz nawias kwadratowy start podkreślnik xy otwórz nawias kwadratowy 0 zamknij nawias kwadratowy zamknij nawias kwadratowy. Linia 13. lista podkreślnik y dwukropek list znak równości otwórz nawias kwadratowy start podkreślnik xy otwórz nawias kwadratowy 1 zamknij nawias kwadratowy zamknij nawias kwadratowy. Linia 15. for krok in range otwórz nawias okrągły 1 przecinek n zamknij nawias okrągły dwukropek. Linia 16. fi dwukropek float znak równości randint otwórz nawias okrągły 0 przecinek k minus 1 zamknij nawias okrągły asterysk 2 asterysk pi prawy ukośnik k kratka losowa wielokrotność kąta 360 prawy ukośnik k kąt. Linia 17. kratka wyliczenie nowych współrzędnych. Linia 18. dx znak równości lista podkreślnik x otwórz nawias kwadratowy krok minus 1 zamknij nawias kwadratowy plus wektor asterysk cos otwórz nawias okrągły fi zamknij nawias okrągły. Linia 19. dy znak równości lista podkreślnik y otwórz nawias kwadratowy krok minus 1 zamknij nawias kwadratowy plus wektor asterysk sin otwórz nawias okrągły fi zamknij nawias okrągły. Linia 20. kratka zapisanie nowego położenia w listach. Linia 21. lista podkreślnik x kropka append otwórz nawias okrągły dx zamknij nawias okrągły. Linia 22. lista podkreślnik y kropka append otwórz nawias okrągły dy zamknij nawias okrągły. Linia 24. return otwórz nawias okrągły lista podkreślnik x przecinek lista podkreślnik y zamknij nawias okrągły. Linia 27. def wizualizuj podkreślnik ruchy podkreślnik browna otwórz nawias okrągły wspolrzedne dwukropek tuple zamknij nawias okrągły minus zamknij nawias ostrokątny None dwukropek. Linia 28. import matplotlib kropka pyplot as plt. Linia 30. X dwukropek list znak równości wspolrzedne otwórz nawias kwadratowy 0 zamknij nawias kwadratowy. Linia 31. Y dwukropek list znak równości wspolrzedne otwórz nawias kwadratowy 1 zamknij nawias kwadratowy. Linia 32. plt kropka plot otwórz nawias okrągły X przecinek Y przecinek apostrof o minus apostrof przecinek color znak równości apostrof red apostrof przecinek linewidth znak równości 1 zamknij nawias okrągły. Linia 33. plt kropka xlabel otwórz nawias okrągły apostrof Współrzędne X apostrof zamknij nawias okrągły. Linia 34. plt kropka ylabel otwórz nawias okrągły apostrof Współrzędne Y apostrof zamknij nawias okrągły. Linia 35. plt kropka title otwórz nawias okrągły apostrof Ruchy Browna apostrof zamknij nawias okrągły. Linia 36. plt kropka grid otwórz nawias okrągły True zamknij nawias okrągły. Linia 37. plt kropka show otwórz nawias okrągły zamknij nawias okrągły. Linia 39. kratka wywołanie funkcji generujących współrzędne punktów. Linia 40. punkty znak równości ruchy podkreślnik browna otwórz nawias okrągły start podkreślnik xy przecinek n przecinek wektor przecinek start podkreślnik seed zamknij nawias okrągły. Linia 42. kratka wypisanie wyliczonych współrzędnych. Linia 43. print otwórz nawias okrągły punkty zamknij nawias okrągły. Linia 45. kratka wywołanie funkcji rysującej punkty. Linia 46. wizualizuj podkreślnik ruchy podkreślnik browna otwórz nawias okrągły punkty zamknij nawias okrągły.
Dla zainteresowanych

W definicjach funkcji po liście parametrów umieściliśmy kod -> tuple oraz -> None. To przykład użycia type hintstype hintstype hints, które omówiono w dokumencie PEP‑484 języka Python. Sposób ten – opcjonalny – umożliwia określenie oczekiwanych typów danych, co jest pomocne podczas przeglądania i modyfikowania kodu. W omawianym wypadku wskazujemy, że pierwsza funkcja zwraca tuplę (krotkę), natomiast druga nie zwraca niczego.

Symulacje utworzone za pomocą powyższego skryptu będą wyglądały następująco:

  • dla k = 3

RxBZkGgNOK02X
  • dla k = 4

R1UKJSB25fOQ2
  • dla k = 5

R1aTBo45Hq1Nr

Wzrost parametru k zwiększa losowość kierunku przemieszczania się cząsteczki. Wynika to z faktu, że im większa wartość k, tym więcej możliwych wielokrotności kąta wyznaczających kierunek ruchu. Dążąc do jak największej losowości, możemy zwiększać parametr k aż do wartości 360.

Dla k = 360 uzyskamy następujący wykres:

RKOJyCnxnzlsT

W przedstawiony sposób możemy symulować ruchy Browna dla różnych danych początkowych, tzn. liczby kroków, liczby wielokrotności oraz długości przemieszczenia.

1
Polecenie 1

Bazując na poprzednich przykładach, spróbuj wywołać program w taki sposób, aby uzyskać na wykresie tylko punkty początkowy i końcowy podczas wykonywania 300 kolejnych kroków. Jako wektor przyjmij wartość 0.73, a start_seed niech wynosi 432432. Dobierz eksperymentalnie wartość dla parametru k.

Symulacja ruchów Browna metodą Monte Carlo z wykorzystaniem procesu Wienera

Przypomnijmy uproszczoną definicję procesu Wienera:

W ( 0 ) = 0
W ( t + d t ) = W ( t ) + k N ( 0 ,   d t )

gdzie N ( μ , σ 2 ) jest zmienną losową zgodną z rozkładem normalnym Gaussarozkład normalny Gaussarozkładem normalnym Gaussa, z parametrami: wartość oczekiwana  μ oraz odchylenie standardowe  σ , przy czym k jest współczynnikiem długości wykonywanych kroków.

Jak możemy wywnioskować z tej definicji, jest to proces uzależniony nie od liczby iteracji, a od czasu. Odległość cząsteczki od położenia początkowego po wykonaniu stu kroków z parametrem d t = 0.1 jest uzależniona od prawdopodobieństwa w taki sam sposób, co odległość po wykonaniu jednego kroku z parametrem d t = 10. Innymi słowy, wynik końcowy stu kroków z parametrem d t = 0.1 jest statystycznie równy wynikowi końcowemu jednego kroku z parametrem d t = 10.

Aby stworzyć symulację wykorzystującą proces Wienera, musimy zastanowić się nad jej warunkami początkowymi. Zaczynamy od umieszczenia cząsteczki w punkcie ( 0, 0 ) . Przed narysowaniem kolejnej pozycji cząsteczki będziemy czekali 0,05 sekundy, czyli przyrost czasu wyniesie dt = 0.05. Dla tak dobranego parametru dt ustawimy współczynnik k = 0.2, co wpłynie na skalę przemieszczenia.

W języku Python, aby uzyskać losową liczbę zgodną z rozkładem normalnym, wykorzystamy moduł random i zawartą w nim funkcję gauss(mu, sigma). Pierwszy parametr oznacza wartość oczekiwaną μ, drugi odchylenie standardowe σ. W omawianym przykładzie przyjmujemy wartość oczekiwaną równą 0 i odchylenie równe 1.

W ten sposób otrzymujemy liczbę N(0,1). Jak możemy uzyskać N ( 0 ,   d t ) ? Skorzystamy z własności wynikającej ze wzoru funkcji rozkładu normalnego, mianowicie:

N ( μ , σ 2 ) = 1 σ N ( 0 ,   1 )   +   μ

Do wyznaczenia N ( 0 ,   d t ) możemy więc użyć kodu:

Linia 1. from random import gauss. Linia 2. from math import sqrt. Linia 3. kratka N otwórz nawias okrągły 0 przecinek dt zamknij nawias okrągły. Linia 4. N podkreślnik 0 podkreślnik dt znak równości 1 prawy ukośnik sqrt otwórz nawias okrągły dt zamknij nawias okrągły asterysk gauss otwórz nawias okrągły 0 przecinek 1 zamknij nawias okrągły.

Wykorzystując tę metodę, kolejne położenia cząsteczek będziemy wyznaczać następująco:

Linia 1. x plus znak równości k asterysk sqrt otwórz nawias okrągły dt zamknij nawias okrągły asterysk gauss otwórz nawias okrągły 0 przecinek 1 zamknij nawias okrągły. Linia 2. y plus znak równości k asterysk sqrt otwórz nawias okrągły dt zamknij nawias okrągły asterysk gauss otwórz nawias okrągły 0 przecinek 1 zamknij nawias okrągły.

Operację tę będziemy powtarzać do momentu, kiedy cząsteczka wykona n kroków. Dodatkowo po każdym wyliczeniu położenia cząsteczki uwzględnimy przyrost czasu wskazany w parametrze dt, wstrzymując na ten czas wykonywanie programu. Parametr dt wykorzystamy również do opóźnienia rysowania ruchu cząsteczki. Pełny kod skryptu zawierający zmienioną funkcję rysującą animuj_ruchy_browna() zamieszczamy poniżej:

Linia 1. start podkreślnik xy znak równości otwórz nawias okrągły 0 przecinek 0 zamknij nawias okrągły kratka początkowe położenie cząsteczki. Linia 2. n znak równości 100 kratka liczba kroków. Linia 3. k znak równości 0 kropka 2 kratka współczynnik długości kroków. Linia 4. dt znak równości 0 kropka 05 kratka przyrost czasu. Linia 6. def ruchy podkreślnik browna otwórz nawias okrągły zamknij nawias okrągły minus zamknij nawias ostrokątny None dwukropek. Linia 7. from math import sqrt. Linia 8. from random import gauss. Linia 9. from time import sleep. Linia 11. lista podkreślnik x dwukropek list znak równości otwórz nawias kwadratowy start podkreślnik xy otwórz nawias kwadratowy 0 zamknij nawias kwadratowy zamknij nawias kwadratowy. Linia 12. lista podkreślnik y dwukropek list znak równości otwórz nawias kwadratowy start podkreślnik xy otwórz nawias kwadratowy 1 zamknij nawias kwadratowy zamknij nawias kwadratowy. Linia 14. for krok in range otwórz nawias okrągły 1 przecinek n zamknij nawias okrągły dwukropek. Linia 15. lista podkreślnik x kropka append otwórz nawias okrągły lista podkreślnik x otwórz nawias kwadratowy krok minus 1 zamknij nawias kwadratowy plus k asterysk sqrt otwórz nawias okrągły dt zamknij nawias okrągły asterysk gauss otwórz nawias okrągły 0 przecinek 1 zamknij nawias okrągły zamknij nawias okrągły. Linia 16. lista podkreślnik y kropka append otwórz nawias okrągły lista podkreślnik y otwórz nawias kwadratowy krok minus 1 zamknij nawias kwadratowy plus k asterysk sqrt otwórz nawias okrągły dt zamknij nawias okrągły asterysk gauss otwórz nawias okrągły 0 przecinek 1 zamknij nawias okrągły zamknij nawias okrągły. Linia 17. sleep otwórz nawias okrągły dt zamknij nawias okrągły. Linia 19. return otwórz nawias okrągły lista podkreślnik x przecinek lista podkreślnik y zamknij nawias okrągły. Linia 22. def animuj podkreślnik ruchy podkreślnik browna otwórz nawias okrągły wspolrzedne dwukropek tuple zamknij nawias okrągły minus zamknij nawias ostrokątny None dwukropek. Linia 23. import matplotlib kropka pyplot as plt. Linia 24. from matplotlib kropka animation import FuncAnimation. Linia 26. X dwukropek list znak równości wspolrzedne otwórz nawias kwadratowy 0 zamknij nawias kwadratowy. Linia 27. Y dwukropek list znak równości wspolrzedne otwórz nawias kwadratowy 1 zamknij nawias kwadratowy. Linia 28. fig przecinek ax znak równości plt kropka subplots otwórz nawias okrągły zamknij nawias okrągły. Linia 29. ax kropka figure kropka set podkreślnik size podkreślnik inches otwórz nawias okrągły 7 przecinek 7 zamknij nawias okrągły. Linia 30. ax kropka set podkreślnik xlabel otwórz nawias okrągły apostrof Współrzędne X apostrof zamknij nawias okrągły. Linia 31. ax kropka set podkreślnik ylabel otwórz nawias okrągły apostrof Współrzędne Y apostrof zamknij nawias okrągły. Linia 32. ax kropka set podkreślnik title otwórz nawias okrągły apostrof Ruchy Browna apostrof zamknij nawias okrągły. Linia 34. line przecinek znak równości ax kropka plot otwórz nawias okrągły X przecinek Y przecinek color znak równości apostrof r apostrof przecinek marker znak równości apostrof o apostrof zamknij nawias okrągły. Linia 36. def animate otwórz nawias okrągły i zamknij nawias okrągły dwukropek. Linia 37. line kropka set podkreślnik data otwórz nawias okrągły X otwórz nawias kwadratowy dwukropek i plus 1 zamknij nawias kwadratowy przecinek Y otwórz nawias kwadratowy dwukropek i plus 1 zamknij nawias kwadratowy zamknij nawias okrągły. Linia 38. return line. Linia 40. ani znak równości FuncAnimation otwórz nawias okrągły fig przecinek animate przecinek interval znak równości dt asterysk 1000 przecinek repeat znak równości False zamknij nawias okrągły. Linia 41. plt kropka show otwórz nawias okrągły zamknij nawias okrągły. Linia 43. punkty znak równości ruchy podkreślnik browna otwórz nawias okrągły zamknij nawias okrągły. Linia 44. animuj podkreślnik ruchy podkreślnik browna otwórz nawias okrągły punkty zamknij nawias okrągły.

Prześledźmy przykładowe wyniki symulacji po uruchomieniu przedstawionego kodu z wybranymi parametrami.

Końcowy wynik symulacji ruchu cząsteczki dla parametrów n = 100, dt = 0.05, k = 0.2:

R1ZAoJw2igPUb

Wizualizacja ruchu cząsteczki dla parametrów n = 10, dt = 0.5, k = 0.2:

R1GmcFZP95QaR

Rezultat symulacji przeprowadzonej dla danych n = 1000, dt = 0.005, k = 0.2:

RuHBPoFOzF4T7

Zwróć uwagę, że przedstawione wyniki nie różnią się czasem symulacji, a jedynie jej dokładnością.

Już wiesz

Podsumujmy najważniejsze wiadomości:

  • metoda Monte Carlo może być przydatna przy szacowaniu wyników na podstawie dużej liczby zmiennych losowych,

  • w języku Python możemy stosować tak zwane type hintstype hintstype hints,

  • funkcja seed() pozwala ustawić wartość początkową (ziarno) dla generatora liczb pseudolosowych.

Słownik

importowanie
importowanie

operacja pozwalająca w programie (funkcji) używać dodatkowych funkcji lub obiektów, dostępnych w zewnętrznych bibliotekach

matplotlib
matplotlib

biblioteka służąca do przedstawienia obrazów złożonych z punktów o współrzędnych x oraz y (np. wykresów, histogramów, rozkładów); moduł matplotlib nie jest dostępny w standardowej instalacji języka Python; należy go zainstalować, korzystając z mechanizmu pip

rozkład normalny Gaussa
rozkład normalny Gaussa

ciągły rozkład prawdopodobieństwa, zdefiniowany następującą funkcją gęstości prawdopodobieństwa:

f μ , σ ( x ) = 1 σ 2 π exp ( ( x μ ) 2 2 σ 2 )

gdzie μ to wartość oczekiwana zmiennej losowej x, a σ to odchylenie standardowe zmiennej losowej x

seed(x)
seed(x)

(ang. seed – ziarno) funkcja pozwalająca ustalić punkt początkowy obliczeń, których rezultatem jest wygenerowanie liczby pseudolosowej przez funkcje takie jak random(), randint() i podobne; po określeniu wartości zmiennej x wyniki działania funkcji losowych są powtarzalne; jeśli nie podano wartości x lub ustalono ją jako None, za wartość x przyjmuje się aktualny czas

type hints
type hints

informacja o typie danego obiektu zapisana literalnie w kodzie programu; z uwagi na dynamiczną naturę języka Python pomaga w zrozumieniu kodu programu i pozwala na sprawdzanie kodu przez różne środowiska IDE (np. PyCharm, SublimeText lub Atom); przykłady użycia: zmienna: int = 10; type hints są dokładnie opisane w dokumencie PEP 484 – ich inna nazwa to annotations (opisane w dokumencie PEP 3107)

zmienna losowa
zmienna losowa

zmienna, której wartość pozostaje nieznana do momentu losowania; w języku Python można uzyskać ją między innymi dzięki funkcjom należącym do modułu random:

  • randint(a, b) – zwraca liczbę całkowitą z przedziału <a, b>

  • randrange(a, b, c) – zwraca liczbę całkowitą z przedziału <a, b‑1>, z krokiem c

  • random() – zwraca liczbę rzeczywistą z przedziału <0.0, 1.0)