domingo , 5 julio 2020
Estima Pi con random() en Python
Estima Pi con random() en Python

¿Podemos calcular $\pi$ con números aleatorios?

En este tutorial demostraremos una manera de aproximar el valor de $\pi$ empleando el método de Montecarlo. Como es bien sabido $\pi$ representa el cociente del perímetro de una circunferencia dividido entre su diámetro y es uno de los números más importantes de la matemática ya que aparece espontáneamente en multitud de ecuaciones y cálculos. A menudo $\pi$ aparece como resultado de ecuaciones que en apariencia no tienen relación alguna con la circunferencia como por ejemplo: $$\sum_{n=1}^\infty \frac{1}{n^2}=\frac{\pi^2}{6}$$ lo que convierte a $\pi$ en un número enigmático del que sabemos que no solo no es racional, sino que tampoco es irracional, perteneciendo a la categoría de números conocidos como trascendentes (números que no pueden expresarse como la raíz de una ecuación algebraica con coeficientes enteros). Aunque parezca extraño, se sabe que la recta numérica está poblada por más números trascendentes que racionales e irracionales, sin embargo hasta ahora nuestro conocimiento de los números trascendentes es limitado ya que son muy pocos los ejemplos de números que sabemos a ciencia cierta que son trascendentes.

Método de Montecarlo

Llamamos método de Montecarlo a la derivación de soluciones aproximadas a problemas numéricos basándonos en la generación de grandes cantidades de números pseudoaleatorios. El método se desarrollo originalmente durante la segunda guerra mundial como parte del proyecto Manhattan para calcular la difusión de neutrones en materiales fisibles, cálculo que era necesario para el diseño de la bomba atómica. Con el advenimiento de la computación moderna, esta herramienta a encontrado aplicaciones mucho más pacíficas dentro de las ciencias y la ingeniería y es parte integral de muchas herramientas de cálculo que se utilizan para medir volúmenes de superficies complejas o aproximar propiedades de materiales.

El caso de $\pi$.

Consideremos el caso de una circunferencia de radio 1 inscrita en un cuadrado como la que se muestra en la siguiente imagen. Si consideramos que las áreas del círculo y del cuadrado están dadas por $A_\bigcirc=\pi r^2$ y $A_\square=4$ respectivamente, entonces el cociente de ambas áreas nos dice que: $$\frac{A_\bigcirc}{A_\square}=\frac{\pi r^2}{4}$$ despejando $\pi$ de esta ecuación y considerando que $r=1$ encontramos: $$\pi=4 \frac{A_\bigcirc}{A_\square}$$ es decir, $\pi$ puede estimarse como el cociente de las áreas del cuadrado y del círculo multiplicado por una constante.

Puntos aleatorios lanzados al círculo unitario tienen una probabilidad proporcional al área de caer en el interior del círculo.
Si lanzamos puntos al azar sobre el cuadrado gris, la probabilidad de que caigan dentro del círculo es proporcional al área del círculo.

Ahora bien, determinar al área del círculo no sería en si mismo una tarea trivial, pero el cociente del área del círculo respecto al área del cuadrado es algo que podemos estimar con la precisión que deseemos empleando los trucos adecuados. Supongamos que lanzamos puntos aleatorios sobre toda el área del cuadrado, y que la distribución de los puntos es uniforme, es decir, cualquier parte del cuadrado tiene la misma probabilidad de recibir un punto. Si este es el caso, podemos asumir que la probabilidad de que los puntos caigan en el interior del círculo es proporcional al área que cubre el círculo, de modo que cuando el número de puntos $num\ puntos \to \infty$, el cociente de los puntos dentro del círculo dividido entre los puntos en toda el área del cuadrado converge al cociente de las áreas: $$\frac{A_\bigcirc}{A_\square}=\lim_{num\ puntos \to \infty} \frac{Puntos_\bigcirc}{Puntos_\square}$$ de donde el valor de $\pi$ puede expresarse como: $$\pi= \lim_{num\ puntos \to \infty} 4 \frac{Puntos_\bigcirc}{Puntos_\square} $$

