Modelamento e simulação de processos

Modelamento

e simulação de

processos 4. Método de Monte Carlo

Prof. Dr. André Carlos Silva

1. INTRODUÇÃO

O Método de Monte Carlo (MMC) é um

método estatístico utilizado em simulações

estocásticas com diversas aplicações em

áreas como a física, matemática e biologia.

O método de Monte Carlo tem sido utilizado

há bastante tempo como forma de obter

aproximações numéricas de funções

complexas.

1. INTRODUÇÃO

Este método tipicamente envolve a geração de observações de alguma distribuição de probabilidades e o uso da amostra obtida para aproximar a função de interesse.

As aplicações mais comuns são em computação numérica para avaliar integrais.

1. INTRODUÇÃO

1. INTRODUÇÃO

A idéia do método é escrever a integral

que se deseja calcular como um valor

esperado.

De acordo com (HAMMERSELEY,1964) o

nome "Monte Carlo" surgiu durante o

projeto Manhattan na Segunda Guerra

Mundial.

1. INTRODUÇÃO

No projeto e de construção da bomba

atómica, Ulam, von Neumann e Fermi

consideraram a possibilidade de utilizar o

método, que envolvia a simulação direta

de problemas de natureza probabilística

relacionados com o coeficiente de

difusão do nêutron em certos materiais.

1. INTRODUÇÃO

Apesar de ter despertado a atenção desses

cientistas em 1948, a lógica do método já era

conhecida há bastante tempo.

Por exemplo, existe um registro de um artigo

escrito por Lord Kelvin dezenas de anos antes

que já utilizava técnicas de Monte Carlo em

uma discussão das equações de Boltzmann.

1. INTRODUÇÃO

Monte Carlo é uma cidade de Mônaco,

a nordeste da capital do Principado

(Monaco-Ville).

Conhecida estância luxuosa, conhecida

pelo seu glamour, celebridades que

enxameiam as revistas cor de rosa, praias

e cassinos.

1. INTRODUÇÃO

É em Monte Carlo que se situa o Circuito

de Mônaco, onde ocorre o Grande

Prêmio de Mônaco de Fórmula Um.

É palco, ainda, de competições de boxe,

apresentações de moda e outros eventos

de grande repercussão cultural.

1. INTRODUÇÃO

1. INTRODUÇÃO

De uma maneira geral o MMC pode ser

sumarizado como:

O evento irá ocorrer com probabilidade p.

0 232-1

1.0 0.0 p

ocorre

Rp

não ocorre

R >p

2. Classes do MME

Existem três classes de algoritmos Monte

Carlo:

Erro-Unilateral;

Erro-Bilateral e

Erro-Não-Limitado.

2.1. MME de Erro-Unilateral

Seja P um problema e A um algoritmo aleatório, A é um algoritmo Monte Carlo de Erro-Unilateral que resolve P se, e somente se:

i) para toda configuração x que é solução de P:

ii) para toda configuração x que não é solução de P:

2

1 SIMxAprob

1 NÃOxAprob

2.1. MME de Erro-Unilateral

Ou seja, sempre que a resposta é NÃO, o

algoritmo garante a certeza da resposta.

Contudo, se a resposta for SIM, o

algoritmo não garante que a resposta

está correta.

2.2. MME de Erro-Bilateral

Um algoritmo aleatório A é um algoritmo

de Monte Carlo de Erro-Bilateral que

computa o problema F se existe um

número real ε, tal que para toda instância

x de F:

2

1xFxAprob

2.3. MME de Erro-Não-Limitado

Os algoritmos de Monte Carlo de Erro-

Não-Limitado são comumente chamados

de Algoritmos Monte Carlo.

Um algoritmo aleatório A é um algoritmo

de Monte Carlo se para qualquer

entrada x do problema F:

2

1 xFxAprob

3. Algoritmo de Metropolis

O algoritmo de Metropolis, também conhecido por Algoritmo de Metropolis-Hastings, foi apresentado inicialmente em 1953 num artigo de Nicholas Metropolis, Arianna Rosenbluth, Marshall Rosenbluth, Augusta Teller e Edward Teller.

Em 1970 este algoritmo foi generalizado por W. K. Hastings.


Este é provavelmente o método de Monte

Carlo mais utilizado na Física, e tem como

