Programação para as
Ciências Experimentais
2007/8
Teórica 9
Ludwig Krippahl, 2008
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.
Ludwig Krippahl, 2008
2
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
Ludwig Krippahl, 2008
3
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);
Ludwig Krippahl, 2008
4
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
Ludwig Krippahl, 2008
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
• Converter logo tudo...
• Calcular o zero
Ludwig Krippahl, 2008
6
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.
Ludwig Krippahl, 2008
7
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
Ludwig Krippahl, 2008
Percorre a tabela
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
Ludwig Krippahl, 2008
Se o aminoácido
nesta posição da
tabela é o primeiro
da sequência,
acrescenta NH3
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
Ludwig Krippahl, 2008
Se o aminoácido
nesta posição da
tabela é o último da
sequência,
acrescenta COOH
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
Ludwig Krippahl, 2008
Processa as
cadeias laterais
11
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
Se tem cadeia
endif
lateral
Ludwig Krippahl, 2008
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
Conta os aminoácidos
endif
destes na sequência
Ludwig Krippahl, 2008
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
Se há acrescenta esse
endif
número de grupos à matriz
dos pks
Ludwig Krippahl, 2008
14
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
Ludwig Krippahl, 2008
15
Trabalho 1
function pH=pisoelec(pks,prec)
x1=1;
x2=13;
Recebe a matriz com os
y1=carga(pks,x1);
dados dos grupos e a
y2=carga(pks,x2);
precisão.
while abs(x1-x2)>prec
(bipartição...)
endwhile
pH=xm;
endfunction
Ludwig Krippahl, 2008
16
Trabalho 1
function pH=pisoelec(pks,prec)
x1=1;
x2=13;
Intervalo inicial (x1 e x2
y1=carga(pks,x1);
podiam vir como
y2=carga(pks,x2);
argumentos)
while abs(x1-x2)>prec
(bipartição...)
endwhile
pH=xm;
endfunction
Ludwig Krippahl, 2008
17
Trabalho 1
function pH=pisoelec(pks,prec)
x1=1;
x2=13;
A carga total é calculada
y1=carga(pks,x1);
em função do pH e da
y2=carga(pks,x2);
matriz com os pks e
while abs(x1-x2)>prec cargas.
(bipartição...)
endwhile
pH=xm;
endfunction
Ludwig Krippahl, 2008
18
Trabalho 1
function pH=pisoelec(pks,prec)
x1=1;
x2=13;
O mesmo do zeropol...
y1=carga(pks,x1);
y2=carga(pks,x2);
while abs(x1-x2)>prec
(bipartição...)
endwhile
pH=xm;
endfunction
Ludwig Krippahl, 2008
19
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
Ludwig Krippahl, 2008
20
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
Ludwig Krippahl, 2008
21
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
Ludwig Krippahl, 2008
22
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
Ludwig Krippahl, 2008
23
Trabalho 1

calculapis
• lefasta
• letabela
• seq2pks
• pisoelec
• carga
• hendhas
Ludwig Krippahl, 2008
24
Integração numérica.

y = exp(-x3)
Ludwig Krippahl, 2008
25
Integração numérica


Aproximar a função considerando cada
rectângulo
dx*y
Ludwig Krippahl, 2008
26
Integração numérica


Quanto menor dx mais preciso
dx*y
Ludwig Krippahl, 2008
27
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
Ludwig Krippahl, 2008
28
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
Ludwig Krippahl, 2008
29
Integração numérica

Aproximar melhor pela regra do trapézio
Ludwig Krippahl, 2008
30
Integração numérica

Àrea: base*(a+b)/2
a
b
base
Ludwig Krippahl, 2008
31
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
Ludwig Krippahl, 2008
y2
32
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
Ludwig Krippahl, 2008
33
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
Ludwig Krippahl, 2008
34
Integração numérica

Implementação mais genérica:
• Separar a função que calcula o y da função
que integra.
Ludwig Krippahl, 2008
35
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
Ludwig Krippahl, 2008
36
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
Ludwig Krippahl, 2008
37
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!
Ludwig Krippahl, 2008
38
E se não podemos traçar a
função?

Exemplo:

Não podemos calcular a área
geometricamente pelo gráfico da função.
• A+BC
• d[C]/dt = K [A] [B]
• d[A]/dt = -K [A] [B]
• d[B]/dt = -K [A] [B]
Ludwig Krippahl, 2008
39
E se não podemos traçar a
função?

Exemplo:

Mas podemos usar um método
semelhante: método de Euler
• A+BC
• d[C]/dt = K [A] [B]
• d[A]/dt = -K [A] [B]
• d[B]/dt = -K [A] [B]
Ludwig Krippahl, 2008
40
Integrar um sistema de
equações diferenciais.

Inicio, t0

Passo 1
• [A]0 [B]0 [C]0
• Usar valores em t0 para calcular derivada em
•
t0
Usar derivada para extrapolar t1:
• [A]1 = [A]0 + d[A]/dt * passo
Ludwig Krippahl, 2008
41
Integrar um sistema de
equações diferenciais.

