Co to jest?

Monte Carlo, spacery i polimery

Z czym kojarzy się Monte Carlo? Z kasynami, hazardem, ruletką. A więc Monte Carlo jest symbolem potęgi przypadku, który jednym pozwala zbijać fortuny, innych rujnuje. Są jednak ludzie, którzy przypadkowość, losowość potrafią okiełznać i wykorzystać z pożytkiem (przynajmniej dla siebie). Oczywiście, takimi ludźmi są właściciele kasyn gry, ale nie tylko. Również matematycy zaprzęgają losowość do pożytecznej roboty. Metody obliczeniowe oparte na generowaniu (symulowaniu komputerowym) zdarzeń losowych nazywają się... oczywiście, metodami Monte Carlo.

Liczenie przedmiotów nie wydaje się zadaniem wymagającym wyrafinowanej matematyki, wykraczającej ponad program przedszkolny. Nie widać też, jaką rolę w czynności zliczania miałby odgrywać przypadek. A jednak! Rozważymy pozornie proste zadanie zliczania dróg, w którym ujawnią się nieoczekiwane trudności i pokażemy, jak pomysłowe metody Monte Carlo wynaleziono, aby te trudności pokonać.

Wyobraźmy sobie, że idziemy na spacer. Załóżmy, że mieszkamy w mieście, w którym ulice tworzą regularną, kwadratową siatkę, tak jak na rysunku 1. Ruszamy z punktu początkowego, wybierając jedną z czterech ulic. Idziemy zatem na północ, na zachód, na południe lub na wschód. Przechodzimy do sąsiedniego skrzyżowania i znowu wybieramy jeden z 4 kierunków. Powtarzamy całą procedurę, powiedzmy, $k$ razy. Przebyta droga jest łamaną składającą się z $k$ odcinków o jednostkowej długości. Na rysunku 1 mamy przykład drogi o długości $| k = 15.$ Używając układu współrzędnych "zaczepionego w punkcie początkowym", możemy tę drogę zapisać jako ciąg kolejno odwiedzanych punktów:

(0,0) (0,1) (1,1) (1,2) (0,2) ⋯ (2,0) (2,−1).

