Càlcul de pi mitjançant el mètode de Monte Carlo

(actualisé le ) per  Sàdur

El mètode de Monte Carlo és un mètode estadístic que permet aproximar expressions matemàtiques complexes. Aquest mètode usa experiments aleatoris —normalment simulacions fetes amb ordinadors— per a avaluar el valor desitjat i s’usa principalment en casos en què no és possible (o és molt costós) fer un càlcul determinista (per exemple en física de fluids, on hi ha un nombre de variables acoblades).

Ací, exposarem un exemple senzill d’ús del mètode de Monte Carlo: el càlcul del valor de π.

Imaginem un cercle de radi 1 (fa igual que siga 1 metre que un centímetre) inscrit dins d’un quadrat (que, lògicament, tindrà el costat de longitud 2).

Circunferència de radi unitat inscrita en un quadrat.

Les àrees respectives, com sabem, vénen donades per les expressions
L'àrea d'un cercle és igual a pi pel quadrat del seu radi.
L'àrea d'un quadrat és igual al quadrat de la longitud del costat.

Imaginem també que hi llancem dards amb els ulls tancats (per tal que els llançaments siguen completament aleatoris).

En tal cas, la probabilitat que un llançament concret caiga dins del cercle és
Probabilitat que el dard caiga dins del cercle inscrit dins del quadrat.
que, en el nostre cas, és
Probabilitat que un dard caiga dins del cercle.

De manera que, si podem computar aquesta probabilitat, podrem calcular el valor de π multiplicant la dita probabilitat per quatre. Per exemple, si haguérem calculat que la probabilitat és de 0,8 [1], el valor de ens donaria 4·0,8=3,2.

En comptes de fer els llançaments realment, podem simular-los mitjançant un programa d’ordinador [2].

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Aquest programa s'ha escrit per a ser executat amb python2.x
import random, time

nre_llancaments = long(raw_input("Introduïu el nombre de llançaments: "))
nre_dianes = 0L
inici = time.time()

for i in xrange(nre_llancaments):
        x = random.random()
        y = random.random()
        if i % 100000 == 0:
                print ".",
        if (x**2 + y**2) <= 1:
                nre_dianes += 1

final = time.time()

print
print "RESULTATS"
print "Temps invertit:", final - inici, "segons."
print "Llançaments:", nre_llancaments
print "Dianes:", nre_dianes
print "Probabilitat:", float(nre_dianes)/nre_llancaments
print "Estimació de π:", float(4*nre_dianes)/nre_llancaments

D’acord amb la llei dels grans nombres, la probabilitat teòrica s’assoleix quan el nombre de llançaments tendisca a infinit. Per tant, el càlcul serà tant més precís quant més llançaments es facen.

Hem fet l’experiment amb 1000000000 llançaments. Les dades obtingudes cada 100000000 de llançaments vénen reflectides en la gràfica següent.

Evolució del valor de l'estimació de π en funció del nombre de llançaments.

Com s’hi pot observar, el valor calculat (vermell) s’acosta al valor de pi (blau) a mesura que augmenta el nombre de repeticions.

Una característica del mètode és que si fem l’experiment més d’una vegada, el resultat és lleugerament diferent en cada cas.

Amb 1000000000 repeticions el resultat obtingut ha estat:

  • Llançaments: 1000000000
  • Dianes: 785412122
  • Probabilitat calculada: 0.785412122
  • Estimació de π: 3.141648488

El valor obtingut té un error menor que una deumil·lèsima (comparat amb el valor estàndard de pi: 3,14159265).

Notes

[1Suposant que haguérem fet 10 llançaments i 8 hagueren caigut dins del cercle.

[2Per a executar aquest programa heu d’instal·lar python en l’ordinador. Si preferiu veure’n una versió web podeu visitar cálculo de pi o pi by Monte Carlo method.