Inicio, t0

Passo 1
• [A]0 [B]0 [C]0
• 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
Ludwig Krippahl, 2008
42
Integrar um sistema de
equações diferenciais.

Inicio, t0

Passo 1
• [A]0 [B]0 [C]0
• Usar valores em t0 para calcular derivada em
•
t0
Usar derivada para extrapolar t1:
• [A]1 = [A]0 + d[A]/dt * passo
Valor anterior
Ludwig Krippahl, 2008
43
Integrar um sistema de
equações diferenciais.

Inicio, t0

Passo 1
• [A]0 [B]0 [C]0
• 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.
Ludwig Krippahl, 2008
44
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
%AeB
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
Ludwig Krippahl, 2008
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
Ludwig Krippahl, 2008
46
Integrar um sistema de
equações diferenciais.
Ludwig Krippahl, 2008
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
Ludwig Krippahl, 2008
48
Generalizando: cinética com
método de Euler

Recebe

Caso geral
• Estequiometria,
• Concentrações iniciais
• Kd, Ki, dt e tmax.
• Veloc. = kd*reagentesesteq - ki*produtosesteq
• Alterar as concentrações
• d[A]/dt= velocidade*esteq(A)
• cs=cs+derivada*dt
Ludwig Krippahl, 2008
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]
Ludwig Krippahl, 2008
50
Cinética com método de Euler

Função
function
tconcs=cinetica(er,ep,cis,kd,ki,dt,tmax)
Ludwig Krippahl, 2008
51
Cinética com método de Euler

Integração:
for t=0:dt:tmax
v=prod(cis.^er)*kd-prod(cis.^ep)*ki;
deriv=v*ep-v*er;
Velocidade da reacção
cis=cis+deriv*dt;
tconcs=[tconcs;t,cis];
endfor
Ludwig Krippahl, 2008
52
Cinética com método de Euler

Integração:
for t=0:dt:tmax
v=prod(cis.^er)*kd-prod(cis.^ep)*ki;
deriv=v*ep-v*er;
Calcula derivadas e
cis=cis+deriv*dt;
tconcs=[tconcs;t,cis]; actualiza concentrações
endfor
Ludwig Krippahl, 2008
53
Cinética com método de Euler

Integração:
for t=0:dt:tmax
v=prod(cis.^er)*kd-prod(cis.^ep)*ki;
deriv=v*ep-v*er;
Guarda valores na
cis=cis+deriv*dt;
tconcs=[tconcs;t,cis]; matriz
endfor
Ludwig Krippahl, 2008
54
Cinética com método de Euler

Reacção do tipo
• A+BC
kd=1;
constante directa
ki=0.5;
constante inversa
cis=[1,2,0];
concentrações t0
er=[1,1,0]
esteq. reagentes
ep=[0,0,1]
esteq. produtos
pontos=cinetica(er,ep,cis,kd,ki,0.1,5);
Ludwig Krippahl, 2008
55
Cinética com método de Euler

Reacção do tipo
• A+BC
hold off
axis
eixo automático
plot(pontos(:,1),pontos(:,2:columns(pontos)))
Ludwig Krippahl, 2008
56
Cinética com método de Euler
Ludwig Krippahl, 2008
57
Sistema de reacções.

Mesma abordagem, mas a estequiometria são
matrizes (uma linha por reacção) e cada
constante cinética um vector (um valor por
reacção)
Ludwig Krippahl, 2008
58
Sistema de reacções.

Alterações:
• Na iteração é preciso calcular primeiro todas
•
as derivadas tendo em conta todas as
reacções (um ciclo pelas linhas das matrizes).
Só depois de ter todas as derivadas é que se
actualiza todas as concentrações.
Ludwig Krippahl, 2008
59
Sistema de reacções.

Função
function tconcs=cineticas(ers,eps,cis,kds,kis,dt,tmax)
(para fazer na aula)
Ludwig Krippahl, 2008
60
Sistema de reacções.

Exemplo:
• Oscilador químico (Lotka)
• Não se conhece nenhuma assim, mas esta é
simples de modelar.
A + X  2X
X + Y  2Y
YQ
Ludwig Krippahl, 2008
61
Sistema de reacções.
ers=[1 1 0 0;0 1 1 0;0 0 1 0];
A+X
eps=[0 2 0 0;0 0 2 0;0 0 0 1];
2X
X+Y
kds=[0.05,1,1];
2Y
kis=[0,0,0];
YQ
cis=[50,1,1,0];
t=cineticas(ers,eps,cis,kds,kis,0.01,50);
plot(t(:,1),t(:,2:columns(t)));
Ludwig Krippahl, 2008
62
Sistema de reacções.
Ludwig Krippahl, 2008
63
Sistema de reacções.
Ludwig Krippahl, 2008
64
Dúvidas
Ludwig Krippahl, 2008
65