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\pi )$, 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:

Wielokrotność kąta podstawowego $\frac{2\pi }{k}$ uzyskamy w następujący sposób:

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

• podzielimy kąt $2\pi$ przez k i pomnożymy przez wylosowaną liczbę.

start_xy = (0, 0)  # początkowe położenie cząsteczki
n = 150  # liczba punktów
k = 3  # obrót nastąpi o wielokrotność kąta 360/k stopni
wektor = 1.0  # długość przemieszczenia
start_seed = 234567  # wartość inicjująca generator liczb

def ruchy_browna(start_xy: tuple, n: int, wektor: float, start_seed: int = 234567) -> tuple:
    from math import sin, cos, pi
    from random import randint, seed

    seed(start_seed)
    lista_x: list = [start_xy[0]]
    lista_y: list = [start_xy[1]]

    for krok in range(1, n):
        fi:float = randint(0, k-1) * 2 * pi / k  # losowa wielokrotność kąta 360/k kąt
        # wyliczenie nowych współrzędnych
        dx = lista_x[krok - 1] + wektor * cos(fi)
        dy = lista_y[krok - 1] + wektor * sin(fi)
        # zapisanie nowego położenia w listach
        lista_x.append(dx)
        lista_y.append(dy)

    return (lista_x, lista_y)


def wizualizuj_ruchy_browna(wspolrzedne: tuple) -> None:
    import matplotlib.pyplot as plt

    X: list = wspolrzedne[0]
    Y: list = wspolrzedne[1]
    plt.plot(X, Y, 'o-', color='red', linewidth=1)
    plt.xlabel('Współrzędne X')
    plt.ylabel('Współrzędne Y')
    plt.title('Ruchy Browna')
    plt.grid(True)
    plt.show()

# wywołanie funkcji generujących współrzędne punktów
punkty = ruchy_browna(start_xy, n, wektor, start_seed)

# wypisanie wyliczonych współrzędnych
print(punkty)

# wywołanie funkcji rysującej punkty
wizualizuj_ruchy_browna(punkty)

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

• dla k = 3

RxBZkGgNOK02X

Symulacja ruchów Browna w języku Python – przykład 1.

• dla k = 4

R1UKJSB25fOQ2

Symulacja ruchów Browna w języku Python – przykład 2.

• dla k = 5

R1aTBo45Hq1Nr

Symulacja ruchów Browna w języku Python – przykład 3.

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

Symulacja ruchów Browna w języku Python – przykład 4.

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.

Pokaż rozwiązanie

Linia 1. start podkreślnik xy znak równości otwórz nawias okrągły 0 przecinek 0 zamknij nawias okrągły. Linia 2. n znak równości 300. Linia 3. k znak równości 360. Linia 4. wektor znak równości 0 kropka 73. Linia 5. start podkreślnik seed znak równości 432432. Linia 7. 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 8. X znak równości otwórz nawias kwadratowy punkty otwórz nawias kwadratowy 0 zamknij nawias kwadratowy otwórz nawias kwadratowy 0 zamknij nawias kwadratowy przecinek punkty otwórz nawias kwadratowy 0 zamknij nawias kwadratowy otwórz nawias kwadratowy minus 1 zamknij nawias kwadratowy zamknij nawias kwadratowy. Linia 9. Y znak równości otwórz nawias kwadratowy punkty otwórz nawias kwadratowy 1 zamknij nawias kwadratowy otwórz nawias kwadratowy 0 zamknij nawias kwadratowy przecinek punkty otwórz nawias kwadratowy 1 zamknij nawias kwadratowy otwórz nawias kwadratowy minus 1 zamknij nawias kwadratowy zamknij nawias kwadratowy. Linia 10. wizualizuj podkreślnik ruchy podkreślnik browna otwórz nawias okrągły otwórz nawias okrągły X przecinek Y zamknij nawias okrągły zamknij nawias okrągły. Linia 12. kratka wypisanie współrzędnych punktu końcowego. Linia 13. print otwórz nawias okrągły cudzysłów x znak równości cudzysłów przecinek X otwórz nawias kwadratowy 1 zamknij nawias kwadratowy przecinek cudzysłów y znak równości cudzysłów przecinek Y otwórz nawias kwadratowy 1 zamknij nawias kwadratowy zamknij nawias okrągły. Linia 14. kratka x znak równości 3 kropka 7971201935166907 przecinek y znak równości 14 kropka 897028925526966.

