Nicholas Metropolis Amerykanin pochodzenia greckiego, podobnie jak bohaterowie poprzedniej opowieści pracował w Las Alamos przy bombie atomowej i też grzebał przy metodzie Monte Carlo. Być może miał jakieś złe doświadczenia z jej zastosowaniem do termodynamiki, bo opracował swój algorytm sprawdzający się w przypadkach takich jak ten z poprzedniej notki, o niebo lepiej od zwykłej Monte Carlo.
Opis algorytmu
Przypomnijmy sobie sytuację z poprzedniego odcinka: chodzi o to, żeby losować stany (konfiguracje) głównie w okolicach silnego choć niezwykle wąskiego maksimum funkcji exp{–Estanu/kT}. A jak się to robi?
-
Losujemy jakiś stan początkowy i obliczamy jego energię.
-
Losujemy jakiś inny stan, gdzieś w pobliżu tego początkowego.
-
Jeśli ten nowy stan ma mniejszą energię, to zapominamy o starym – nowy zajmuje miejsce starego. Jeśli natomiast nowy ma większą energię, to czy zajmie on miejsce starego określa prawdopodobieństwo:
p(Enowy, Estary) = exp{–(Enowy–Estary)/kT}
W szczegółach wygląda to tak: losujemy liczbę z zakresu (0,1):
los = double(rand())/RAND_MAX
Jeśli los<p(Enowy, Estary) dokonujemy zamiany, w przeciwnym razie pozostawiamy stary stan.
- Wracamy do punktu 2, to znaczy losujemy kolejny stan.
W ten sposób obliczamy sobie swoistą, skokową ewolucję stanu. Jeśli model nie jest jakiś pokręcony, kolejne losowania przybliżą nasz stan do miejsca gdzie obliczone wagi Boltzmana będą duże. Nawiązując do schematycznego rysunku z poprzedniej notki znajdowane kolejne rozkłady powinny wyglądać następująco:
Jak wygląda „ruch” losowanego stanu? Zasadniczo, ze względu na pkt. 3 układ zmierza w stronę najniższej energii. Możliwe są jednak „wzbudzenia” do stanów o wyższych energiach. Prawdopodobieństwo wzbudzeń „reguluje” temperatura T. Im jest większa, tym większa jest szansa, że układ na skutek losowania „przeskoczy” do stanu o wyższej energii. Te dwa procesu powodują, że w pewnym momencie ustala się swoista równowaga.
I teraz zupełnie bez dowodu podam twierdzenie w formie wyjątkowo nieścisłej i może nawet nie do końca prawdziwej: Przy dostatecznie długim czasie „błądzenia”, średnie zachowanie układu odtwarza zespół kanoniczny (pod warunkiem, że model ma dobre własności). Innymi słowy gdy uruchomimy sobie proces i przy obliczaniu odpowiednich średnich będziemy brali wagi Boltzmana z iluśtam ostatnich kroków procesu, dostaniemy poszukiwaną reprezentatywną próbkę.
Przykład niskiej temperatury
Nie będę od nowa opisywał modelu, bo zrobiłem to w poprzednich odcinkach. Tu napiszę, jak zrealizowałem metropolisowskie błądzenie. Wystartowałem od losowego stanu (konfiguracji). W każdym kroku losowane były jakieś dwa poziomy „jednocząstkowe” – jeden zajęty, drugi wolny i następna konfiguracja miała zajętość zamienioną miejscami na tych poziomach. Potem według procedury opisanej powyżej podejmowana była decyzja czy nowy stan zajmie miejsce starego.
Dla małych temperatur wszystkie cząstki powinny zapełnić stany o najniżej energii. Przy okazji sprawdzę sobie ile mniej więcej kroków potrzeba, żeby na skutek losowań układ znalazł się w stanie „równowagi”. Model identyczny z poprzednim, parametry: n=80, k=40, T=0,001; Ei = 0,2∙i (poziomy zajęte to gwiazdki, I to numer iteracji):
* * * * * * * ** ***** * * *** * ** * ** *** ** * **** * ** ** I=0
******* ************* * **** * * ** ** * * * ** * * * * I=100
********* ****************** * ***** * * * * * * * I=200
**************************** **** ** * * *** * I=300
********************************* ** *** * * I=400
********************************* ** *** * * I=500
********************************* **** * ** I=600
************************************* * ** I=700
************************************* * ** I=800
************************************* * ** I=900
************************************* *** I=1000
**************************************** I=1100
**************************************** I=1200
**************************************** I=1300
**************************************** I=1400
Jak widać, dla tak niewielkiej temperatury stan „równowagi” ustala się dość szybko. Niemniej gdy układ zostanie już dobrze „schłodzony” to znaczy zbliży się do stanu najniższej energii, potrafi się „zamrozić” – widać to na przykład pomiędzy 600-tną a 900-tną iteracją. Takie „zastopowanie” procesu wynika z tego, że dla małych T, decydująca o skoku wielkość exp{–(Enowy–Estary)/kT} jest niewielka i rzadko da się wylosować wystarczająco małą wartość zmiennej los. Jeśli chodzi o ruch „w drugą stronę”, to w stany o mniejszej energii znów nie tak łatwo trafić, bo tych o większej energii jest o wiele, wiele więcej.
Czasami, czego akurat w powyżej namalowanej historii nie widać, wylosowane prawdopodobieństwo pozwoli czasami na „przeskok” którejś „cząstki” na jakiś wyższy poziom, bo temperatura choć mała, to jednak jest większa od zera.
Przybliżony rozkład Fermiego-Diraca
No i pora na właściwą próbę. Dobrałem następujące parametry: n=200, k=80, T=2, Ei = 0,5∙i. Wykonałem 300 tysięcy iteracji, przy obliczaniu średniego obsadzenia stanu „jednocząstkowego” wziąłem pod uwagę wagi z ostatnich 200 tysięcy iteracji.
Jak widać wykres całkiem, całkiem. Nie jest może jakoś super gładki, ale dość udolnie przybliża się do wyniku ścisłego. W obszarze największej zmienności poszczególne wartości różnią się o ok. 1-2%. W kilku punktach niedokładność sięga kilku procent. Rakiety w kosmos to może bym nie wysyłał w oparciu o te dane, ale gdyby chodziło o szybkie szacowanie to powinno wystarczyć. Szczególnie, że gdyby ktoś chciał policzyć wagi wszystkich możliwych stanów czyli konfiguracji, będzie tego 1,64728x1057 kombinacji.
* * *
Seria notek o obliczaniu wag Boltzmana, wbrew pozorom, nie miała na celu zreferowania zasad obowiązujących przy prowadzeniu symulacji komputerowych. Nie jestem bowiem specjalistą i trudno mi powiedzieć gdzie na potencjalnego amatora programowania czekają zdradzieckie rafy. No bo jak wiadomo diabeł tkwi w szczegółach i chodzi o to, żeby takie szczegóły móc rozpoznawać i pokonywać lub omijać. Przy pisaniu zależało mi raczej na tym, żeby pokazać, że mimo wszystko nie jest to dziedzina przeznaczona dla wybrańców, ale każdy kto ma trochę obycia komputerowego może spróbować swoich sił.
Każdą serię tematyczną staram się zakończyć jakimś lżejszym materiałem. I taka też ma być kolejna notka na którą już dziś zapraszam
