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

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

2. Ficha 7 • Não é para entregar • Mas é para fazer... • Aula de compensação? • Turnos P3-6

3. Na aula de hoje... • Resolução do trabalho 1 • Integração de funções de uma variável. • Integração de equações diferenciais. • Osciladores químicos.

4. Trabalho 1 • Discussões do trabalho • Não entregaram intermédia, ou muito incompleta • Não tivemos oportunidade de acompanhar • Muito diferente do que foi dado • ... • Lista a afixar, depois contactem o docente das práticas. • Se possível, será na aula prática

5. Objectivo • Acertar reacções químicas ?H2 + ?O2 ?H2O

6. Objectivo • Acertar reacções químicas ?H2 + ?O2 ?H2O • Reacções simples (não redox, etc) • Nenhum termo com parêntesis • Ca(NO3)2 fica CaN2O6.

7. Objectivo • Acertar reacções químicas ?H2 + ?O2 ?H2O 2H2 + O2 2H2O

8. Partir o problema • Ler o ficheiro • Estruturar cada reacção • Procurar coeficientes • Testar se acertou • Gravar resultados

9. Tirar os espaços function t=semespacos(s) t=""; for f=1:length(s) if s(f)>" " t=[t,s(f)]; endif endfor endfunction

10. Tirar os espaços • Menos intuitivo, mas mais sucinto: function t=semespacos(s) t=s(s>""); endfunction

11. Testar coeficientes ?H2 + ?O2 ?H2O • Vector de coeficientes: • [1 1 1] • Codificar em vectores tb • H: [2 0 -2] • O: [0 2 -1] (produtos a negativo) • Testar se sum(co.*el)==0

12. Testar coeficientes Estrutura com .esteq .el A estequiometria para o teste O elemento para podermos organizar a estequiometria

13. Testar coeficientes function r=testa(esteq,lista) r=true; for f=1:length(lista) if sum(esteq.*lista(f).esteq)!=0 r=false; break endif endfor endfunction

14. Procurar coeficientes • 0: [1 1 1] • 1: [2 1 1] [1 2 1] [1 1 2] • 2: [3 1 1] [2 2 1] [2 1 2] [2 2 1][1 3 1] [1 2 2] [2 1 2] [1 2 2][1 1 3] • 3: ...

15. Procurar coeficientes • testatodos.m

16. Criar os coficientes • Lista de termos • .termo • .produto • Cada termo, decompor • Cada elemento, acrescentar a lista de elementos • .el • .esteq

17. Criar os coficientes • acrescentar.m • Acrescentar cada elemento à lista • crialista.m • Percorre lista de termos • Decompõe e acrescenta

18. Processar cada reacção • Partir em termos • Tirar os coeficientes • Gerar o vector com • .termo • .produto

19. Processar cada reacção • tiranumero.m • Devolve o número (ou 1) e o resto do termo • separatermos.m • Recebe a reacção e devolve o vector com • .termo • .produto • Devolve também a estequiometria inicial

20. Processar cada reacção • Com uma reacção • Separar os termos • [ltermos,esteq]=separatermos(reac); • Criar o vector com as estequiometrias • lelems=crialista(ltermos); • Testar o inicial e depois todos • [encontra,esteq]=testatodos ... • acerta.m

21. Juntar tudo • Percorre ficheiro de entrada • Cada reacção • acerta(r,maximo); • Se vazio, não conseguiu, • Caso contrário, escreve reacção acertada • acertada.m • acertatudo.m

22. Trabalhos práticos • Começar com tempo • Aproveitar as aulas • práticas e teóricas

23. Integração numérica. • y = exp(-x3)

24. Integração numérica • Aproximar a função considerando cada rectângulo • dx*y

25. Integração numérica • Quanto menor dx mais aproximado • dx*y

26. Integração numérica function int=intexpxcubo(dx,x0,x1); int=0; for x=x0:dx:x1-dx int=int+dx*exp(-x^3); endfor endfunction

27. Integração numérica octave:114> intexpxcubo(0.2,0,2) ans = 0.99296 octave:115> intexpxcubo(0.02,0,2) ans = 0.90296 octave:116> intexpxcubo(0.002,0,2) ans = 0.89395 octave:117> intexpxcubo(0.0002,0,2) ans = 0.89305

28. Integração numérica • Aproximar melhor pela regra do trapézio

29. Integração numérica • Àrea: base*(a+b)/2 a b base

30. Integração numérica • Implementação: • Método ingénuo: calcular os dois y em cada iteração. • Método mais eficiente: calcular o y2 e guardá-lo no y1 para a próxima y1 y2

31. Integração numérica function int=intexcubot(dx,x0,x1) int=0; y1=exp(-x0^3); for x=x0+dx:dx:x1 y2=exp(-x^3); int=int+dx*(y1+y2)/2; y1=y2; endfor endfunction

