1 / 65

Programação para as Ciências Experimentais 2007/8

Programação para as Ciências Experimentais 2007/8. Teórica 9. 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. Trabalho 1. Problema: Ler sequências (lefasta) Devolve matrizes com nomes e sequências.

lotta
Download Presentation

Programação para as Ciências Experimentais 2007/8

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Programação para as Ciências Experimentais2007/8 Teórica 9 Ludwig Krippahl, 2008

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

  3. Trabalho 1 • Problema: • Ler sequências (lefasta) • Devolve matrizes com nomes e sequências. • Ler a tabela com pks. • Encontrar grupos na sequência • Calcular o zero

  4. Trabalho 1 function t=letabela(fich) id=fopen(fich,"r"); ix=1; fgetl(id); while !feof(id) m=split(fgetl(id),"\t"); t(ix).cod=m(1,1); t(ix).coh=str2num(m(2,:)); t(ix).nh=str2num(m(3,:)); t(ix).cl=str2num(m(4,:)); t(ix).cd=str2num(m(5,:)); ix=ix+1; endwhile fclose(id);

  5. Trabalho 1 • Problema: • Ler sequências (lefasta) • Devolve matrizes com nomes e sequências. • Ler a tabela com pks (letabela). • Devolve vector de estruturas • Encontrar grupos na sequência • Calcular o zero

  6. Trabalho 1 • Problema: • Ler sequências (lefasta) • Devolve matrizes com nomes e sequências. • Ler a tabela com pks (letabela). • Devolve vector de estruturas • Encontrar grupos na sequência • Converter logo tudo... • Calcular o zero

  7. Trabalho 1 • seq2pks • Recebe sequência e tabela com os pks • Devolve matriz de 3 colunas e uma linha para cada grupo • 1ª coluna: o número de vezes que o grupo ocorre • 2ª coluna: o pka do grupo • 3ª coluna: a carga desprotonada do grupo • Desta forma só é preciso consultar a tabela uma vez para cada sequência.

  8. Trabalho 1 function pks=seq2pks(seq,tabela) ultimo=length(seq); pks=[]; for f=1:length(tabela) if strcmp(tabela(f).cod,seq(1)) pks=[pks;1,tabela(f).nh,0]; endif if strcmp(tabela(f).cod,seq(ultimo)) pks=[pks;1,tabela(f).coh,-1]; endif ...(cadeias laterais) endfor endfunction Percorre a tabela

  9. Trabalho 1 function pks=seq2pks(seq,tabela) ultimo=length(seq); pks=[]; for f=1:length(tabela) if strcmp(tabela(f).cod,seq(1)) pks=[pks;1,tabela(f).nh,0]; endif if strcmp(tabela(f).cod,seq(ultimo)) pks=[pks;1,tabela(f).coh,-1]; endif ...(cadeias laterais) endfor endfunction Se o aminoácido nesta posição da tabela é o primeiro da sequência, acrescenta NH3

  10. Trabalho 1 function pks=seq2pks(seq,tabela) ultimo=length(seq); pks=[]; for f=1:length(tabela) if strcmp(tabela(f).cod,seq(1)) pks=[pks;1,tabela(f).nh,0]; endif if strcmp(tabela(f).cod,seq(ultimo)) pks=[pks;1,tabela(f).coh,-1]; endif ...(cadeias laterais) endfor endfunction Se o aminoácido nesta posição da tabela é o último da sequência, acrescenta COOH

  11. Trabalho 1 function pks=seq2pks(seq,tabela) ultimo=length(seq); pks=[]; for f=1:length(tabela) if strcmp(tabela(f).cod,seq(1)) pks=[pks;1,tabela(f).nh,0]; endif if strcmp(tabela(f).cod,seq(ultimo)) pks=[pks;1,tabela(f).coh,-1]; endif ...(cadeias laterais) endfor endfunction Processa as cadeias laterais

  12. Trabalho 1 if !isempty(tabela(f).cl) v=findstr(seq,tabela(f).cod); if !isempty(v) pks=[pks;length(v), tabela(f).cl,tabela(f).cd]; endif endif Se tem cadeia lateral

  13. Trabalho 1 if !isempty(tabela(f).cl) v=findstr(seq,tabela(f).cod); if !isempty(v) pks=[pks;length(v), tabela(f).cl,tabela(f).cd]; endif endif Conta os aminoácidos destes na sequência

  14. Trabalho 1 if !isempty(tabela(f).cl) v=findstr(seq,tabela(f).cod); if !isempty(v) pks=[pks;length(v), tabela(f).cl,tabela(f).cd]; endif endif Se há acrescenta esse número de grupos à matriz dos pks

  15. Trabalho 1 • Problema: • Ler sequências (lefasta) • Devolve matrizes com nomes e sequências. • Ler a tabela com pks (letabela). • Devolve vector de estruturas • Encontrar grupos na sequência • Converter logo tudo... • Calcular o zero • Usa a função zeropol modificada

  16. Trabalho 1 function pH=pisoelec(pks,prec) x1=1; x2=13; y1=carga(pks,x1); y2=carga(pks,x2); while abs(x1-x2)>prec (bipartição...) endwhile pH=xm; endfunction Recebe a matriz com os dados dos grupos e a precisão.

  17. Trabalho 1 function pH=pisoelec(pks,prec) x1=1; x2=13; y1=carga(pks,x1); y2=carga(pks,x2); while abs(x1-x2)>prec (bipartição...) endwhile pH=xm; endfunction Intervalo inicial (x1 e x2 podiam vir como argumentos)

  18. Trabalho 1 function pH=pisoelec(pks,prec) x1=1; x2=13; y1=carga(pks,x1); y2=carga(pks,x2); while abs(x1-x2)>prec (bipartição...) endwhile pH=xm; endfunction A carga total é calculada em função do pH e da matriz com os pks e cargas.

  19. Trabalho 1 function pH=pisoelec(pks,prec) x1=1; x2=13; y1=carga(pks,x1); y2=carga(pks,x2); while abs(x1-x2)>prec (bipartição...) endwhile pH=xm; endfunction O mesmo do zeropol...

  20. Trabalho 1 • Calcular a carga de um grupo: function c=hendhas(pH,pKa,ch) lf=pH-pKa;%log 10 da proporção [A-]/[AH] f=10.^lf; % proporção [A-]/[AH]; A=f./(1+f); % fracção [A-] c=ch.*A+(ch+1).*(1-A); endfunction

  21. Trabalho 1 • Calcular a carga da sequência: function c=carga(pks,pH) c=0; for f=1:rows(pks) c=c+pks(f,1)*hendhas(pH,pks(f,2),pks(f,3)); endfor endfunction

  22. Trabalho 1 • Calcular a carga da sequência: • Mas como hendhas está preparada para vectores, podia ser: function c=carga(pks,pH) c=sum(pks(:,1).*hendhas(pH,pks(:,2),pks(:,3))); endfunction

  23. Trabalho 1 function vpis=calculapis(fich) [ns,seqs]=lefasta(fich); tabela=letabela("pKas.txt"); vpis=[]; for f=1:rows(seqs) pks=seq2pks(deblank(seqs(f,:)),tabela); vpis=[vpis;pisoelec(pks,0.1)]; endfor endfunction

  24. Trabalho 1 • calculapis • lefasta • letabela • seq2pks • pisoelec • carga • hendhas

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

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

  27. Integração numérica • Quanto menor dx mais preciso • dx*y

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

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

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

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

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

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

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

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

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

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

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

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

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

  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

  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 Próximo valor

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

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

  45. Integrar um sistema de equações diferenciais. function tcs=reacabc(cis,dt,tmax,k) tcs=[0,cis];%valores em t0 for t=0: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

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

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

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

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

  50. 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]

More Related