Inny sposób zakodowania tej samej drogi to podanie kierunków ruchu w kolejnych "krokach": $,E,N,W,W,...,E,S. |N$ Ile jest wszystkich dróg o długości $k?$ Odpowiedź jest łatwa: $|4k.$ Musimy tu wyjaśnić, że licząc drogi wychodzące z tego samego punktu, utożsamiamy drogi z ciągami literek z alfabetu $,W,S,E}. {N$

Zmieńmy teraz nasze "reguły spacerowania". Przypuśćmy, że nie chcemy powtórnie przechodzić przez punkt, który już wcześniej odwiedziliśmy. W ciągu opisującym drogę nie może dwukrotnie pojawić się ten sam punkt. Tę własność ma np. droga na rysunku 1. Po angielsku taka droga nazywa się Self Avoiding Walk, w skrócie SAW. Ile jest SAW-ów o długości $k?$ Wystarczy je wszystkie narysować i policzyć, prawda? Tą metodą można zbudować następującą tabelkę. Liczbę SAW-ów o długości $| k$ oznaczymy przez $| Zk.$ Ale co dalej? Im dalej, tym gorzej. Ile jest SAW-ów o długości $|k = 50?$ W istocie, liczenie SAW-ów jest na tyle trudnym zadaniem, że matematycy (i fizycy, o czym dalej) uciekają się do przybliżonego zliczania metodą Monte Carlo. Najprostszy algorytm tego rodzaju jest niemal oczywisty. Ustalmy $| k.$ Oznaczmy przez $|𝒳$ zbiór wszystkich dróg o długości $k,$ zaś przez $𝒵 ⊂ 𝒳$ zbiór SAW-ów o długości $k.$ Oczywiście $𝒵 = Zk$ i $| 𝒳 = 4k$ (kreseczki $|$ oznaczają liczbę elementów zbioru; dla uproszczenia przy symbolach zbiorów zrezygnowaliśmy z indeksu $| k$ ). Jeśli wybierzemy losowo jedną z dróg ze zbioru $𝒳,$ to prawdopodobieństwo otrzymania "przypadkiem" SAW-a jest równe $𝒵 / 𝒳 .$ Jeśli wybierzemy losowo i niezależnie $n$ dróg, to spośród nich znajdzie się pewna liczba, powiedzmy $| n1,$ SAW-ów. Widać, że dla dużych $| n$ powinniśmy mieć

n1 ≈ P(𝒵) = 𝒵 . n 𝒳

(1)

Przybliżenie staje się coraz lepsze w miarę zwiększania $n.$ Ścisłe sformułowanie tego stwierdzenia nazywa się Prawem Wielkich Liczb (PWL) i jest podstawą rachunku prawdopodobieństwa. Nie trzeba jednak wielkiej wiedzy z tej dziedziny, wystarczy trochę zdrowego rozsądku, aby w przybliżoną równość (1) uwierzyć i ją zrozumieć. Z naszego punktu widzenia ważne jest wykorzystanie PWL do oszacowania nieznanej liczby SAW-ów. Wiemy, że $| 𝒳 = 4k.$

$n1 k |n 4 ≈ 𝒵$ jest oszacowaniem liczby $Zk$ prostą metodą Monte Carlo.

W zasadzie możemy osiągnąć dowolną dokładność dla odpowiednio dużych $|n.$ Zauważmy, że bardzo łatwo wylosować drogę ze zbioru $𝒳$ w taki sposób, aby dla każdego $x ∈𝒳$ prawdopodobieństwo wybrania $x$ było jednakowe, równe $P(x) = 1/4k.$ Na każdym kroku wybieramy jeden z 4 kierunków losowo. Możemy dwa razy rzucić monetą i umówić się: jeśli otrzymamy OO, to $, N$ jeśli OR, to $W,$ jeśli RO, to $|S,$ jeśli RR, to $E.$ W ten sposób generujemy ciąg kodujący drogę "błądzenia przypadkowego". Komputer pozwala szybko "rzucić monetą" np. $|(2k)⋅10000$ razy i wygenerować $n = 10000$ takich dróg.

Niestety, prosta metoda Monte Carlo jest bardzo nieefektywna, bo dla dużych $k$ prawdopodobieństwo $P (𝒵) = 𝒵 / 𝒳$ szybko zbliża się do zera. Oczekiwanie na przypadkowe trafienie w zbiór $𝒵,$ czyli wylosowanie SAW-a, zaczyna przypominać poszukiwanie igły w stogu siana. Możliwym rozwiązaniem tego problemu jest metoda "wzrostu" zaproponowana przez Ariannę i Marshalla Rosenbluthów polega na losowaniu kolejnych kroków błądzenia przypadkowego tylko spośród "dopuszczalnych punktów", to znaczy punktów wcześniej nieodwiedzonych. W każdym kroku, z wyjątkiem pierwszego, mamy co najwyżej $3$ możliwości. Rysunek 2 (następna strona) pokazuje kolejne kroki prowadzące do zbudowania SAW-a z rysunku 1. Widać, że kolejne kroki wybieraliśmy spośród

4,3,3,3, 2,3,2,2, 3,2,3,3, 2,1,3

możliwych. Nasz SAW został zatem wylosowany z prawdopodobieństwem

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 --⋅-⋅ -⋅--⋅ -⋅--⋅--⋅--⋅ --⋅--⋅--⋅--⋅ --⋅- ⋅-. 4 3 3 3 2 3 2 2 3 2 3 3 2 1 3

Powiedzmy ogólniej, że przy budowaniu SAW-a $x$ o długości $l$ mamy kolejno

m1

możliwości. Nie jest przy tym wykluczone, że w pewnym kroku $l < k$ nie mamy żadnej dalszej możliwości, $ml+1$ Mówimy wtedy, że powstał nieprzedłużalny SAW o długości $l.$ Niech $𝒴$ oznacza zbiór wszystkich SAW-ów o długości $l = k$ oraz nieprzedłużalnych SAW-ów o długości $l < k.$ Oczywiście, $|𝒵 ⊂𝒴.$

Prawdopodobieństwo wylosowania SAW-a $|x ∈𝒴$ jest równe

-1- -1- 1-- P (x) = m ⋅ m ⋯ m .

Niech, dla $x ∈𝒴,$ będzie

⎧ W(x) = ⎪⎪⎨m1m⋅2 jeżeli x∈ 𝒵(a więc l = k); ⎪⎪0 jeżeli x /∈ 𝒵(a więc l < k i ml+1 . ⎩

Można interpretować $W(x)$ jako "wagę" SAW-a $x.$ Powtórzmy teraz losowanie metodą wzrostu $| n$ razy. Powstaje $| n$ losowych SAW-ów $1,X2,...,XnX$ ze zbioru $𝒴.$ Obliczamy "średnią wagę" wylosowanych SAW-ów i zauważamy, że

n )≈P(x)W(x)=1= 𝒵 . ¯W = 1- W(X iQx>𝒴Qx>𝒵 nQi 1

(2)

W rezultacie,

$W¯ ≈ 𝒵$ jest to oszacowaniem liczby $|Z k$ "metodą wzrostu".

Dodajmy komentarz na temat wzoru (2). Po lewej stronie mamy średnią wagę obliczoną dla wylosowanych $n$ elementów, czyli "średnią z próbki". Po prawej mamy też średnią, ale "w całej populacji $| 𝒴$ ". Wagi $| W$ zostały specjalnie dobrane tak, aby ta średnia populacyjna była równa $𝒵 .$ Jeśli $|n$ jest dużo, dużo mniejsze od $𝒵 ,$ to obliczenie $W¯$ jest możliwe, a bezpośrednie policzenie elementów $|𝒵$ jest niemożliwe. Kluczowy jest fakt, że średnia próbkowa przybliża średnią populacyjną. To jest odrobinę ogólniejsza wersja PWL. Poniższy rachunek jest próbą przekonania Czytelnika, że aproksymacja we wzorze (2) jest równie intuicyjna jak ta we wzorze (1). Niech $|ns$ oznacza liczbę tych spośród wylosowanych elementów $| i,X$ dla których $| i)=s W(X$ i $| 𝒲s = {x ∈𝒴 W(x) = s}.$ Podobnie jak we wzorze (1), na mocy PWL, $ns/n ≈P (𝒲s).$ Mamy zatem

ns 1- n i)=Qsn≈QsP(𝒲s)=QsQP(x)=QP(x)W(x). n Q W(X sssx>𝒲sx i 1

Poniższa tabelka podaje dokładne wartości $| Zk,$ oszacowania tych liczb prostą metodą Monte Carlo (MC) i "metodą wzrostu". Dla obu metod, dla każdego $k,$ liczba symulowanych doświadczeń losowych była równa $|n = 10000.$

$|---|--------------------------|----------|----------------| | k | Zk |proste MC |metoda wzrostu | |---|--------------------------|----------|----------------| |-1-|------------4-------------|----4-----|-------4--------| |-2-|------------12-------------|--11,99----|-----12,00------| | 3 | 36 | 36,27 | 36,00 | |-4-|-----------100------------|--98,74---|-----100,18------| |---|--------------------------|----------|----------------| |-5-|-----------284------------|--282,73---|-----283,21-----| |10-|----------44100-----------|-41628,47--|----44041,61----| |20 | 897697164 |659706977 | 893264097 | |---|--------------------------|----------|------------13--| |30-|-----16741957935348-------|----0-----|--1,66208-⋅10----| |40-|---300798249248474268-----|----0-----|-3,067392-⋅1017--| -50--5292794668724837206644---------0--------5,261991-⋅1021--|$

Skąd wzięły się niepokojące zera w dolnej części tabelki? To się natychmiast stanie jasne, jeśli Czytelnik zechce obliczyć

k |Zk/4

dla

|k = 30,40,50.

Z tabelki widać, że metoda wzrostu radzi sobie nie najgorzej nawet dla sporych wielkości

k.

Trzeba jednak dodać, że ta metoda, wynaleziona w latach pięćdziesiątych ubiegłego wieku, może być znacznie ulepszona. Istnieją bardziej efektywne algorytmy Monte Carlo, przeznaczone do zliczania SAW-ów. Niektóre z nich wykorzystują zupełnie inne pomysły niż metoda wzrostu. Niestety, opowieść o tym przekracza ramy tego artykułu.

SAW-y pojawiają się w fizyce jako najprostszy model budowy polimerów. Są to duże cząsteczki mające postać łańcucha złożonego z monomerów. Znajdują się w najróżniejszych materiałach: od włókien tkanin, tworzyw sztucznych, gumy i celulozy, aż po białka w żywych organizmach. Struktura przestrzenna łańcucha wpływa na własności fizyczne polimeru. W ogromnym uproszczeniu tę strukturę reprezentuje trójwymiarowy SAW (o długości $|k$ od ok. $3 |10$ do $5 10$ ). W przestrzeni trójwymiarowej SAW definiuje się bardzo podobnie jak dwuwymiarowy: jako ciąg punktów o trzech współrzędnych całkowitych, w którym nie ma powtórzeń. Okazuje się, że pewne istotne cechy polimerów są związane np. z odległością (euklidesową) między końcami łańcucha. Badanie tej wielkości było zasadniczą motywacją przytoczonych prac. Zadanie zliczenia SAW-ów pojawiło się "przy okazji" i zafrapowało matematyków. Na zakończenie przytoczę przypuszczenie, którego dotąd nie udało się ani udowodnić, ani obalić. Powróćmy do dwuwymiarowych SAW-ów o długości $k$ i liczb $Zk.$ Przypuszcza się, że istnieje granica

lim --Zk---= A, k ∞ k11~32µk

gdzie $|0 < A$ i $2,6256 < µ ⩽ 2,6792.$ Coś, co zaczyna się jak przedszkolna zabawa, prowadzi do takich tajemnic.