Metody Numeryczne - Metoda Monte Carlo

Metoda Monte Carlo (MC) to jedna z najbardziej fascynujących metod numerycznych w historii informatyki i fizyki. W przeciwieństwie do klasycznych, sztywnych kwadratur (jak metoda trapezów czy Simpsona), Monte Carlo opiera się na losowości i statystyce.

Jej narodziny są bezpośrednio powiązane z Projektem Manhattan (budową amerykańskiej bomby atomowej) oraz powstaniem pierwszych komputerów na świecie.

Oto historia tego algorytmu:

Pierwsze ślady: Igła Buffona (1777 r.)
Zanim pojawiły się komputery, francuski filozof i matematyk Georges-Louis Leclerc, hrabia Buffon, przeprowadził słynny eksperyment probabilistyczny. Rzucał igłą na deskowaną podłogę i liczył, jak często igła przecina linie styku desek. Udowodnił, że na podstawie liczby rzutów i trafień można wyznaczyć przybliżoną wartość liczby \(\pi \). Był to pierwszy w historii udokumentowany przypadek użycia losowości do rozwiązania problemu matematycznego.
Narodziny w Los Alamos: Pasjans Stanisłwa Ulama (1946 r.)
Prawdziwym twórcą nowoczesnej metody Monte Carlo jest Stanisław Ulam – genialny polski matematyk należący do lwowskiej szkoły matematycznej, który po emigracji do USA pracował w ośrodku tajnych badań wojskowych w Los Alamos.
Anegdota o pasjansie: W 1946 roku, chorując w szpitalu, Ulam grał w pasjansa. Zastanawiał się, jakie jest matematyczne prawdopodobieństwo, że układ kart pozwoli na wygranie całej partii. Kombinatoryczne obliczenie tego było potwornie trudne. Ulam wpadł na rewolucyjny pomysł: zamiast liczyć skomplikowane wzory, lepiej... po prostu rozegrać pasjansa 100 razy, zapisać liczbę wygranych i wyciągnąć średnią statystyczną. Natychmiast przełożył tę intuicję na problem fizyczny, nad którym głowił się w Los Alamos – ruchu i dyfuzji neutronów w uranie podczas reakcji łańcuchowej.
Współpraca z von Neumannem i kryptonim „Monte Carlo”
Ulam podzielił się swoim pomysłem z Johnem von Neumannem, jednym z ojców współczesnej informatyki. Von Neumann od razu dostrzegł potęgę tego podejścia.
Problem obliczeniowy: Ludzki mózg nie był w stanie wykonać milionów losowań papierem i ołówkiem. Von Neumann zaadaptował pomysł Ulama do uruchomienia na ENIAC-u – pierwszym w pełni elektronicznym komputerze na świecie.
Kryptonim: Ponieważ badania w Los Alamos były ściśle tajne, projekt musiał otrzymać kryptonim. Nicholas Metropolis, inny fizyk z zespołu, zaproponował nazwę „Monte Carlo”. Było to nawiązanie do słynnego kasyna w Monako, w którym wujek Stanisława Ulama regularnie przegrywał pieniądze w ruletkę. Nazwa idealnie oddawała hazardowy, losowy charakter metody.
Oficjalna publikacja (1949 r.)
W 1949 roku Stanisław Ulam i Nicholas Metropolis opublikowali przełomowy artykuł pt. The Monte Carlo Method w czasopiśmie Journal of the American Statistical Association. Od tego momentu metoda weszła do oficjalnego obiegu naukowego, rewolucjonizując fizykę statystyczną, chemię kwantową, a z czasem także ekonomię i grafikę komputerową (tzw. Ray Tracing).

Aby zrozumieć zasadę metody Monte Carlo wyobraź sobie, iż chcesz wyznaczyć pole koła wpisanego w kwadrat o boku równym 2 (pole to co do wartości jest równe liczbie pi, ale na razie udajmy, że o tym nie wiemy). W tym celu wyznaczasz wewnątrz kwadratu dużo losowych punktów. Następnie liczysz poszczególne punkty i określasz, ile z nich jest wewnątrz koła. Pole koła jest w przybliżeniu równe:

\[P_{\circ} \approx \dfrac{n_{\circ}}{n} \cdot P_{\square}\]
\(P_{\circ}\)	: pole koła
\(P_{\square}\)	: pole kwadratu
\(n_{\circ}\)	: liczba punktów w kole
\(n\)	: liczba wszystkich punktów

Przykład:

Oto odpowiedni skrypt JavaScript symulujący obliczanie powierzchni koła dla podanego przykładu:

JavaScript

<html>
<head>
  <title>Wyznaczanie liczby PI
  metodą Monte Carlo</title>
</head>
<body>

