Programação para as
Ciências Experimentais
2007/8
Teórica 7
Ludwig Krippahl, 2008
Na aula de hoje...

Métodos estocásticos (Monte Carlo)

Simulação

Trabalho 1 (dúvidas)
• Calculo de áreas
• Integração de funções
• Contar unidades formadoras de colónias
Ludwig Krippahl, 2008
2
Monte Carlo


Nome cunhado pelo matemático
Nicholas Constantine Metropolis(19151999)
Conjunto de métodos baseados em
números aleatórios.
Ludwig Krippahl, 2008
3
Calcular uma área
y2+x2<4
x>1
y<x-1
Ludwig Krippahl, 2008
4
Calcular uma área
y2+x2<4
x>1
y<x-1
Área?
Ludwig Krippahl, 2008
5
Algoritmo

Pontos ao acaso no quadrado -2, 2.
Área?
Ludwig Krippahl, 2008
6
Algoritmo

Contamos os pontos dentro e fora.
Área?
Ludwig Krippahl, 2008
7
Algoritmo



Contamos os pontos dentro e fora.
A fracção de pontos dentro da área será
a proporção entre a área a medir e a
área do quadrado.
Quanto mais pontos melhor.
Ludwig Krippahl, 2008
8
Implementação

Separar as tarefas:
• Dentro ou fora?
• Uma função que recebe x, y
e devolve true ou
false conforme x, y está dentro da área que
queremos.
• Contar os pontos
• Outra função que recebe o nome da função que
testa, o rectângulo que inclui a área a medir, e o
número de pontos, e devolve a área pretendida.
Ludwig Krippahl, 2008
9
Implementação

Função triangcirc testa se está dentro
do triângulo e circulo
function dentro=triangcirc(x,y)
dentro= ( x^2+y^2<4 && x>1 && y<x-1);
endfunction
Ludwig Krippahl, 2008
10
Implementação

Função areamc:
• devolve área
• recebe
• nome da função teste (string, para o feval)
• mínimos de x e y
• máximos de x e y (definem o rectângulo)
• número de pontos
Ludwig Krippahl, 2008
11
Ludwig Krippahl, 2008
12
Ludwig Krippahl, 2008
13
Ludwig Krippahl, 2008
14
Ludwig Krippahl, 2008
15
Ludwig Krippahl, 2008
16
Implementação

Função areamc:
function
a=areamc(testfn,minx,miny,maxx,maxy,pontos)
•
valor a devolver com a área
Ludwig Krippahl, 2008
17
Implementação

Função areamc:
function
a=areamc(testfn,minx,miny,maxx,maxy,pontos)
•
string com o nome da função que testa cada ponto
Ludwig Krippahl, 2008
18
Implementação

Função areamc:
function
a=areamc(testfn,minx,miny,maxx,maxy,pontos)
•
rectângulo que contêm a área a determinar
Ludwig Krippahl, 2008
19
Implementação

Função areamc:
function
a=areamc(testfn,minx,miny,maxx,maxy,pontos)
•
número de pontos a testar
Ludwig Krippahl, 2008
20
Implementação

Função areamc:
• área (variável a) começa a zero
• calcular a largura e altura do rectângulo
• ciclo para testar o número especificado de
•
pontos.
no final, a área é o número de pontos dentro
a dividir pelo total de pontos e multiplicar pela
área do rectângulo
Ludwig Krippahl, 2008
21
Implementação

Função areamc:
• ciclo para testar o número especificado de
pontos.
• criar coordenadas x,y aleatórias, de minx a maxx,
e miny a maxy respectivamente.
• se feval(testefn, x, y) for verdadeiro então
incrementar a variável a (para contar o número de
pontos dentro da área)
Ludwig Krippahl, 2008
22
Implementação

