Programação para as Ciências Experimentais 2008/9

1. Programação para as Ciências Experimentais2008/9 Teórica 7 Ludwig Krippahl, 2009

2. Na aula de hoje... • Métodos estocásticos (Monte Carlo) • Calculo de áreas • Integração de funções • Simulação • Contar unidades formadoras de colónias • Trabalho 1 (dúvidas)

3. Monte Carlo • Nome cunhado pelo matemático Nicholas Constantine Metropolis • (1915-1999) • Conjunto de métodos baseados em números aleatórios.

4. Calcular uma área x>1 y2+x2<4 y<x-1

5. Calcular uma área x>1 y2+x2<4 y<x-1 Área?

6. Algoritmo • Pontos ao acaso no quadrado -2, 2. Área?

7. Algoritmo • Contamos os pontos dentro e fora. Área?

8. 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.

9. 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.

10. 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

11. 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

17. Implementação • Função areamc: functiona=areamc(testfn,minx,miny,maxx,maxy,pontos) • valor a devolver com a área

18. 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

19. Implementação • Função areamc: function a=areamc(testfn,minx,miny,maxx,maxy,pontos) • rectângulo que contêm a área a determinar

20. Implementação • Função areamc: function a=areamc(testfn,minx,miny,maxx,maxy,pontos) • número de pontos a testar

21. 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

22. 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)

23. 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)

24. 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.

25. 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

26. Dicas • Mais pontos, mais rigor:

27. 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)

28. Dicas • Em vez de (-2,-2) a (2,2) usar (1, -2) a (2,1) • Área=1.7

29. Integrar função

30. Integrar função • O integral de f(x)=exp(-x3) não tem solução analítica. • Mas o integral é a área:

31. 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

32. 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).

33. 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;;");

35. Contar microorganismos no ar • Bomba aspira ar. • Orifícios sobre placa. • Contar colónias. • Estimar UFCs.

36. Contar microorganismos no ar • Problema: • Podem entrar vários esporos ou bactérias pelo mesmo orifício, resultando numa só colónia. Ar

37. 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?

38. 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.

39. 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

40. 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.

41. Implementação • Função functioncs=colonias(buracos,ufcs,tentativas) • Devolve o número de colónias

42. Implementação • Função function cs=colonias(buracos,ufcs,tentativas) • Número de orifícios

43. Implementação • Função function cs=colonias(buracos,ufcs,tentativas) • Número de UFCs no ar

44. Implementação • Função function cs=colonias(buracos,ufcs,tentativas) • Número de vezes que repete a simulação para calcular a média

45. Implementação • Dois ciclos for: • número de tentativas e, para cada tentativa • número de ufcs. • A cada tentativa somar a um contador o número de orifícios marcados.

46. Implementação • Dois ciclos for: for a=1:4 for b=1:3 disp([a,b]) endfor endfor

47. Implementação • Dois ciclos for: for a=1:4 for b=1:3 [a,b] endfor endfor Repetido 4 vezes

48. Implementação • Dois ciclos for: for a=1:4 for b=1:3 [a,b] endfor endfor Repetido 3 vezes cada uma das 4 do ciclo de fora (12 vezes no total)

49. 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.

50. Implementação • Como seleccionar o orifício aleatoriamente. • É preciso arredondar: • round(x) inteiro mais próximo de x • floor(x) maior inteiro menor que x ix=floor(rand*buracos)+1; (de 1 a buracos) Nota: o rand nunca devolve 1. Ver help.