objetivo determinar valores esperados de

propriedades do sistema simulado, através

de uma média sobre uma amostra.

O algoritmo é concebido de modo a se

obter uma amostra que siga a distribuição de

Boltzmann.

3. Algoritmo de Metropolis Para se determinar a probabilidade de uma

dada configuração, seria necessário conhecer a chance de ocorrência de todas as outras configurações.

No caso de variáveis contínuas, seria necessário uma integração da densidade de probabilidade sobre todo o espaço de configurações, mas esse procedimento fica muito custoso quando se utiliza um número de variáveis da ordem de centenas.


A eficiência do algoritmo de Metropolis

está diretamente ligada ao fato de não

levar em conta a probabilidade das

configurações em si, mas sim a razão

entre elas, pois a razão entre as

probabilidades de duas dadas

configurações pode ser determinada

independentemente das outras.


Dadas duas configurações m e n

quaisquer, a razão entre a probabilidade

da configuração m, Pm, e a

probabilidade da configuração n, Pn,

pode ser escrita como:

Tk

UU

Tk

U

Tk

U

P

P

b

nm

b

n

b

m

n

m exp

exp

exp


A partir dessa igualdade, o algoritmo de Metropolis pode ser implementado através do seguinte conjunto de regras:

(a) Geração de uma configuração inicial aleatória, ou seja, com valores aleatórios para todos os graus de liberdade do sistema, respeitando as suas restrições. Vamos atribuir o índice m a essa configuração, que é aceita para a amostra.

3. Algoritmo de Metropolis (b) Geração de uma nova configuração-

tentativa de índice n, resultado de pequenas alterações nas coordenadas da configuração m.

(c) Se a energia da configuração n for menor que a da configuração m, inclui-se a configuração n na nossa amostra, e se atribui a ela o índice m a partir desse momento. Caso contrário, realizam-se os passos descritos nos subitems (c1) e (c2) abaixo:

3. Algoritmo de Metropolis (c1) Gera-se um número aleatório entre 0 e 1;

(c2) Se esse número aleatório for menor que Pn/Pm, aceita-se na amostra a configuração n, e se atribui a ela o índice m. Caso contrário, o índice m permanece designando a configuração original.

(d) Repete-se os passos (b) e (c) até que

algum critério de parada seja satisfeito. Cada uma dessas repetições é dita um passo (ou iteração) de Monte Carlo (MC).

Algoritmo Metropolis

4. Exemplo 1 – π

Neste exemplo um algoritmo de Monte Carlo simples é usado para calcular o valor de π usando uma seqüência de números aleatórios.

Considere o quadrado colocado no plano x-y com a extremidade inferior esquerda na origem como mostrado na figura abaixo.

4. Exemplo 1 – π

x

y

r

Área A

4. Exemplo 1 – π

A área do quadrado é r2, onde r é o

comprimento do seu lado.

Um quarto de um círculo é inscrito no

quadrado.

Seu raio é r e seu centro está na origem

do plano x-y.

4. Exemplo 1 – π

A área da quarta parte do círculo A é

dada por:

Pelo Método de Monte Carlo:

2

4

1rA

total

in

P

PlA 2

4. Exemplo 1 – π

Onde Pin é o número de pontos sorteados

aleatoriamente que ficaram dentro da

área em amarelo, e o Ptotal é o número

total de pontos gerados.

Considerando que o lado = raio = 1,

temos:

total

in

P

P4

4. Exemplo 1 – π

O algoritmo de cálculo consiste então

em gerar pontos P(x, y), sendo x e y [0,

1] e, de posse deste ponto verificar se

este se encontra (ou não) dentro da área

do círculo.

Lembrando que equação da

circunferência é:

1222 ryx

4. Exemplo 1 – π

4. Exemplo 1 – π

../Softwares/Integral via Monte Carlo/IntegralMonteCarlo.exe

4. Exemplo 2 –

Decaimento nuclear

Os núcleos radioativos possuem uma

probabilidade p de decair em cada

intervalo de tempo.

Desta forma o processo de decaimento

nuclear é um processo determinístico

dado por:

4. Exemplo 2 –

Decaimento nuclear

N1 = (1 – p) N0

N2 = (1 – p)2 N0

N3 = (1 – p)3 N0