Para fazer os gráficos:
• Além da área devolver também duas matrizes
com as coordenadas x,y dos pontos dentro e
fora.
function
[a,dentros,foras]=areamc(testfn,minx,miny,
maxx,maxy,pontos)
Ludwig Krippahl, 2008
23
Exemplo de utilização
octave:13> areamc("triangcirc",-2,-2,2,2,1000)
ans = 1.7440
Nota: chamando assim ignora os outros valores
devolvidos, dentros e foras.
Ludwig Krippahl, 2008
24
Exemplo de utilização
pontos=1000;
[a,ds,fs]=areamc("triangcirc",-2,-2,2,2,pontos);
hold off
axis("equal")
eixos iguais
title([num2str(pontos)," pontos"]);
plot(ds(:,1),ds(:,2),"og;;");
hold on;
plot(fs(:,1),fs(:,2),"or;;");
“or;;” – circulo, red, ;; indica que não tem legenda
Ludwig Krippahl, 2008
25
Dicas

Mais pontos, mais rigor:
Ludwig Krippahl, 2008
26
Dicas


Mais pontos, mais rigor.
O rectângulo deve estar o mais próximo
possível da área que queremos.
• Em vez de (-2,-2) a (2,2), usar (1, -2) a (2,1)
Ludwig Krippahl, 2008
27
Dicas
• Em vez de (-2,-2) a (2,2)
usar (1, -2) a (2,1)
• Área=1.7
Ludwig Krippahl, 2008
28
Integrar função
Ludwig Krippahl, 2008
29
Integrar função


O integral de f(x)=exp(-x3) não tem
solução analítica.
Mas o integral é a área:
Ludwig Krippahl, 2008
30
Integrar função

Podemos usar a areamc, só precisamos
de uma função nova:
function dentro=expxcubo(x,y)
dentro=y<=exp(-x^3);
endfunction
Ludwig Krippahl, 2008
31
Integrar função
Basta usar:
areamc("expxcubo",0,0,2,1.2,5000);
ans= 0.89952

Nota: neste caso só é devolvido o primeiro
valor (a).
Ludwig Krippahl, 2008
32
Integrar função
Para fazer o gráfico:
pontos=5000;
[a,ds,fs]=areamc("expxcubo",0,0,2,1.2,pontos);
clearplot
axis("equal")
title([num2str(pontos)," pontos"]);
plot(ds(:,1),ds(:,2),"og;;");
hold on;
plot(fs(:,1),fs(:,2),"or;;");

Ludwig Krippahl, 2008
33
Ludwig Krippahl, 2008
34
Contar microorganismos no ar




Bomba aspira ar.
Orifícios sobre placa.
Contar colónias.
Estimar UFCs.
Ludwig Krippahl, 2008
35
Contar microorganismos no ar

Problema:
• Podem entrar vários esporos ou bactérias
pelo mesmo orifício, resultando numa só
colónia.
Ar
Ludwig Krippahl, 2008
36
Contar microorganismos no ar

Problema:
• Podem entrar vários esporos ou bactérias
•
pelo mesmo orifício, resultando numa só
colónia.
Sabendo o número de colónias na placa,
quantas UFCs no ar?
Ludwig Krippahl, 2008
37
Contar microorganismos no ar

Dividir em 2 fases
• Simular o processo para calcular quantas
•
colónias sabendo o número de UFCs.
Usar a simulação com diferentes valores de
UFCs até que obter o número de colónias
observado.
Ludwig Krippahl, 2008
38
Simulação



Temos N orifícios e X UFCs.
Cada UFC pode entrar por qualquer dos
N orifícios.
O número de colónias será o número de
orifícios diferentes por onde entraram
UFCs
Ludwig Krippahl, 2008
39
Algoritmo



Para cada UFC seleccionar um orifício
aleatoriamente e marcar esse orifício.
Contar o número de orifícios marcados.
Repetir um número grande de vezes
(50, 100, ...) e tirar o valor médio.
Ludwig Krippahl, 2008
40
Implementação
Função
function
cs=colonias(buracos,ufcs,tentativas)

• Devolve o número de colónias
Ludwig Krippahl, 2008
41
Implementação
Função
function
cs=colonias(buracos,ufcs,tentativas)

• Número de orifícios
Ludwig Krippahl, 2008
42
Implementação
Função
function
cs=colonias(buracos,ufcs,tentativas)

• Número de UFCs no ar
Ludwig Krippahl, 2008
43
Implementação
Função
function
cs=colonias(buracos,ufcs,tentativas)