<div style="overflow-x: auto;"
     align="center">
  <table
  border="0"
  cellpadding="4"
  style="border-collapse:
         collapse">
    <tr>
      <td nowrap>
        <form
        name="frm1"
        style="text-align:
               center;
               background-color:
               #E7E7DA">
          <b>
          Obliczanie pola
          koła wpisanego
          <br/>
          w kwadrat o boku 2
          <br/>&nbsp;&nbsp;
          za pomocą metody
          Monte Carlo
          &nbsp;&nbsp;</b><br/>
          (C)2026 mgr
          Jerzy Wałaszek<hr>
          Podaj liczbę
          punktów<br/>
          do wygenerowania:<br/>
          <input
          type="text"
          name="inp_n"
          size="16"
          value="10000"
          style="text-align:
                 right">
          <hr>
          <input
          onclick="js_p();"
          value="Oblicz pole koła"
          name="B1"
          type="button">
          <hr>
          Pole koła wynosi:
          <div id="out1">.</div>
        </form>
      </td>
    </tr>
  </table>
</div>

<script language="javascript">

// Przykładowa aplikacja
// obliczająca pole koła
// wpisanego w kwadrat
// za pomocą metody
// Monte Carlo
//----------------------
// (C)2004 mgr Jerzy Wałaszek

function js_p()
{
  var n = parseInt(document
          .frm1.inp_n.value);
  var nk,s,x,y,i;

  if(isNaN(n) || (n < 2))
    s = "<font color=red><b>" +
        "Popraw dane" +
        "</b></font>";
  else
  {
    nk = 0
    for(i = 0; i < n; i++)
    {
      x = Math.random() * 2;
      y = Math.random() * 2;
      if(Math.sqrt((x-1) *
        (x-1) + (y-1) *
        (y-1)) <= 1) nk++;
    }
    s = 4 * nk / n;
    s = Math.round(s * 100000) /
        100000;
  }
  document.getElementById("out1")
  .innerHTML = s
}

</script>

</body>
</html>

Oczywiście wynik jest bliski liczbie \({\pi = 3{,}1415926535\dots}\) dopiero dla dużych wartości \(n\). Zaczynają wtedy obowiązywać prawa dużych liczb i pomimo przypadkowości wyboru punktów, pojawia się ogólna prawidłowość rozkładu równomiernego. Ponieważ punkty rozkładają się równomiernie w obrębie pola kwadratu, to stosunek liczby punktów wewnątrz koła do liczby wszystkich punktów w kwadracie przybliża się do stosunku pola koła do pola kwadratu. Stąd właśnie pochodzi nasz wzór:

\[P_{\circ} \approx \frac{n_{\circ}}{n} \cdot P_{\square}\]

Wzór ten jest podstawą wyznaczania wartości całki oznaczonej za pomocą metody Monte Carlo, czyli losowania punktów. Zasada jest następująca:

Dla danej funkcji \(f(x)\), której całkę oznaczoną chcemy obliczyć w przedziale całkowania \({[x_p ; x_k]}\), wyznaczamy prostokąt obejmujący pole pod wykresem tej funkcji o wysokości \(h\) i długości podstawy \({(x_k - x_p)}\). W dalszej kolejności losujemy \(n\) punktów i zliczamy te punkty \(n_w\), które wpadają w pole pod wykresem funkcji. Wartość całki wyraża się wzorem przybliżonym:

\[
\int\limits_{x_p}^{x_k} f(x) \, dx \approx \frac{n_w}{n} \cdot h \cdot (x_k - x_p)
\]

Otrzymany wzór ma kilka wad. Na przykład w ogólnym przypadku trudno wyznaczyć wysokość \(h\). Również kłopoty pojawiają się, gdy funkcja zmienia znak w przedziale całkowania. Dlatego częściej jako metodę Monte Carlo przyjmuje się metodę, która wyznacza średnią z wartości funkcji w przedziale całkowania na podstawie serii losowo wybranych współrzędnych \(x\). Następnie średnia ta jest mnożona przez długość przedziału całkowania i otrzymujemy przybliżoną wartość całki oznaczonej. Wzór ma następującą postać:

\[
\int\limits_{x_p}^{x_k} f(x) \, dx \approx (x_k - x_p) \cdot \frac{\sum \limits_{i=1}^{n} f(x_{los})}{n}
\]

do podrozdziału do strony

Specyfikacja problemu

Dane wejściowe

\(x_p\)	: początek przedziału całkowania; \({x_p \in \mathbb{R}}\).
\(x_k\)	: koniec przedziału całkowania; \({x_k \in \mathbb{R}}\).
\(n\)	: liczba losowanych punktów; \({n \in \mathbb{N}}\).
\(f(x)\)	: funkcja rzeczywista, której całkę liczymy

Dane wyjściowe

\(s\)

: przybliżona wartość całki oznaczonej funkcji \(f(x)\) w przedziale \({[x_p ; x_k]}\); \({s \in \mathbb{R}}\).

Zmienne pomocnicze

\(h\)	: szerokość przedziału całkowania; \({h \in \mathbb{R}}\).
\(i\)	: licznik punktów; \({i \in \mathbb{N}}\).
\(x_{los}\)	: punkt wybierany losowo z przedziału całkowania; \({x_{los} \in \mathbb{R}}\).

Lista kroków