N4 = (1 – p)4 N0

.

.

.

Ni = (1 – p)i N0

4. Exemplo 2 –

Decaimento nuclear

Se agora modelarmos o processo cmo

sendo um processo estocástico podemos

verificar individualmente quantos núcleos

decaem em cada intervalo de tempo,

da seguinte forma:

1.0 0.0 p decai não decai

p=0.03

p=0.06

4. Exemplo 3 –

Partículas em uma caixa

Uma caixa é separada em dois compartimentos.

Inicialmente todas as partículas estão no

compartimento da esquerda, quando uma

passagem é aberta.

4. Exemplo 3 –


A cada intervalo de tempo uma partícula

é escolhida para trocar de

compartimento.

Todas as partículas tem a mesma

probabilidade de serem selecionadas.

4. Exemplo 3 –


Sorteia-se uma das N partículas e troca a

partícula de compartimento.

Como as partículas são idênticas

(indistinguíveis), temos:

1.0 0.0

esquerda/N

esquerda -1 esquerda +1

4. Exemplo 3 –


Para t >> 1 teremos (em média) a mesma

quantidade de partículas nos dois

compartimentos.

4. Exemplo 4 –

Caminhada aleatória

Uma dada partícula escolherá

aleatoriamente para qual direção irá se

mover.

Tanto o tempo quanto o espaço são

discretos.

4. Exemplo 4 –


Caminhada aleatória em 1D

A partícula escolherá aleatoriamente para

qual direção irá se mover (direita ou

esquerda).

-5 -4 -3 -2 -1 0 1 2 3 4 5

4. Exemplo 4 –





esquerda).

-5 -4 -3 -2 -1 0 1 2 3 4 5

4. Exemplo 4 –





esquerda).

-5 -4 -3 -2 -1 0 1 2 3 4 5

4. Exemplo 4 –





esquerda).

-5 -4 -3 -2 -1 0 1 2 3 4 5

4. Exemplo 4 –





esquerda).

-5 -4 -3 -2 -1 0 1 2 3 4 5

4. Exemplo 4 –


Caminhada aleatória com tendência

Persistente – prefere se mover na mesma direção do passo anterior.

Probabilidade p de manter a direção.

p = ½ caminho aleatório padrão

p = 1 movimento uniforme

p =0 vai e volta

Anti-persistente – prefere se mover direção contrária ao passo anterior.

4. Exemplo 4 –




qual direção irá se mover (cima, baixo,

direita ou esquerda).

R

4. Exemplo 4 –



4. Exemplo 4 –



4. Exemplo 4 –



../Softwares/Caminhada aleatória/Caminhada aleatória 360/CaminhaAleatoria360.exe

4. Exemplo 4 –


Caminhada aleatória em 2D - Self-avoid

walker

Ele não passa por uma posição que já

tenha sido ocupada anteriormente.

Desta forma deve-se verificar quais são as

direções permitidas antes de dar o passo.

4. Exemplo 4 –



walker

A simulação termina quando não existe

mais opção de movimentação.

Esta técnica é amplamente utilizada para

representar crescimento de polímeros.

4. Exemplo 4 –



walker

Fim da simulação

4. Exemplo 4 –



walker

Existem caminhantes mais inteligentes, com

previsão de dois passos, para evitar

armadilhas.

Simulação ainda continua...

4. Exemplo 4 –



walker

Pode-se incluir campos que interferem na

direção ou na probabilidade de

crescimento do polímero, simulando

experimentos reais.

4. Exemplo 4 –



walker

../Softwares/Caminhada aleatória/Ver 2.3 - R250/CaminhadaAleatoria.exe

5. Integração pelo MMC

A integração deve ocorrer em um

intervalo limitado (cálculo de uma área).

As coordenadas (x,y) dos pontos são

escolhidas aleatoriamente, sorteando-se

dois números aleatórios independentes.


a b x

y

f(x)

N

nAdxxf

b

a


a b x

y

f(x)

g(x)

N

nAdxxgxf

b

