Liczba pi w Monte Carlo

Z pozoru temat tego artykułu może nie mieć sensu. Jak przecież powiązać ze sobą najbardziej znaną stałą matematyczną z dzielnicą Monako? Tłumaczę: nie chodzi tu o miejscowość, ale o metodę symulowania, która może pomagać w obliczaniu przybliżonych wartości dla procesów, które trudno byłoby wyliczyć sposobem tradycyjnym. Jak się to robi? Pokażę to, wyliczając przybliżenie liczby pi w sposób, jakiego nikt w szkole raczej nie widział.

Metoda Monte Carlo

Za jej twórcę uznaje się Polaka – Stanisława Ulama, który pracował również przy projekcie Manhattan. Jak już napisałem, została ona zaprojektowana w celu przyspieszenia obliczeń, które wykonywane tradycyjnie, trwałyby bardzo długo. Korzystając z mocy obliczeniowych komputerów, możemy określać wyniki, zwyczajnie poprzez dokonywanie symulacji danych doświadczeń.

Za przykład może posłużyć prosty skrypt, który losujący literę o lub r (czyli symulację rzutu monetą). Dokonując wielu powtórzeń (np. 1 000 000) zaobserwujemy, że około połowę wyników będą stanowiły orły, a połowę reszki. Stąd możemy stwierdzić, że prawdopodobieństwo wyrzucenia którejkolwiek strony monety będzie wynosiło $\frac{1}{2}$.

Ciekawostka: Nazwa metody Monte Carlo pochodzi od dzielnicy Monako znanej z hazardu. Kości, black jack, poker, czy ruletka mają w końcu wiele wspólnego z losowością.

Gdzie ta liczba pi?

Ale jak się ma do tego wspomniana wcześniej liczba $\pi$? No cóż – korzystając z metody Monte Carlo spróbujemy dokonać jej przybliżenia i zobrazujemy to na diagramach, by lepiej można było poznać działanie metody Monte Carlo i przy okazji zobaczyć inny sposób na wyliczenie $\pi$. Za narzędzie posłuży nam język Python (i jego biblioteka PyPlot).

Teraz trochę matmy

Jak wiadomo ze szkoły, pole koła wyraża się jako $P = \pi r^2$. Lekko przekształcając ten wzór możemy wyliczyć, że liczba $\pi$ będzie równa

$\pi = \frac{P}{r^2}$

A dodatkowo obierając za promień liczbę $r=1$, zauważymy, że

$\pi = P$

Mogliśmy oczywiście inaczej zdefiniować promień, ale z takim chyba będzie nam najwygodniej się liczyło.

Troszkę bardziej zaawansowaną wiedzą jest już równanie koła (obecnie to poziom rozszerzony w szkole ponadgimnazjalnej). Zgodnie z nim, do koła o środku $(0,0)$ i naszym promieniu $1$ należą punkty $(x,y)$ spełniające nierówność $x^2 + y^2 \leq 1$.

To koło będzie wyglądało właśnie tak:

Jak widać koło zamyka się (jest wpisane) w kwadrat o boku długości 2 i wierzchołkach $(-1,-1)$, $(1,-1)$, $(1,1)$ i $(-1,1)$.

Losowanie czas zacząć

Program, który stworzyłem będzie działał w następujący sposób:

  • wylosuje punkt należący do kwadratu opisanego na naszym kole,
  • sprawdzi, czy ten punkt spełnia nierówność koła (czyli czy należy do koła),
    • jeśli tak, doda ten punkt do zbioru punktów koła ($P_k$)
    • jeśli nie, doda ten punkt do zbioru punktów spoza koła.

W ten sposób będziemy dokładnie wiedzieli, ile udało nam się wylosować punktów należących do koła, a ile takich, które tam nie pasują.

Zobaczmy, co się stanie dla:

A. Stu losowań:


B. Tysiąca losowań:

C. Dziesięciu tysięcy losowań:

Nietrudno zauważyć, że różowe punkty zaczynają wypełniać kształt koła, a łącznie z niebieskimi będą tworzyły cały kwadrat – oczywiście tylko pozornie, gdyż punkty są nieskończenie małe.

Czym będzie zatem pole koła? W przybliżeniu będzie to stosunek wylosowanych punktów wewnątrz koła – $P_k$, do liczby losowań i przemnożony przez pole kwadratu, czyli 4. Zapisując to bardziej elegancko i pamiętając, że pole naszego koła to pi, otrzymamy, że:

$\pi = \frac{P_k}{n} * 4$

Pora na miliony

Obliczmy w takim razie przybliżenie $\pi$, korzystając ze skryptu napisanego w języku Python. Dla dobrego przybliżenia użyjemy miliona prób i oto co dostaniemy:

Otrzymaliśmy wynik $\pi = 3.141808$, czyli całkiem nieźle. Mając odrobinę więcej czasu można pokusić się o większe liczby prób – miliard, bilion … – każda jednak będzie trwała 1000 razy dłużej. Warto też zaznaczyć, że generator liczb pseudolosowych może po pewnym czasie zacząć powtarzać liczby, co będzie fałszowało wyniki.

Podsumowanie

Jak widać metoda Monte Carlo to całkiem ciekawy sposób, pozwalający na ominięcie niepotrzebnych przeszkód w postaci skomplikowanych obliczeń i zastąpieniem ich wielokrotnymi próbami. Jest ona używana na giełdzie, w przewidywaniu ile energii zwróci farma wiatraków, czy w tworzeniu promieni światła w grach komputerowych.

Bonus: Skrypt Pythona

A gdyby ktoś chciał samodzielnie pobawić się symulacjami, to publikuję mój skrypt poniżej:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
from random import random
from random import seed
import matplotlib.pyplot as plt

x_in = []
y_in = []

x_out = []
y_out =[]

n = 1000000 #Tu trzeba wpisać liczbę prób do symulacji

for i in range(n):
  seed()
  x_tmp = 2*random() -1
  y_tmp = 2*random() -1

  if x_tmp**2 + y_tmp**2 <= 1:
    x_in.append(x_tmp)
    y_in.append(y_tmp)
  else:
    x_out.append(x_tmp)
    y_out.append(y_tmp)

pi = len(x_in) / n * 4
label = 'Pole koła to ' + str(pi)
#Tu robi się wykres
plt.figure(num=None, figsize=(10, 10), dpi=72)
plt.axis([-1, 1, -1, 1])
plt.grid(color='#707070', linestyle='-')
plt.scatter(x_in, y_in, c="#fe5fc8", alpha=0.5)
plt.scatter(x_out, y_out, c="#7a00d8", alpha=0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.title(label)
plt.show()

Nauczyciel matematyki i przedmiotów informatycznych w Dwujęzycznych Szkołach im. Władysława Kopalińskiego w Bielsku-Białej. Autor wielu publikacji poświęconych nauce podstaw programowania, w tym „Scratch. Nauka programowania przez zabawę”, wydanej drukiem przez wydawnictwo Komputer Świat. Entuzjasta wykorzystywania nowych nowoczesnych narzędzi w czasie lekcji. Prywatnie miłośnik koszykówki i aktywnego wykorzystywania czasu.

Dodaj komentarz

Twój adres email nie zostanie opublikowany. Pola, których wypełnienie jest wymagane, są oznaczone symbolem *

This site uses Akismet to reduce spam. Learn how your comment data is processed.