Números Aleatórios

Introdução

Existem na natureza processos de caráter estocástico, cujo resultado não pode ser predito a priori. Os exemplos típicos são o lançamento de uma moeda, o jogo dos dados ou a roleta. Em cada um desses casos sabemos quais são os resultados possíveis, porem não ha como predizer o resultado da próxima jogada. Entretanto podemos sim dizer qual será resultado em termos estatísticos de um grande número de jogadas: se lançarmos 1000 vezes uma moeda aproximadamente 50% das vezes o resultado será cara, e a percentagem será cada vez mais próxima de 50% quanto maior o número de lançamentos. O mesmo pode ser dito da roleta: num grande número de jogadas aproximadamente 1/37 das vezes terá dado o 17 por exemplo. Embora não seja suficiente para ganhar no jogo é um resultado previsível do ponto de vista estatístico. De uma outra forma, se nos fosse permitido colocar uma ficha nos pares e no zero ao mesmo tempo repetindo a jogada um grande número de vezes, teríamos lucro no final das contas.

Por outra parte existem na física processos cuja natureza é em principio determinística, cujo resultado é, em princípio, previsível conhecendo as condições inicias, mas que na prática resulta impossível como por exemplo a exata trajetória das moléculas de um gás.

Mesmo assim poderíamos integrar as equações de Newton de N partículas (moléculas ou planetas) com um dos algoritmos estáveis usados na unidade anterior: o Euler-Cromer ou de preferência o Verlet. O principio é o mesmo que no caso de uma única partícula como foi o caso do oscilador harmônico unidimensional. No caso de N partículas existem algumas complicações adicionais pois é necessário somar a força que as N-1 partículas fazem em uma dada partícula. À integração numérica das equações de movimento de N partículas é dada o nome de Dinâmica Molecular. O método é muito usado é viável num PC até 100000 ou 1 milhão de partículas. Porem depois de um determinado tempo a trajetória obtida para uma dada partícula pode estar muito afastada da trajetória teórica. No fundo acontece algo parecido com o problema da moeda. Ele é em princípio um problema determinístico mas na prática resulta randômico. Se diz que são problema sensíveis as condições iniciais: uma pequena variação na posição ou velocidade inicial da moeda se traduz numa trajetória muito diferente depois de um certo tempo.

Mas no caso de um gás, no fim não interessa pois nada importa a posição exata de uma dada molécula, o que importa são as propriedades estatísticas do sistema com um todo. Fazendo o análogo com o problema da moeda, é saber que proporção de cara ou coroa teremos depois de N lançamentos. No gás as medidas que interessam são a pressão, energia, temperatura, etc, propriedades médias do sistema em equilíbrio.

Dependendo das propriedades de interesse e dada essa caraterística randômica do movimento individual complexo de uma molécula qualquer, podemos adotar um modelo mais simples para simular o comportamento do sistema que o método de Dinâmica Molecular, relativamente difícil de programar. Tomemos o exemplo de um gás que está inicialmente contido numa metade de uma caixa, que está dividida ao meio por uma parede com um fim pequeno furo pelo qual pode passar uma molécula por vez. Sabemos que depois de um tempo o sistema alcançará o equilíbrio e que o número médio de partículas em cada metade da caixa será N/2.

Uma maneira simples de simular esse processo e supor que as partículas não interagem entre elas e que a probabilidade (por unidade de tempo) de uma dada partícula trocar de lado e a mesma para todas independente do lado da caixa na qual se encontre. Podemos implementar esse modelo escolhendo aleatoriamente uma partícula e trocá-la de lado. Acompanhando este processo no tempo deveríamos observar a evolução do sistema ao equilíbrio.

Números aleatórios

Para implementar esse modelo precisamos números aleatórios. Os verdadeiros números aleatórios são aqueles que resultam de um experimento estocástico, como os exemplos citados ao começo: moeda, dados, roleta. Mas para poder usar em cálculos no computador devemos achar um método mais rápido de gerar um seqüência de números aleatórios, deve ser o própio computador quem forneça essa seqüência. Isso é possível como veremos em seguida, porem eles não são estritamente aleatórios, eles parecem não obedecer a nenhum padrão, nenhuma seqüência lógica parece estar por trás desses números, mas na verdade se trata de uma seqüência completamente previsível, só para quem conhece a regra de geração da seqüência, ou seja ninguém poderia acertar o próximo número a menos que conheça a receita. Ou seja, que podem se comportar como números aleatórios tão bons quanto os verdadeiros. Por isso são as vezes chamados de pseudo-aleatórios.

