Método de Monte Carlo e transformações: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
Sem resumo de edição
 
Linha 77: Linha 77:
<pre>
<pre>
#Cálculo da área de um círculo
#Cálculo da área de um círculo
import matplotlib.pyplot as plt #Plotar gráfico
import matplotlib.pyplot as plt   #Plotar gráfico
import numpy as np             #Funções matemáticas
import numpy as np               #Funções matemáticas
import random                   #Números aleatórios
import random                     #Números aleatórios


r=1 #Raio do círculo
r=1                               #Raio do círculo
Np=1000000 #Número de pontos
Np=1000000                       #Número de pontos
hitx=[];hity=[];missx=[];missy=[] #Resultados
hitx=[];hity=[];missx=[];missy=[] #Gurdar resultados
for n in range(Np):
for n in range(Np):
   x=random.random()*2-1           #Pontos aleatórios gerados entre -1 e +1  
   x=random.random()*2-1           #Pontos aleatórios gerados entre -1 e +1  
   y=random.random()*2-1           #Pontos aleatórios gerados entre -1 e +1
   y=random.random()*2-1           #Pontos aleatórios gerados entre -1 e +1
   if(np.sqrt(x**2+y**2)<=r):       #Se o ponto caiu dentro do círculo
   if(np.sqrt(x**2+y**2)<=r):     #Se o ponto caiu dentro do círculo
     hitx.append(x);hity.append(y)
     hitx.append(x);hity.append(y)
   else:
   else:
     missx.append(x);missy.append(y)
     missx.append(x);missy.append(y)


#Plotar curvas
#Plotar o círculo
x=np.arange(-1,1,1E-5)
x=np.arange(-1,1,1E-5)
plt.plot(x,-np.sqrt(r-x**2),'-k')
plt.plot(x,-np.sqrt(r-x**2),'-k')
Linha 100: Linha 100:
plt.plot(missx,missy,'ob')
plt.plot(missx,missy,'ob')


a_ret=2*2 #Área do retângulo
#Cálcular a área
a_cir=a_ret*len(hitx)/Np     #A_ret *(hit/N)
a_ret=2*2                         #Área do retângulo: geramos pontos entr -1 e +1 na duas dimensões, então é um quadrado de lado 2
a_cir=a_ret*len(hitx)/Np         #A_ret *(hit/N)
print('A área do círculo estimada é: {:.2f}'.format(a_cir))
print('A área do círculo estimada é: {:.2f}'.format(a_cir))
print('A área do círculo exata    é: {:.2f}'.format(np.pi*r**2))
print('A área do círculo exata    é: {:.2f}'.format(np.pi*r**2))
</pre>
</pre>

Edição atual tal como às 17h09min de 29 de abril de 2022

Transformação linear

Sorteando um número aleatório então fazemos uma transformação para obter um número . Isto é, obtemos a seguinte transformação da seguinte forma:

Se todos os números entre e tinham igual probabilidade de serem sorteados, após a transformação todos os números entre e também possuem igual probabilidade, pois varia com de forma linear, isto é, a distribuição uniforme de números é mantida.Por sua vez, a distribuição uniforme significa que a probabiliade de obter um número entre e é dada pela função densidade de probabilidade da seguinte forma:

Sendo constante, temos que:

Pela normalização. Mas se ampliarmos o intervalo dos números possíveis para entre e :

Então agora . Isto é distribuição de probabilidade continua constante, mas com uma menor probabilidade de sortear um número qualquer, quando em comparação de sortear um qualquer.

Transformação não-linear

O mesmo não ocorre com uma transformação não linear. Por exemplo se , derivando temos que:

Diferente do caso anterior que tínhamos apenas a transformação linear . Podemos ver ainda usando a própria definição de derivada:

E sendo um diferencial então , logo . Agora a distribuição de probabilidade é alterada com a transformação. Para uma transformação qualquer, como a probabilidade se conserva, ainda temos:

Considrando que tem inversa, então então:

Sendo então reescrevendo , logo Então temos:

E sendo nossa ditribuição em uniforme com como vimos anteriormente, ficamos com:

Para a transformação . O mais comum é que saibamos a distribuição de probabilidade que queremos, e uma vez que:

E integramos então para encontrar .

Método da rejeição

Nem sempre o desejado é fácil de definir matematicamente. O método da rejeição é um método rústico para obtermos .

  • Desenhar a função desejada dentro dos limites e .
  • Geramos um ponto aleatório , se estiver abaixo da curva desejada é aceito.

Se gerarmos pontos aleatórios em grande quantidade, a razão de pontos aceitos para cada em relação à todos o pontos aleatórios gerados neste , nos dá uma estimativa de neste ponto, em relação do . Isto é, se para um qualquer, metade dos pontos gerado aleatoriamente foram válidos, então .

Cálulo de integrais definidas

Em uma ideia bastante análoga à anterior, aqui utilizamos a ideia que a integral definida é cálculo da área sob a curva. Definindo novamente limites e no qual a região que queremos calcular está contida, geramos uma série de pontos aleatórios, se gerado pontos suficientes, a razão de pontos que foram gerados dentro da área que queremos calcular em relação ao total de pontos gerados, será a mesma razão da área que queremos calcular em relação à área total definida pelos limites.

#Cálculo da área de um círculo
import matplotlib.pyplot as plt   #Plotar gráfico
import numpy as np                #Funções matemáticas
import random                     #Números aleatórios

r=1                               #Raio do círculo
Np=1000000                        #Número de pontos
hitx=[];hity=[];missx=[];missy=[] #Gurdar resultados
for n in range(Np):
  x=random.random()*2-1           #Pontos aleatórios gerados entre -1 e +1 
  y=random.random()*2-1           #Pontos aleatórios gerados entre -1 e +1
  if(np.sqrt(x**2+y**2)<=r):      #Se o ponto caiu dentro do círculo
    hitx.append(x);hity.append(y)
  else:
    missx.append(x);missy.append(y)

#Plotar o círculo
x=np.arange(-1,1,1E-5)
plt.plot(x,-np.sqrt(r-x**2),'-k')
plt.plot(x,np.sqrt(r-x**2),'-k')
#Plotar pontos
plt.plot(hitx,hity,'or')
plt.plot(missx,missy,'ob')

#Cálcular a área
a_ret=2*2                         #Área do retângulo: geramos pontos entr -1 e +1 na duas dimensões, então é um quadrado de lado 2
a_cir=a_ret*len(hitx)/Np          #A_ret *(hit/N)
print('A área do círculo estimada é: {:.2f}'.format(a_cir))
print('A área do círculo exata    é: {:.2f}'.format(np.pi*r**2))