Modelowanie ruchów Browna

Ruchy Browna (o których możesz przeczytać m.in. w e‑materiale Modelowanie ruchów BrownaP14JSHgfzModelowanie ruchów Browna) zostały opisane — niezależnie — przez Alberta Einsteina, w $1905$ roku oraz przez polskiego fizyka Mariana Smoluchowskiego, w $1906$ roku. Ich nazwa pochodzi od botanika Roberta Browna, który pod mikroskopem zaobserwował chaotyczne ruchy pyłku kwiatowego. Ruchy Browna powodowane są nieustannym zderzaniem się cząsteczek płynu w danym ośrodku. Przykładowy model takiego procesu może polegać na umieszczeniu wirtualnej cząsteczki w wirtualnym płynie, w którym — co pewien określony czas — cząsteczka zostaje przemieszczona w losowym kierunku, w wyniku wirtualnej kolizji.

Cały proces ma charakter uśredniony. Dzięki temu każda kolizja powoduje przesunięcie cząsteczki o określoną z góry odległość.

Metoda Monte Carlo

Metoda Monte Carlo została opracowana między innymi przez polskiego naukowca Stanisława Ulama oraz węgiersko — amerykańskiego naukowca Johna von Neumanna. Wykorzystano ją podczas prac nad bombą jądrową w celu symulacji procesu rozpadu jąder cząsteczek. Proces był na tyle skomplikowany, że zastosowanie klasycznych modeli obliczeniowych nie zdawało egzaminu.

Metoda polega na przeprowadzeniu procesu w losowy sposób, zakładając, że po wystarczająco dużej liczbie prób otrzymamy wynik zbliżony do rzeczywistego.

Symulacja ruchów Browna metodą Monte Carlo

Symulacja ruchów Browna metodą Monte Carlo wykorzystywana jest zazwyczaj w celu poznania możliwej trasy, zgodnie z którą będzie poruszać się cząsteczka. Tego rodzaju symulacje mogą generować sztuczne wykresy giełdowe, dźwięki (np. szum czerwonyszum czerwonyszum czerwony), trasy losowego błądzenia oraz różnego rodzaju kształty (np. płatki śniegu).

W celu utworzenia szkicu trasy, którą poruszała się cząsteczka, możemy wykorzystać graficzne możliwości języków programowania. W tej części lekcji omówimy jednak najważniejsze aspekty symulacji ruchów Browna metodą Monte Carlo.

Aby zwizualizować proces losowego poruszania się cząsteczki w płynie, zakładamy, że każda kolizja powoduje przemieszczenie się cząsteczki o określoną długość w dowolnym kierunku. W tym celu wybieramy losowy kąt z zakresu $⟨0,2\pi )$, a następnie dokonujemy jej przesunięcia o wybraną długość pod wylosowanym kątem.

W poprzednich e—materiałach dotyczących modelowania ruchów Browna posługiwaliśmy się algorytmem, w którym cząsteczka reprezentowana była przez dwie zmienne zmiennoprzecinkowe. W każdej kolizji przesuwaliśmy cząsteczkę o wektor długości r pod losowo wybranym kątem. Pozycję po każdym przesunięciu zapisywaliśmy w tablicy. Po zakończeniu algorytmu wypisywaliśmy ostatni punkt.

Jeżeli chcemy utworzyć regularniejsze trasy, to możemy pozwolić sobie na odstępstwo od wyboru każdego losowego kąta z przedziału $⟨0,2\pi )$. Możemy na przykład postanowić, że losujemy tylko spośród wielokrotności wybranego kąta. W tym celu należy wylosować liczbę całkowitą z przedziału od 0 (włącznie) do wybranej liczby k (nie włączając k). Wówczas, po podzieleniu kąta $2\pi$ przez k, a następnie pomnożeniu przez wylosowaną liczbę, uzyskamy losową wielokrotność kąta $\frac{2\pi }{k}$.

Przykładem pseudokodu uwzględniającego powyższą zmianę (random_int(k)) jest funkcja generująca losową liczbę całkowitą od 0 (włącznie) do k (włącznie):

// liczba kolizji
n = 10000

// obrót nastąpi o wielokrotność kąt 360/k stopni
k = 3

// długość wektora
r = 0.1

// początkowe położenie cząsteczki
x = 0.0
y = 0.0

// tablice przechowujące kolejne pozycje cząsteczki
deklaracja tablicy ruchy_x o wielkości n+1
deklaracja tablicy ruchy_y o wielkości n+1