32. Integração numérica octave:180> intexcubot(0.2,0,2) ans = 0.89293 octave:181> intexpxcubo(0.2,0,2) ans = 0.99296 octave:182> intexcubot(0.002,0,2) ans = 0.89295 octave:183> intexpxcubo(0.002,0,2) ans = 0.89395

33. Integração numérica • Implementação mais genérica: • Separar a função que calcula o y da função que integra.

34. Integração numérica • Implementação mais genérica: • Cálculo y em função de x • function y = nome(x) function y=expxcubo(x) y=exp(-x.^3); endfunction

35. Integração numérica • Implementação mais genérica: • Integração, chamando a função com feval function int=trapezio(funcao,dx,x0,x1) int=0; y1=feval(funcao,x0); for x=x0+dx:dx:x1 y2=feval(funcao,x); int=int+dx*(y1+y2)/2; y1=y2; endfor endfunction

36. Integração numérica • Implementação mais genérica: octave:185> trapezio("expxcubo",0.2,0,2) ans = 0.89293 Nota: “expxcubo” em vez de expxcubo!

37. E se não podemos traçar a função? • Exemplo: • A + B  C • d[C]/dt = K [A] [B] • d[A]/dt = -K [A] [B] • d[B]/dt = -K [A] [B] • Não podemos calcular a área geometricamente pelo gráfico da função.

38. E se não podemos traçar a função? • Exemplo: • A + B  C • d[C]/dt = K [A] [B] • d[A]/dt = -K [A] [B] • d[B]/dt = -K [A] [B] • Mas podemos usar um método semelhante: método de Euler

39. Integrar um sistema de equações diferenciais. • Inicio, t0 • [A]0 [B]0 [C]0 • Passo 1 • Usar valores em t0 para calcular derivada em t0 • Usar derivada para extrapolar t1: • [A]1 = [A]0 + d[A]/dt * passo

40. Integrar um sistema de equações diferenciais. • Inicio, t0 • [A]0 [B]0 [C]0 • Passo 1 • Usar valores em t0 para calcular derivada em t0 • Usar derivada para extrapolar t1: • [A]1 = [A]0 + d[A]/dt * passo Próximo valor

41. Integrar um sistema de equações diferenciais. • Inicio, t0 • [A]0 [B]0 [C]0 • Passo 1 • Usar valores em t0 para calcular derivada em t0 • Usar derivada para extrapolar t1: • [A]1 = [A]0 + d[A]/dt * passo Valor anterior

42. Integrar um sistema de equações diferenciais. • Inicio, t0 • [A]0 [B]0 [C]0 • Passo 1 • Usar valores em t0 para calcular derivada em t0 • Usar derivada para extrapolar t1: • [A]1 = [A]0 + d[A]/dt * passo Derivada vezes passo.

43. Integrar um sistema de equações diferenciais. function tcs=reacabc(cis,dt,tmax,k) tcs=[0,cis];%valores em t0 for t=dt:dt:tmax abk=cis(1)*cis(2)*k*dt; calcular a derivada % A e B cis(1)=cis(1)-abk; cis(2)=cis(2)-abk; % C cis(3)=cis(3)+abk; actualizar concentrações % guarda o novo ponto tcs=[tcs;t,cis]; endfor endfunction

44. Integrar um sistema de equações diferenciais. function tcs=reacabc(cis,dt,tmax,k) tcs=[0,cis];%valores em t0 for t=dt:dt:tmax abk=cis(1)*cis(2)*k*dt; calcular a derivada % A e B cis(1)=cis(1)-abk; cis(2)=cis(2)-abk; % C cis(3)=cis(3)+abk; actualizar concentrações % guarda o novo ponto tcs=[tcs;t,cis]; endfor endfunction

45. Integrar um sistema de equações diferenciais. cis=[1,1.5,0] k=1; pontos=reacabc(cis,0.1,10,k); hold off axis plot(pontos(:,1),pontos(:,2:columns(pontos))) Plot: primeira coluna é o tempo, as restantes colunas são as concentrações

46. Integrar um sistema de equações diferenciais.

47. Generalizando: cinética com método de Euler • Exemplo • A + B  C • Reacção reversível: kd e ki • Velocidade = [C]ki - [A][B]kd

48. Generalizando: cinética com método de Euler • Recebe • Estequiometria, • Concentrações iniciais • Kd, Ki, dt e tmax. • Caso geral • Veloc. = kd*reagentesesteq -ki*produtosesteq • Alterar as concentrações • d[A]/dt= velocidade*esteq(A) • cs=cs+derivada*dt

49. Generalizando: cinética com método de Euler • Devolve • Matriz com tempo (1ª coluna) e concentração por iteração (para fazer o gráfico). • Estequiometria, • 2 vectores, reagentes e produtos • Exemplo: • A + B  2C • er=[1, 1, 0] • ep=[0, 0, 2]

50. Generalizando: cinética com método de Euler • Exemplo: • A + B  2C • er=[1, 1, 0] • ep=[0, 0, 2] • A1*B1*kd v directa • C2*ki v inversa