Uma maneira de gerar uma seqüência assim é por meio da operação de congruência (resto da divisão inteira):

xn+1 = resto { (a xn + b)/m }

onde a, b, m e os xn são todos números inteiros. Dependendo da escolha dos inteiros a, b e m, a sequencia de números {x1, x2, ..., xn} é uma sequencia randómica dos números entre 1 e m. A melhor maneira de entender como isso funciona é com um exemplo: Seja a=4, b=3 e m=1, e seja o x inicial x0= 2

Em python, o método da congruência pode ser implementado da seguinte forma:

# gerador de num aleatorios de congruencia
def randi():
    global x
    a = 1029; b= 221591; m = 1048576
    x = (a*x + b) % m
    return x

# entrar aqui um interiro como semente inicial
x = int(input('semente (entrar inteiro)? '))
for i in range(10):
    n = randi()
    print(i,n)

Exemplos

Como uma solução nativa em Python temos a biblioteca random, para obter um número aleatório entre 0 e 1 podemos utilizar random(), ou randint(a,b) para inteiros entre (incluindo) a e b:

import random
print("Número aleatório: ",random.random()) # Frações no conjunto [0,1)
print("Roleta:",random.randint(0,37))       # Inteiros no conjunto [0,37]
print("Dado:",int(random.random()*6)+1)     # Inteiros no conjunto [1,6]
print("Moeda: ", 2*random.randint(0,1)-1)   # Inteiros no conjunto {-1,1}

O caminhante aleatório, pode ser facilmente implementado, tanto para se ter a mesma probabilidade de se mover em qualquer direção quanto para probabilidades diferentes.

import random
import matplotlib.pyplot as plt

n = 1000000 #Número de passos
p =0.5    #probabilidade de ir para o lado direito
x=[0]     #Posição inicial
c=0       #contador
for j in range(n):                #Com n passos em cada repetição
  dx = (+1) if (random.random()<=p) else (-1)
  x.append(x[j]+dx)               #Atualizamos a posição
  c = (c+1) if (dx>0) else (c)    #Atualizamos o contador
plt.plot(x)       
print("A particula se moveu para a direita {:.2f}% das vezes".format(100*c/n))

Também podemos implementar uma "apuração" de uma eleição. Por exemplo supondo que determinado candidato tem p% de chance de receber um voto, então:

import random
import matplotlib.pyplot as plt

p =0.5    #probabilidade do candidato ser votado
v=[0]     #Quantidade de votos em %
T=0       #Quantidade de votos totais
for i in range(1,1601): #Quantidade de candidatos
  x=(1) if (random.random()<=p) else (0)
  T+=x
  v.append(T/i)
plt.plot(v)       

Gás na caixa

E por fim, podemos implementar a simulação do gás na caixa em Python da seguinte forma:

import random
import matplotlib.pyplot as plt

N = 1000     #Quantidade de partículas
cx =N*[0]    #Consideramos 0 à esquerda, e 1 à direita
np = 10000   #Número de pasos
t = [sum(cx)]#Quantidade de partícula a direita

for i in range(np):
  j=random.randint(0,N-1) #Pegamos um índice
  cx[j] = (1) if (cx[j]==0) else (0)
  t.append(sum(cx))

plt.plot(t)

A continução temos um exemplo de código FORTRAN também implementando a simulação do exemplo do gas na caixa. No exemplo simulamos 1000 partículas e fica como exercício decidir que variáveis imprimir e como calculá-las. Também se é necessário imprimir a cada passo da simulação.

!Aula5: Programa Gas, simulação com 1000 partículas
Program Gas
Implicit None
Real*8 Randi
Integer N, i, j, t, x
Parameter (N = 1000)
Integer Box(N)

Read(*,*) t, x  !t:tempo de simulação, x:semente de Randi

Do i = 1, N
   Box(i) = -1
End Do