En la práctica, calcular el valor exacto del límite es imposible, pero podemos aproximarnos al valor límite tanto como deseemos lanzando los suficientes puntos.

Implementación

Para implementarlo, comenzamos por definir una función punto() que genera las coordenadas de un único punto dentro del cuadrado y determina si tales coordenadas se encuentran dentro o fuera del círculo de radio 1. La función retorna una tupla con True o False indicando si las coordenadas están dentro o fuera del círculo además de las coordenadas del punto. Aunque el cálculo de $\pi$ solo requiere saber si el punto cayo dentro de la circunferencia, en este ejemplo devolvemos las coordenadas para poder graficar los puntos y servirnos de guía para averiguar que tan completa fue nuestra cobertura del círculo.

import numpy as np
import matplotlib.pyplot as plt

def punto():
    coordenadas=[np.random.uniform(-1,1), np.random.uniform(-1,1)]
    return coordenadas[0]**2+coordenadas[1]**2<1, coordenadas

Con esto podemos definir una nueva función que llamamos aproxPi(n) que recibe como argumento un valor n, que representa el número de puntos totales a lanzar sobre el cuadrado. La función calcula el valor de $\pi$ haciendo varias llamadas a punto() y entrega varias listas con las coordenadas de los puntos que cayeron dentro o fuera del círculo. En un principio las listas de coordenadas se inicializan como listas vacías y se van poblando punto por punto con append agregando cada coordenada al final de la lista que corresponda. Al final la función retorna un valor estimado de $\pi$ basado en el número de puntos que cayeron dentro del círculo junto con las listas de coordenadas.

def aproxPi(n):
    dentro=0
    Xdentro=[]
    Ydentro=[]
    Xfuera=[]
    Yfuera=[]
    for i in range(n):
        estado, coordenadas = punto()
        if estado:
            dentro+=1
            Xdentro.append(coordenadas[0])
            Ydentro.append(coordenadas[1])
        else:
            Xfuera.append(coordenadas[0])
            Yfuera.append(coordenadas[1])
    return 4*dentro/n, Xdentro, Ydentro, Xfuera, Yfuera

Ahora podemos graficar directamente la respuesta de esta función, en este ejemplo estaremos usando 10,000 puntos. Dentro de Jupyter podemos usar la combinación de teclas \pi + tab para generar el simbolo griego de $\pi$. Lo mismo se aplica a cualquier otra letra griega, por ejemplo, escribir \alpha + tab genera $\alpha$. Adviértase que estos solo son símbolos, para Python $\pi$ es una letra como cualquier otra que podemos usar para designar una variable.

n=10000
π=aproxPi(n)

fig=plt.figure(figsize=(7,7))
ax=fig.add_subplot(111)
ax.set_aspect('equal')

circulo = plt.Circle((0, 0), 1, color='0.9', fill=True, zorder=1)
ax.add_artist(circulo)
ax.scatter(π[1],π[2], s=0.2, c='r', marker="o", zorder=2)
ax.scatter(π[3],π[4], s=0.2, c='b', marker="o", zorder=3)

plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.title(f'$n={n:0>9,}$, $\pi \sim {π[0]:0<7.5}$')
plt.show()

Aquí el comando plt.figure(figsize=(7,7)) crea una figura a la que llamamos fig que mide 7×7 pulgadas. Con fig.add_subplot(111) creamos unos ejes coordenados que llamamos ax sobre los que vamos a graficar, ax.set_aspect('equal') le dice a Python que emplee la misma relación de aspecto para ambos ejes. Aquí usamos circulo para designar a una circunferencia que nos ayuda a visualizar los límites para aprobar o rechazar los puntos. La circunferencia está centrada en (0, 0) tiene un radio de 1, un color gris claro dado por color='0.9', la circunferencia está coloreada por dentro (fill=True) y es el primer objeto en agregarse a la gráfica (zorder=1), con lo que nos aseguramos que los puntos que vendrán después se grafiquen encima del disco y no al revés. El círculo se añade a los ejes con ax.add_artist.