start_xy = (0, 0)
n = 300
k = 360
wektor = 0.73
start_seed = 432432

punkty = ruchy_browna(start_xy, n, wektor, start_seed)
X = [punkty[0][0], punkty[0][-1]]
Y = [punkty[1][0], punkty[1][-1]]
wizualizuj_ruchy_browna((X, Y))

# wypisanie współrzędnych punktu końcowego
print("x = ",X[1]," y = ",Y[1])
# x = 3.7971201935166907, y = 14.897028925526966

R1CK6ecBeF3cW

Symulacja ruchów Browna w języku Python – punkt początkowy i końcowy.

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

Przypomnijmy uproszczoną definicję procesu Wienera:

gdzie $N(\mu ,{\sigma }^2)$ jest zmienną losową zgodną z rozkładem normalnym Gaussarozkład normalny Gaussarozkładem normalnym Gaussa, z parametrami: wartość oczekiwana $\mu$ oraz odchylenie standardowe $\sigma$, 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 $dt=0.1$ jest uzależniona od prawdopodobieństwa w taki sam sposób, co odległość po wykonaniu jednego kroku z parametrem $dt=10$. Innymi słowy, wynik końcowy stu kroków z parametrem $dt=0.1$ jest statystycznie równy wynikowi końcowemu jednego kroku z parametrem $dt=10$.

Aby stworzyć symulację wykorzystującą proces Wienera, musimy zastanowić się nad jej warunkami początkowymi. Zaczynamy od umieszczenia cząsteczki w punkcie $\left(\mathrm{0, 0}\right)$. 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ą $\mu$, drugi odchylenie standardowe $\sigma$. 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,dt)$? Skorzystamy z własności wynikającej ze wzoru funkcji rozkładu normalnego, mianowicie:

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

from random import gauss
from math import sqrt
# N(0, dt)
N_0_dt = 1 / sqrt(dt) * gauss(0, 1)

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

x += k * sqrt(dt) * gauss(0, 1)
y += k * sqrt(dt) * gauss(0, 1)

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:

start_xy = (0, 0)  # początkowe położenie cząsteczki
n = 100  # liczba kroków
k = 0.2  # współczynnik długości kroków
dt = 0.05  # przyrost czasu

def ruchy_browna() -> None:
    from math import sqrt
    from random import gauss
    from time import sleep

    lista_x: list = [start_xy[0]]
    lista_y: list = [start_xy[1]]

    for krok in range(1, n):
        lista_x.append(lista_x[krok - 1] + k * sqrt(dt) * gauss(0, 1))
        lista_y.append(lista_y[krok - 1] + k * sqrt(dt) * gauss(0, 1))
        sleep(dt)

    return (lista_x, lista_y)


def animuj_ruchy_browna(wspolrzedne: tuple) -> None:
    import matplotlib.pyplot as plt
    from matplotlib.animation import FuncAnimation

    X: list = wspolrzedne[0]
    Y: list = wspolrzedne[1]
    fig, ax = plt.subplots()
    ax.figure.set_size_inches(7, 7)
    ax.set_xlabel('Współrzędne X')
    ax.set_ylabel('Współrzędne Y')
    ax.set_title('Ruchy Browna')

    line, = ax.plot(X, Y, color='r', marker='o')

    def animate(i):
        line.set_data(X[:i+1], Y[:i+1])
        return line

    ani = FuncAnimation(fig, animate, interval=dt*1000, repeat=False)
    plt.show()

punkty = ruchy_browna()
animuj_ruchy_browna(punkty)

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

Zrzut ekranu przedstawia program: Symulacja ruchów Browna dla n = 100, dt=0,05, k=0,2.

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

R1GmcFZP95QaR

Zrzut ekranu przedstawia program: Symulacja ruchów Browna dla n = 10, dt=0,5, k=0,2.

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

RuHBPoFOzF4T7

Zrzut ekranu przedstawia program: Symulacja ruchów Browna dla n = 1000, dt=0,005, k=0,2.

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

Słownik