Sprawdzanie rozkładu k fermionów na n poziomach jest, jak okazało się w poprzednim odcinku dość trudne dla całkiem niewielkich liczb n i k. Od czego jednak metody przybliżone? Może zamiast przeliczać wszystkie możliwe konfiguracje, wystarczy wziąć „reprezentatywną” próbkę?
Metoda Monte-Carlo
Nazwa pochodzi od losowania, bo Monte Carlo to nazwa kasyna, gdzie grając w ruletkę losuje się numerki[1]. W opisywanej metodzie też się losuje numerki. Jej szkolnym przykładem zastosowania jest „policzenie” liczby π. Oczywiście nie chodzi o to, żeby się dowiedzieć ile ta liczba wynosi, bo to wiemy skądinąd, ale żeby przetestować metodę.
Żeby obliczyć przybliżone π, program losuje sobie dwie liczby – współrzędne punktu na płaszczyźnie. W C/C++ istnieje nawet funkcja losująca rand(), która zwraca liczby całkowite z zakresu od 0 do 32767 czyli RAND_MAX – podobno istnieją środowiska gdzie stała RAND_MAX jest inna, ale z tego co się zorientowałem, to powszechnie stosowana wartość. Im więcej punktów wylosujemy, tym gęściej pokryjemy kwadrat jak na rysunku poniżej. Gdybyśmy policzyli, ile punktów jest wewnątrz pola zakreślonego przez ćwiartkę okręgu, to powinno ich być mniej więcej π/4 wszystkich losowań.
Dla porządku przedstawiam kod w C++ „liczący” liczbę π[2], dla tych, co chcieliby spróbować:
srand( time( NULL ) );
const int IlProb = 100000;
int rKwadrat = RAND_MAX*RAND_MAX;
int suma=0, x, y;
for( int i=0; i<=IlProb; i++ )
{
x = rand();
y = rand();
if( x*x+y*y < rKwadrat ) suma++;
}
double pi = (4.0*suma)/IlProb;
Strzelamy w konfiguracje
Nie będę się skupiał nad szczegółami rozwiązania, bo i trudno je nazwać rozwiązaniem: wyniki okazały wysoce niezadowalające. Model jest identyczny z tym rozważanym w poprzedniej notce – układ złożony z k cząstek, n dostępnych poziomów, energia rosnąca liniowo z numerem poziomu, liczę średnią zajętość każdego poziomu – tyle, że zamiast brania pod uwagę wszystkich stanów układu złożonego, losuję kolejne konfiguracje.
Reszta programu działa bez zmian: Dla każdego losowania obliczana była całkowita energia konfiguracji Ec (trzeba wysumować zajęte poziomy), a z energii, waga według wzoru:
pT ~ exp{–Ec/kT}
Wyliczony rozkład zajętości poziomów powinien być podobny do rozkładu Fermiego-Diraca, ale nie był. Oto uzyskany rozkład dla k=20 cząstek na n=50 poziomach – wykres niebieski to uśrednienie 2 milionów losowań, zielony odpowiada już 10 milionom losowań. Linia czerwona pokazuje co powinno wyjść na podstawie wzoru książkowego.
Wykresy wyszły niefajne[3]. Jak widać wylosowana próbka jakoś nie chce być reprezentatywna. Dlaczego?
Problem w tym, że najczęściej (prawie zawsze) trafiałem na tragicznie malutkie wagi. Te konfiguracje, „odległe” od stanu o najniższej energii, miały dość dużą energię i eksponenta, która brała ją z minusem, dawała praktycznie zero. Od czasu do czasu „wstrzeliwałem się” w miarę blisko, ale były to rzadkie przypadki i z powodu małej ilości nie dawały się ładnego uśredniania.
Główny wkład do sumy statystycznej dają konfiguracje postaci:
- poziomy o małej energii – wszystkie zajęte,
- poziomy o średniej energii – rozpiętość energii rzędu kilku kT – pozajmowane mniej więcej pół na pół,
- poziomy o wyższej energii – puste.
Trafić na taki stan było bardzo ciężko. Schematycznie sytuację można przedstawić symbolicznie jak na wykresie poniżej:
Dziedziną funkcji są możliwe konfiguracje. Wybieramy kolejne na chybił-trafił (to te czerwone kropki), ale istotny wkład może dać jedynie bardzo wąskie i bardzo wysokie maksimum, które niezwykle trudno wylosować.
Co dalej?
A może jest jakaś „ulepszona metoda Monte Carlo”, która wybiera tylko takie rozkłady dla których wagi są znaczące, a opuszcza te nieważne? Jest, i nazywa się metodą Metropolisa, ale o tym w następnym odcinku.
[1] MC powstało przy okazji tworzenia amerykańskiej bomby atomowej, kiedy trzeba było symulować zachowanie neutronów. Pomysłodawcą metody był Stanisław Ulam – matematyk rodem z Polski. Pomysł wdrożyli w życie Ulam i von Neumann, który wymyślił kryptonim „Monte Carlo” – wujek Ulama grał w kasynie „Monte Carlo” w Monako. Kryptonim był potrzebny, bo w Los Alamos wszystko było tajne.
[2] Uzyskiwana dokładność dla 100 tysięcy prób jest rzędu promila.
[3] Zwykle w takich przypadkach rodzą się w człowieku wątpliwości: Może generator liczb pseudolosowych jest za mało losowy? A może sumowanie miliona liczb daje nieprawidłowe wyniki ze względu na błędy zaokrągleń? Żeby uspokoić sumienie wynalazłem w sieci inny generator oraz posortowałem uzyskane losowania ze względu na wagi (kiedy sumuje się liczby dodatnie „od najmniejszej do największej” błędy zaokrągleń mniej bolą). Niestety nie zmieniło to niedobroci uzyskanego rozwiązania.