• Número de vezes que repete a simulação
para calcular a média
Ludwig Krippahl, 2008
44
Implementação

Dois ciclos for:

A cada tentativa somar a um contador o
número de orifícios marcados.
• número de tentativas e, para cada tentativa
• número de ufcs.
Ludwig Krippahl, 2008
45
Implementação
Dois ciclos for:
for a=1:4
for b=1:3
[a,b]
endfor
endfor

Ludwig Krippahl, 2008
46
Implementação
Dois ciclos for:
for a=1:4
for b=1:3
[a,b]
endfor
endfor

Ludwig Krippahl, 2008
Repetido 4 vezes
47
Implementação
Dois ciclos for:
for a=1:4
for b=1:3
[a,b]
endfor
endfor

Ludwig Krippahl, 2008
Repetido 3 vezes
cada uma das 4
do ciclo de fora
(12 vezes no total)
48
Implementação


Os orifícios podem ser representados
como um vector de zeros, e marcados
com 1.
O total de orifícios marcados é o
somatório do vector.
Ludwig Krippahl, 2008
49
Implementação

Como seleccionar o orifício
aleatoriamente.
• É preciso arredondar:
• round(x)
• floor(x)
inteiro mais próximo de x
maior inteiro menor que x
ix=floor(rand*buracos)+1; (de 1 a buracos)
Nota: o rand nunca devolve 1. Ver help.
Ludwig Krippahl, 2008
50
Implementação
Exemplo: 10 orifícios e 3 UFCs,
[0 0 0 0 0 0 0 0 0 0] inicio, todos vazios
[0 0 1 0 0 0 0 0 0 0] 3
[0 0 1 0 0 0 1 0 0 0] 5
[0 0 1 0 0 0 1 0 0 0] 3

2 Colónias
Ludwig Krippahl, 2008
51
Problema real: saber UFCs


A função colonias dá-nos o número de
colónias formadas na placa sabendo as
UFCs no ar.
O problema real é o inverso: as colónias
na placa é o que se observa, e o que
queremos saber é as UFCs do ar.
Ludwig Krippahl, 2008
52
Algoritmo


O número de colónias será sempre igual
ou inferior ao número de UFCs no ar.
Começar por UFCs= colónias, e ir
incrementando os UFCs até obter da
função colonias o número certo de
colónias.
Ludwig Krippahl, 2008
53
Implementação

Função contaufcs
function
u=contaufcs(buracos,cs,tentativas)
• Recebe o número de orifícios, colónias
observadas, e o número de tentativas para
estimar as colónias para cada valor de u.
Ludwig Krippahl, 2008
54
Implementação

Função contaufcs
• u = cs
• Enquanto o número estimado de colónias for
•
inferior a cs, incrementar u e recalcular a
estimativa
Para estimar o número de colónias em função
de u usar a função colonias, com o número
de tentativas no argumento do contaufcs
Ludwig Krippahl, 2008
55
Ficha 7

Implementar:
• areamc
• colonias
• contaufcs
Ludwig Krippahl, 2008
56
Trabalho 1: Dúvidas
Ludwig Krippahl, 2008
57
Trabalho 1: Estrutura
Ponto Isoeléctrico
Ler sequências
Processar cada uma
Calcular o zero
Função carga
Ludwig Krippahl, 2008
Polinómios
Ler polinómios
Processar cada um
Calcular o zero
Função polinómio
58
Trabalho 1: Cálculo das cargas
Ponto Isoeléctrico
Polinómios
Função carga
Função polinómio
pH
X
Ler pKas.txt
Multiplicar coefs.
Consultar a tabela
Calcular cadeias laterais
Adicionar NH3 e COOH
Ludwig Krippahl, 2008
59
Trabalho 1: Cálculo da carga
pH = pKa + log (A/AH)
A + AH = 1
Resolvendo:
x = 10^(pH-pKa)
A= x/(1+x)
AH= 1-A
Carga média = A*CD + AH * (CD+1)
(CD é a carga da forma desprotonada)
Ludwig Krippahl, 2008
60
Dúvidas
Ludwig Krippahl, 2008
61