a

)(


O gráfico abaixo apresenta o resultado

da integral de f(x)=1/(1+x2) no intervalo

[0,1] analiticamente e via MMC.


Como exemplo seja estimar a área

ocupada por cada espécie em um

coral.

1. Identificar a região de área A;

2. Escolher aleatoriamente as coordenadas

do ponto;

3. Identificar a espécie que ocupa o local e

fazer a contagem.

5.1. Integração pelo MMC

em n dimensões

No caso de n dimensões, devemos escolher n

coordenadas aleatórias e verificar se o ponto

está no interior da região de interesse, ou

não.

Tomemos como exemplo uma hiperesfera

(S3), sendo esta definida como o lugar

geométrico dos pontos (x,y,z,w) do R4 que

satisfazem a relação x2 + y2 + z2 + w2 = r2.

Introdução

Construção de uma teoria microscópica

da transição ferromagnética;

Modelo simples (mas não trivial);

Interações de curto alcance numa rede d-

dimensional.

6. Modelo de Ising

6. Modelo de Ising Hamiltoniano de spin:

Onde:

σi = ±1;

J é a constante de acoplamento entre os spins;

H é o campo magnético externo.

N

i

i

N

ji

ji HJ1,

H

Interpretações:

Spin “aponta” para cima ou para baixo;

Ocupação de um sítio por uma partícula

do tipo A ou B, como numa liga binária do

tipo AB (vizinhos iguais contribuem com

uma energia –J; vizinhos distintos, com uma

energia +J);

6. Modelo de Ising

Interpretações:

Número de ocupação – presença ou

ausência de uma molécula numa

determinada célula de um “gás de rede”.

6. Modelo de Ising

6. Modelo de Ising

Resolver o modelo de Ising significa escrever a função de partição canônica.

Várias técnicas foram desenvolvidas para tentar resolver o modelo de Ising em 2 ou 3 dimensões.

N

N

i

NHTZZ

Hexp,,

Podemos estudar fenômenos magnéticos

a partir do modelo de Ising, que

considera momentos magnéticos (spins)

localizados, distribuídos em uma rede.

No caso mais simples, consideramos que

cada spin interage apenas com os seus

primeiros vizinhos.

6. Modelo de Ising

Interações entre segundos vizinhos podem ser

incluídas no modelo, possibilitando o

aparecimento de interações competitivas,

dependendo de alguns parâmetros físicos.

O modelo de Ising pode ser estudado de

diferentes maneiras, entre elas a técnica de

Monte Carlo e a aproximação de Campo

médio.

6. Modelo de Ising

A técnica de Monte Carlo consiste em

obter uma seqüência de configurações

do sistema de uma maneira estocástica,

que depende de números aleatórios

gerados durante a simulação.

6. Modelo de Ising

6. Modelo de Ising

Para uma quantidade Q, magnetização

ou energia, por exemplo, o seu valor

médio é definido por:

M

i

cE

M

i

cE

i

i

i

e

ecQ

Q

1

1

Onde E(ci) a energia da n-ésima configuração.

Como é praticamente impossível que todas as configurações sejam somadas em um sistema com muitos graus de liberdade, então soma-se as M configurações mais importantes numa dada temperatura.

6. Modelo de Ising

6. Modelo de Ising

Suponhamos que essas configurações

mais importantes sejam escolhidas com

uma certa probabilidade.

Então passa a ser:

iM

i

cE

i

M

i

cE

i

cPe

cPecQ

Qi

i

/

/

1

1

No modelo de Ising com técnica de

Monte Carlo, começa-se com uma rede

de spins Si = +/- 1 dispostos

aleatoriamente.

Escolhe-se um spin e faz-se a mudança

Si = - Si

6. Modelo de Ising

Se a diferença de energia entra a

configuração após a mudança e antes

dela ΔE = E(-Si) – E(Si), for negativa

aceitamos essa mudança em Si.

Se ΔE > 0, aceitamos a nova

configuração com a probabilidade

exp(ΔE / kT).

6. Modelo de Ising

Na simulação isso é feito sorteando um

número aleatório r no intervalo de 0 < r <

1 e se exp(ΔE / kT) > r também aceitamos

a mudança.

Só rejeitamos se exp(ΔE / kT) < r.

6. Modelo de Ising

6. Modelo de Ising

É percorrida toda a rede fazendo essa

mudança, após terminado, pega-se a

configuração final de spins para calcular

a média termodinâmica de acordo com

a equação:

M

i

icQM

Q1

1

Documents

Modelamento e simulação de processos