K01:	\(s \leftarrow 0; \; h \leftarrow x_k - x_p\)
K02:	Dla \(i = 1,2,\dots,n\) wykonuj kroki K03...K04
K03:	\(x_{los} \leftarrow \text{liczba losowa} \in [x_p;x_k]\)
K04:	\(s \leftarrow s + f(x_{los})\)
K05:	\(s \leftarrow \dfrac{s}{n} \cdot h\)
K06:	Pisz \(s\)
K07:	Zakończ

Obliczenia rozpoczynamy od pobrania informacji o przedziale całkowania oraz o ilości punktów losowych, które należy wygenerować w celu obliczenia wartości średniej funkcji w tym przedziale. Dokładność metody rośnie wraz ze wzrostem \(n\).

W zmiennej \(s\) będziemy obliczać sumy wartości funkcji. Zmienna ta posłuży później do wyliczenia średniej oraz samej całki oznaczonej. Na początku obliczeń ustawiamy ją na 0. W zmiennej \(h\) zapamiętujemy szerokość przedziału całkowania. Wartość ta jest wykorzystywana zwykle przy generacji liczby losowej oraz na końcu przy obliczeniu wartości całki.

Rozpoczynamy pętlę iteracyjną kontrolowaną przez zmienną \(i\). Pętla ta wykona się \(n\) razy. Wewnątrz pętli generujemy liczbę pseudolosową \(x_{los}\) leżącą w przedziale \({[x_p ; x_k]}\). Metoda generacji zależy od wybranego języka programowania, który udostępnia odpowiednie funkcje pseudolosowe. Po wyznaczeniu liczby pseudolosowej \(x_{los}\) obliczamy wartość funkcji w tym punkcie i dodajemy ją do sumy \(s\). Gdy pętla się zakończy, wyliczamy średnią wartość funkcji w przedziale całkowania, mnożymy ją przez długość tego przedziału i otrzymujemy przybliżoną wartość całki oznaczonej. Wypisujemy wyniki i kończymy algorytm.

do podrozdziału do strony

W przedziale \({[0 ; 1]}\) całka funkcji \({f(x) = x^2 + 2x}\) ma dokładną wartość \({1{,}333\dots}\) W naszym programie liczymy \({n = 10000}\) losowych punktów.

C++

// Obliczanie całki oznaczonej
// metodą Monte Carlo
// Metody numeryczne 0084
// ---------------------------
// (C)2026 mgr Jerzy Wałaszek

#include <iomanip>
#include <iostream>
#include <cstdlib>
#include <time.h>

using namespace std;

// Tutaj definiujemy funkcję
double f(double x)
{
  return(x * x + 2 * x);
}

// Program główny

int main()
{
  // liczba punktów losowych
  const int N = 10000;
  double xp,xk,s,dx;
  int i;

  // 3 cyfry po przecinku
  cout << setprecision(3)
  // format stałoprzecinkowy
       << fixed;

  cout << "Obliczanie calki oznaczonej\n"
          "    Metoda Monte Carlo\n"
          "---------------------------\n"
          "(C)2026 mgr J.Walaszek  I LO\n\n"
          "f(x) = x * x + 2 * x\n\n"
          "Poczatek przedzialu calkowania\n\n"
          "xp = ";
  cin >> xp;
  cout << "\nKoniec przedzialu calkowania\n\n"
          "xk = ";
  cin >> xk;
  cout << endl;
  srand(time(NULL));
  s  = 0;
  dx = xk - xp;
  for(i = 1; i <= N; i++)
    s += f(xp+((double)rand()/(double)(RAND_MAX+1)*dx));
  s = dx * s / N;
  cout << "Wartosc calki wynosi : " << s
       << endl << endl;
  system("pause");
  return 0;
}

Wynik:

Obliczanie  całki oznaczonej
za pomocą metody Monte Carlo
----------------------------
(C)2026 mgr J.Wałaszek  I LO

f(x) = x * x + 2 * x

Początek przedziału całkowania

xp = 0

Koniec przedziału całkowania

xk = 1

Wartość całki wynosi : 1.331

Naciśnij Enter...

Python (dodatek)

# Obliczanie całki oznaczonej
# metodą Monte Carlo
# Metody numeryczne 0084
# ---------------------------
# (C)2026 mgr Jerzy Wałaszek

from random import uniform

def f(x):
    return x * x + 2 * x

# Program główny

# Liczba punktów losowych
n = 10000

print("Obliczanie  całki oznaczonej")
print("za pomocą metody Monte Carlo")
print("----------------------------")
print("(C)2026 mgr J.Wałaszek  I LO")
print()
print("f(x) = x * x + 2 * x")
print()
print("Początek przedziału całkowania")
print()
xp = float(input("xp = "))
print()
print("Koniec przedziału całkowania")
print()
xk = float(input("xk = "))
print()

s = 0
dx = xk - xp
for i in range(1, n + 1):
    s += f(uniform(xp, xk))
s *= dx / n
print(f"Wartość całki wynosi : {s:.3f}")
print()
input("Naciśnij Enter...")

do podrozdziału do strony

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone pod warunkiem podania źródła oraz niepobierania za to pieniędzy.
Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl
Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.

Informacje dodatkowe.

Całkowanie

Metoda Monte Carlo

Wstęp

Algorytm

Specyfikacja problemu

Dane wejściowe

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

Przykładowa implementacja