Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como calcular a relação entre as áreas de um quadrado e um círculo que está inscrito nesse quadrado usando Monte Carlo?

+1 voto
312 visitas
perguntada Jul 24, 2015 em Estatística por danielcajueiro (5,186 pontos)  
Compartilhe

1 Resposta

0 votos
respondida Jul 26, 2015 por danielcajueiro (5,186 pontos)  

A idéia é bem simples. Suponha por exemplo que o círculo tem raio \(r=1\) e, consequentemente, o quadrado tem lado \(l=2\).

A relação entre as áreas pode ser calculada "empiricamente" pela razão entre a frequencia que uma partícula aleatoriamente sorteada com coordenadas \((x,y)\) com distribuição uniforme na caixa \([-1,1]\times [-1,1]\) cai dentro do círculo.

Para resolver explicitamente esse problema, construímos o seguinte código em Python:

import random
import matplotlib.pyplot as plt

def quadradoCirculo(n):
    nc=0
    plt.axis([-1.5, 1.5, -1.5, 1.5])        
    plt.hold(True)
    for i in range(0,n):
        x=random.uniform(-1.0,1.0)        
        y=random.uniform(-1.0,1.0)        
        if(x**2+y**2<1.0):
            plt.plot([x],[y],'ro')
            nc=nc+1.0
        else:
            plt.plot([x],[y],'bo')            
    return nc/n

if __name__ == '__main__':
    n=500
    print quadradoCirculo(n)

Esse pequeno programa tem basicamente uma função que sorteia um vetor de variáveis aleatórias \((x,y)\) com distribuição uniforme na caixa \([-1,1]\times [-1,1]\) e dentro do "if" testa-se apenas se a partícula caiu dentro do círculo ou não. A razão é calculada pelo número de vezes que a partícula caiu dentro do círculo e caiu na caixa inteira \([-1,1]\times [-1,1]\). Veja também a figura abaixo que em uma das simulações eu pintei de vermelho as partículas que cairam dentro do círculo e de azul as outras partículas.

A imagem será apresentada aqui.

Se você rodar esse programa, você vai perceber que a solução converge para a relação entre as áreas que é dada por

\[\frac{\pi r^2}{l^2}=\frac{\pi}{4}.\]

Por preferências pessoais eu escolhi Python (e várias vantagens do Python) para escrever a solução. Você poderia resolver o mesmo problema escolhendo uma outra linguagem de programação.

comentou Ago 23, 2017 por MauroRibeiro (1 ponto)  
Daniel, muito boa sua explicação. Vou deixar aqui minha contribuição, um script em R.
O que enfatizo é:  evitar o  for sempre que possível.

simFunc0 <- function(N){
  contador = 0
  x=numeric(N)
  y=numeric(N)
  for(i in (1:N)){
    x[i] = runif(1,-1,1)
    y[i] = runif(1,-1,1)
    if((x[i])^2+(y[i])^2<=1 ){
      contador = contador + 1  
    }
    result0 = contador/N  
  }
  print(result0)
}

simFunc1 <- function(N){
  x = runif(N,-1,1)
  y = runif(N,-1,1)
  z = cbind(x,y)
  logic  = z[,1]^2+z[,2]^2 <= 1
  result1 = sum(logic)/N
  print(result1)
  }

> system.time(simFunc0(9999999))
[1] 0.7853498
  usuário   sistema decorrido
    36.91      0.20     37.14
> system.time(simFunc1(9999999))
[1] 0.7853512
  usuário   sistema decorrido
     1.25      0.17      1.42
> pi/4
[1] 0.7853982
comentou Ago 23, 2017 por danielcajueiro (5,186 pontos)  
Ola Mauro! Que legal! Mas por que vc não adiciona seu comentario como uma nova resposta para essa pergunta?
comentou Ago 23, 2017 por MauroRibeiro (1 ponto)  
Olá Daniel, ainda estou aprendendo a usar o PRorum ....
...