Natodia blog

Każdy, kto spotkał się z całkami na pewno zna kilka sposobów ich liczenia. W najprostszych przypadkach możemy wyznaczyć całkę nieoznaczoną z danej funkcji analitycznie. Co z trudniejszymi przypadkami? Z nimi możemy sobie poradzić o ile chcemy znać wartość całki oznaczonej z pewną skończoną dokładnością.

Jest kilka niezwykle popularnych metod – metoda prostokątów (dzielimy całkowany przedział na n mniejszych podprzedziałów i przybliżamy ich pole prostokątem), metoda trapezów (jak w metodzie poprzedniej tylko tym razem przybliżamy pole trapezem), metoda parabol (to już chyba jasne?). Mają one kilka wspólnych cech: są deterministyczne, teoretycznie możemy przy ich pomocy uzyskać wartość całki z dowolną dokładnością (ogranicza nas „tylko” czas i szybkość procesora). Gdy czas jest dla nas kluczowy (i nie mamy go za dużo) a nasza funkcja jest wielowymiarowa, możemy spróbować szczęścia używając niedeterministycznej metody całkowania numerycznego – metody Monte Carlo.

Aby ją zobrazować posłużę się przykładem zaczerpniętym z Wikipedii. Załóżmy, że chcemy wyznaczyć przybliżoną wartość liczby pi. Możemy wymyślić sobie następujący eksperyment: weźmy kwadrat o boku o znanej długości i wpiszmy w niego koło.

Teraz możemy losować punkty w środku kwadratu i zliczać te, które „trafiły” w koło. Sprawdzamy to licząc odległość punktu od środka kwadratu i sprawdzając, czy jest ona mniejsza (bądź równa) promieniowi naszego koła. Oznaczmy liczbę wylosowanych punktów jako n i liczbę punktów znajdujących się w kole jako t. Możemy policzyć przybliżony stosunek pola koła do pola kwadratu: k = t / n. Wiedząc, że k = pi * r² / 4r² wyznaczamy pi. Możemy zwiększać dokładność naszego szacowania zwiększając liczbę losowanych punktów.

Mamy już właściwie wszystkie niezbędne informacje aby liczyć wartości całek oznaczonych metodą Monte Carlo. Załóżmy, że chcemy policzyć całkę oznaczoną na przedziale <a, b> z funkcji f(x). Na początku losujemy n wartości z przedziału całkowania (x₁, x₂, …, x_i, …, x_n) i liczymy średnią s z wartości f(x_i). Wartość całki wynosi c = s * (b – a).

Skąd bierze się powyższy wzór? Intuicja jest następująca: licząc naszą średnią s wyznaczamy średnią „wysokość” funkcji. Znając ją liczymy pole prostokąta o bokach (b – a) oraz s.

Implementacja metody nie jest trudna. Przygotowałem przykładową implementację w języku C. Zachęcam do zapoznania się z kodem. Porównuje on metody prostokątów, trapezów oraz Monte Carlo. Jak można zauważyć, w przypadku jednowymiarowym metody deterministyczne są szybsze i dokładniejsze.

Niniejszy wpis jest tylko krótkim wstępem do metody Monte Carlo, która jest wykorzystywana do granic możliwości przez fizyków (przykładem jest Wrocławska Grupa Neutrinowa, która używając tej metody symuluje oddziaływanie neutrin z materią). W następnym wpisie pokażę przykładowe zastosowanie tej metody do konkretnego problemu fizycznego.