// zapisz w tablicach początkowe położenie cząsteczki
ruchy_x[0] = x
ruchy_y[0] = y

// dla każdej kolizji
dla i = 1, 2, ..., n wykonuj
	// losowa wielokrotność kąta 360/k kąt
	fi = random_int(k - 1) * 2 * PI / k;

	// zamień postać wektora (długość, kąt) na współrzędne
	dx = r * cosinus(fi)
    dy = r * sinus(fi)
    
    // przesuń cząsteczkę o wektor [dx, dy]
    x = x + dx
    y = y + dy
    
    // zapisz w tablicach obecne położenie cząsteczki
	ruchy_x[i] = x
	ruchy_y[i] = y

// wypisz końcową pozycję
wypisz x
wypisz y

Utworzone w ten sposób trasy mogą wyglądać następująco:

k = 3

RJPw6FbdydsiX

Ilustracja przedstawia trasę dla n = 3. Tworzona jest ona z trójkątów.

k = 4

RXgsB5a8WGu3D

Ilustracja przedstawia trasę dla n = 4. Tworzona jest ona z prostokątów.

n = 5

RRO7pnIOdwbSP

Ilustracja przedstawia trasę dla n = 4. Tworzona jest ona z nieregularnych kształtów.

Opisane powyżej zagadnienie znane jest także jako błądzenie losowe. Spróbujmy przeprowadzić prosty eksperyment za pomocą metody Monte Carlo. Umieśćmy cząsteczkę na osi liczbowej w punkcie 0. W każdej kolejnej iteracji niech cząsteczka porusza się losowo w prawo lub lewo o 1. Po k (np. 20) iteracjach zapiszmy pozycję cząsteczki. Eksperyment powtórzmy n (np. 10 000) razy. Na koniec stwórzmy wykres słupkowy, w którym na osi poziomej znajdzie się odległość cząstki od położenia początkowego, a na osi pionowej liczba prób, które zakończyły się na tej pozycji.

Eksperyment możemy zrealizować zgodnie z poniższym pseudokodem:

wynik[0..2 * k] - tablica wypełniona zerami

// liczba prób
n = 10000
// liczba ruchów
k = 20

dla i = 0, 1, ..., n - 1 wykonuj
	// pozycja na środku
    pozycja = k
    
    dla j = 0, 1, ..., k - 1 wykonuj
        pozycja += random_int(1) * 2 - 1
        wynik[pozycja] += 1
        
rysuj wykres tablicy wynik

Wynik przykładowego eksperymentu:

Rr41KrKjrT4jx

Ilustracja przedstawia wykres słupkowy: Wynik eksperymentu dla 10000 prób. Oś Y opisana jest jako: Liczba zakończona prób na danej pozycji, a oś X jako: Pozycja od punktu początkowego. Wartości rosną od punktu minus 20 do punktu 0, następnie maleją do punktu 20.

Przeprowadzony eksperyment to nic innego jak puszczone kulki na desce Galtonadeska Galtonadesce Galtona. Kulka zaczyna swoją wędrówkę od góry, a następnie natrafia na przeszkodę, która powoduje jej ruch w lewo lub w prawo z jednakowym prawdopodobieństwem.

RWCMyyYhweLuD

Ilustracja przedstawia puszczone kulki na desce Galtona. Po przejściu przeszkód w kształcie czerwonych kół, kula wpada do jednego z koszu umieszczonych pod przeszkodami.

Eksperyment wizualizuje rozkład dwumianowyrozkład dwumianowyrozkład dwumianowy (w Polsce zwany także rozkładem Bernoulliego), który w nieskończoności dąży do rozkładu normalnegorozkład normalnyrozkładu normalnego, znanego jako rozkład Gaussa.

Proces Wienera

Przeprowadzony eksperyment doprowadza nas do procesu Wienera, który w gruncie rzeczy jest dokładnie tym samym, co ruchy Browna, czyli matematycznym opisem chaotycznego ruchu cząsteczki w płynie. Za pomocą tego procesu oraz metody Monte Carlo możemy przeprowadzić symulację ruchów Browna, która dużo lepiej odzwierciedla poruszanie się cząsteczki w naturze. Proces ten w celu wyznaczenia następnego położenia cząsteczki wykorzystuje funkcję rozkładu prawdopodobieństwa Gaussa, który określamy wzorem:

gdzie $\mu$ oznacza wartość oczekiwaną (średnią), a $\sigma$ to odchylenie standardowe.

Przykładowe gaussianygaussiangaussiany dla danych parametrów $\mu$ oraz $\sigma$ wyglądają następująco:

RiYLUdgwrYQ0z

Ilustracja przedstawia wykres z opisanymi prostymi: mikro=0, sigma=0,2 kolor niebieski; mikro=0, sigma=0,5 kolor czerwony; mikro=0, sigma=1 żółty, mikro=2, sigma=0,25 Zielony. Największe wartości przyjmuje linia niebieska, następnie , zielona, czerwona oraz żółta.

Zwróćmy uwagę na podobieństwo rozkładów prawdopodobieństwa do poprzednio utworzonego wykresu.

Proces Wienera, oznaczony symbolem $W$, zdefiniowany jest następująco:

1. Proces $W$ jest krzywą ciągłą o punkcie początkowym równym 0.

   $$W\left(0\right)=0$$

2. Proces $W$ ma niezależne od siebie przyrosty (lub spadki).

3. Rozkład prawdopodobieństwa przyrostów w czasie o długości $dt$ jest rozkładem normalnym Gaussa o wartości średniej 0 z wariancją równą $dt$.

   $$W(t+dt)-W\left(t\right)~N(0, dt)$$

Definicję tę można zatem przedstawić w uproszczeniu jako funkcję błądzenia losowego, w zależności od czasu, w następujący sposób:

gdzie $N(\mu ,{\sigma }^2)$ jest zmienną losową zgodną z rozkładem normalnym Gaussa z parametrami $\mu$ oraz $\sigma$. Przy czym $k$ jest współczynnikiem długości wykonywanych kroków.

Realizacja wyznaczania powyższego procesu powinna być zgodna z następującym pseudokodem (sqrt(d) to funkcja obliczająca pierwiastek danej liczby d):

// liczba kolizji
n = 100

// współczynnik długości kroku
k = 1

// przyrost czasu
dt = 0.1

// tablica przechowujące kolejne pozycje cząsteczki
deklaracja tablicy w o wielkości n+1

// początkowe położenie cząsteczki
w[0] = 0.0

dla i = 1, 2, ..., n wykonuj
    w[i] = w[i - 1] + k * random_n(0, sqrt(dt))

wykreśl trasę x, y

Wykres przykładowej trasy $W$ utworzonej tą metodą w zależności od czasu (n = 1000, k = 1, dt = 0.1):

R1WjKUO3lg86z

Ilustracja przedstawia wykres przykładowej trasy W. Na osi X są wartości od zera do stu, na osi Y od minus trzydziestu do plus czterdziestu. Krzywa wykresu ma trend rosnący. Rośnie od wartości zero na osiach do wartości około dziesięciu na osi Y i stu na osi X. Krzywa ma liczne odchylenia.

Powyższy kształt można zastosować jako losowo wygenerowany profil hipsometrycznyprofil hipsometrycznyprofil hipsometryczny lub wykres giełdowy.

Jeżeli chcemy wygenerować trasę dwuwymiarową, musimy jednocześnie utworzyć 2 trasy, które traktujemy osobno jako współrzędną poziomą oraz pionową.

// liczba kolizji
n = 100

// współczynnik długości kroku
k = 1

// przyrost czasu
dt = 0.1

// tablice przechowujące kolejne pozycje cząsteczki
deklaracja tablicy x o wielkości n+1
deklaracja tablicy y o wielkości n+1

// początkowe położenie cząsteczki
x[0] = 0.0
y[0] = 0.0

dla i = 1, 2, ..., n wykonuj
    x[i] = x[i - 1] + k * random_n(0, sqrt(dt))
    y[i] = y[i - 1] + k * random_n(0, sqrt(dt))

wykreśl trasę x, y

Wygenerowane w ten sposób trasy mogą wyglądać następująco:

RJtEArs4QsWwz

Ilustracja przedstawia trasę stworzoną z nieregularnych kształtów.

Zwróćmy uwagę na różnorodne długości, o które przemieszcza się cząsteczka. Symulacja wykorzystująca proces Wienera ma wyjątkowo przydatną własność. Aby poznać wynik symulacji, nie musimy jej przeprowadzać w całości. W momencie, gdy interesuje nas tylko wynik symulacji, wystarczy, że przeprowadzimy jej pierwszy krok z parametrem $dt$, który ustawimy na długość całej symulacji.

Słownik