Do j = 1, t
   i = N*Randi(x)+1
   Box(i) = -Box(i)
   Print*, j, 'deixo/vc(exercício 2)'
End Do

End Program Gas

Observações sobre FORTRAN:

  1. Parameter é para assignar (em tempo de compilação) um valor fixo a uma variável que não pode mudar durante a execução do programa, evitando o erro decorrente de uma alteração involuntária por distração, etc.
  2. O segundo novo elemento é o vector A(10) ou B(5,20).
     No caso que nos ocupa é Box(N); a utilidade dele fica evidente no própio programa: de uma vez só definimos uma variável a qual, por meio do indice, permite guardar o estado de um grande número de variáveis. Neste caso são as posições das 1000 partículas.

Exercícios:

  1. Graficar a evolução temporal das partículas a esquerda e a direita da parede (ver figura acima).
  2. Observe o "tempo" necessário para entrar em regime estacionário. Como varia esse tempo se mudarmos o número de partículas. Você pode estabelecer uma relação entre t e N?
  3. Implemente uma maneira de calcular a média de partículas num lado da caixa (precissa os dois?) e veja como varia no tempo.
  4. Depende do intervalo usado para calcular a média?
  5. Observe a variação no número de partículas no regime estacionário. Implemente uma forma de quantificar essa variação no seu programa. Ou seja a desvio cuadrático médio ou flutuação
  6. Quantifique a dependencia do desvio no regime estacionário com o número de partículas

Campo médio

Escrevendo então a quantidade de partículas no lado direto da caixa em um tempo ${\textstyle t}$ como ${\textstyle n_{+}\left(t\right)}$ então escrevemos:

$$n_{+}\left(t+1\right)=n_{+}\left(t\right)-p_{+}\left(t\right)+p_{-}\left(t\right)$$

Onde ${\textstyle p_{+}\left(t\right)}$ é a probabilidade de encontrar uma partícula no lado direito no instante ${\textstyle t}$. Ou seja, uma vez que qualquer partícula pode ser selecionada com a mesma probabilidade, então:

$$p_{+}\left(t\right)={\frac {n_{+}\left(t\right)}{N}},\qquad p_{+}\left(t\right)=1-{\frac {n_{+}\left(t\right)}{N}},$$

Substituindo:

$${\begin{aligned}n_{+}\left(t+1\right)&=n_{+}\left(t\right)-{\frac {n_{+}\left(t\right)}{N}}+1-{\frac {n_{+}\left(t\right)}{N}}\\&=n_{+}\left(t\right)+1-2{\frac {n_{+}\left(t\right)}{N}}\\n_{+}\left(t+1\right)-&n_{+}\left(t\right)=1-{\frac {2n_{+}\left(t\right)}{N}}\end{aligned}}$$

Esta é a taxa de variação em um passo ${\textstyle \Delta t=1}$. Então podemos obter o seguinte modelo de campo médio:

$${\frac {d}{dt}}n_{+}\left(t\right)=1-{\frac {2n_{+}\left(t\right)}{N}}\longrightarrow n_{+}\left(t\right)={\frac {N}{2}}\left(1-e^{-{\frac {2t}{N}}}\right)$$

import random
import matplotlib.pyplot as plt
import numpy as np

N = 10000     #Quantidade de partículas
cx =N*[0]     #Consideramos 0 à esquerda, e 1 à direita
al = [0]      #Solução analítica para a quantidade de partículas à direita
NP = 30000    #Número de pasos
t = [sum(cx)] #Quantidade de partícula a direita

for i in range(NP):
  j=random.randint(0,N-1)               #Pegamos um índice
  cx[j] = (1) if (cx[j]==0) else (0)    #Números aleatórios
  t.append(sum(cx))
  al.append((1-np.exp(-2*(i+1)/N))*N/2) #Campo médio

plt.plot(t,label="Número aleatórios")
plt.plot(al,label="Campo-médio")
plt.legend()                 

O campo médio nos dá a média de várias trajetórias no sistema, isto é, se somássemos todas trajetórias do sistema, o que obtemos representa a média de infinitas trajetórias. Logo o campo-médio não só representa a média de várias trajetórias com baixas partículas, mas a medida que aumentamos o número de partículas, ele se aproxima desta trajetória.