Para graficar los puntos a los ejes de manera dispersa usamos ax.scatter, usamos un tamaño de punto pequeño (s=0.2) y asignamos un color rojo a los puntos que caen dentro del círculo y el azul a los puntos que caen fuera con c='r' y c='b'. Se asigna un marcador de tipo punto a todos los datos (marker="o") y se asigna un número de zorder mayor a 1 (2 y 3) para que los puntos se grafiquen por encima del círculo.

Usamos plt.xlim(-1, 1) para limitar el área de la gráfica a la zona de interés y usamos una f-string para desplegar como título de la gráfica el valor de $n$ y nuestra aproximación de $\pi$ usando un formateo personalizado para desplegar los números. Como un ejemplo {π[0]:0<7.5} significa que el valor de $\pi$ debe presentarse colocando 0‘s para rellenar cualquier espacio vacío, que los números estarán alineados hacia la izquierda (<) con lo que cualquier 0 que se agregue se colocará a la derecha, que el número debe medir en total 7 caracteres (incluyendo puntos decimales y comas) y que contaremos con 5 elementos de precisión decimal. Con esto, la gráfica resultante es la siguiente:

Gráfico de puntos en el circulo unitario para n = 10000.
Resultado de evaluar $\pi$ con 10,000 puntos.

Y haciendo una composición de las gráficas obtenidas para valores de $n$ entre 10 y 1,000,000 obtenemos:

Evolución de Pi para un número creciente de puntos n=10..1000000.
Composición de los valores de $\pi$ obtenidos para $n=10…1,000,000$.

Resultados

Una de las desventajas de este método es que requiere grandes cantidades de números para arribar a una conclusión. De la animación de arriba nótese que para arribar a las primeras 2 cifras decimales (3.14) se requieren al menos 100,000 puntos y que aún para 1,000,000 de puntos no se ha ganado ninguna otra cifra decimal. Esto convierte al método en un proceso muy intensivo en su uso de recursos computacionales, con lo que resulta claro que no es el método más eficiente, pero si es una opción cuando no tenemos herramientas más precisas a la mano. En el presente ejemplo la serie convergente mencionada al principio de este artículo sería una manera mucho más eficiente de calcular $\pi$, pero el ejercicio nos brinda una idea clara de como se aplica el Método de Montecarlo para la solución de problemas concretos.

Código fuente

Si conoces GitHub, puedes hallar el notebook de Jupyter con los comandos de este ejercicio en el link:

GitHub

También puedes descargarlo directamente de aquí:

Descarga

Acerca de itsgaraet

Físico especializado en aplicaciones ópticas. En la prepa me llamó la atención aprender lenguaje C y desde entonces comencé a usar programación para resolver problemas que lo ameritaran. Después de aprender Fortran, Matlab y LabVIEW para automatizar mediciones, realizar simulaciones o calcular propiedades físicas, me di cuenta de que a menudo solo buscaba pretextos para resolver problemas empleando programación ya que la experiencia de definir un problema y plantear su solución es una experiencia que disfruto bastante. Esto me llevó a intentar profundizar en el campo del software, principalmente en lo relativo a Ciencia de Datos y Desarrollo Web. También soy amante de la fotografía, planear una nueva imagen y terminar de llevarla a la realidad en Photoshop siempre es algo interesante.

Checa también

¿Sabes que significa? - Diagrama de bifurcacion

Aplicación logística y diagrama de bifurcación

Tutorial en Python explicando el significado de la aplicación logística y desarrollando el código para graficar los diagramas de bifurcación de varias funciones.

Deja un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *