Programação para as
Ciências Experimentais
2008/9
Teórica 10
Ludwig Krippahl, 2009
Ficha 7

Não é para entregar

Aula de compensação?
• Mas é para fazer...
• Turnos P3-6
Ludwig Krippahl, 2009
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.
Ludwig Krippahl, 2009
3
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
Ludwig Krippahl, 2009
4
Objectivo

Acertar reacções químicas
?H2 + ?O2  ?H2O
Ludwig Krippahl, 2009
5
Objectivo

Acertar reacções químicas
?H2 + ?O2  ?H2O

Reacções simples (não redox, etc)
Nenhum termo com parêntesis

• Ca(NO3)2
Ludwig Krippahl, 2009
fica CaN2O6.
6
Objectivo

Acertar reacções químicas
?H2 + ?O2  ?H2O
2H2 + O2  2H2O
Ludwig Krippahl, 2009
7
Partir o problema

Ler o ficheiro
Estruturar cada reacção
Procurar coeficientes

Gravar resultados


• Testar se acertou
Ludwig Krippahl, 2009
8
Tirar os espaços
function t=semespacos(s)
t="";
for f=1:length(s)
if s(f)>" "
t=[t,s(f)];
endif
endfor
endfunction
Ludwig Krippahl, 2009
9
Tirar os espaços

Menos intuitivo, mas mais sucinto:
function t=semespacos(s)
t=s(s>" ");
endfunction
Ludwig Krippahl, 2009
10
Testar coeficientes

?H2 + ?O2  ?H2O
Vector de coeficientes:

Codificar em vectores tb

• [1 1 1]
• H: [2 0 -2]
• O: [0 2 -1]
(produtos a negativo)
Testar se sum(co.*el)==0
Ludwig Krippahl, 2009
11
Testar coeficientes
Estrutura com
.esteq
.el
A estequiometria para o teste
O elemento para podermos organizar a
estequiometria
Ludwig Krippahl, 2009
12
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
Ludwig Krippahl, 2009
13
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: ...
Ludwig Krippahl, 2009
14
Procurar coeficientes

testatodos.m
Ludwig Krippahl, 2009
15
Criar os coficientes

Lista de termos

Cada termo, decompor
Cada elemento, acrescentar a lista de
elementos

• .termo
• .produto
• .el
• .esteq
Ludwig Krippahl, 2009
16
Criar os coficientes

acrescentar.m

crialista.m
• Acrescentar cada elemento à lista
• Percorre lista de termos
• Decompõe e acrescenta
Ludwig Krippahl, 2009
17
Processar cada reacção



Partir em termos
Tirar os coeficientes
Gerar o vector com
• .termo
• .produto
Ludwig Krippahl, 2009
18
Processar cada reacção

tiranumero.m

separatermos.m
• Devolve o número (ou 1) e o resto do termo
• Recebe a reacção e devolve o vector com
• .termo
• .produto
• Devolve também a estequiometria inicial
Ludwig Krippahl, 2009
19
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
Ludwig Krippahl, 2009
20
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
Ludwig Krippahl, 2009
21
Trabalhos práticos


Começar com tempo
Aproveitar as aulas
• práticas e teóricas
Ludwig Krippahl, 2009
22
Integração numérica.

y = exp(-x3)
Ludwig Krippahl, 2009
23
Integração numérica


Aproximar a função considerando cada
rectângulo
dx*y
Ludwig Krippahl, 2009
24
Integração numérica


Quanto menor dx mais aproximado
dx*y
Ludwig Krippahl, 2009
25
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
Ludwig Krippahl, 2009
26
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, 2009
27
Integração numérica

Aproximar melhor pela regra do trapézio
Ludwig Krippahl, 2009
28
Integração numérica

Àrea: base*(a+b)/2
a
b
base
Ludwig Krippahl, 2009
29
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, 2009
y2
30
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, 2009
31
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, 2009
32
Integração numérica

Implementação mais genérica:
• Separar a função que calcula o y da função
que integra.
Ludwig Krippahl, 2009
33
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, 2009
34
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, 2009
35
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, 2009
36
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, 2009
37
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, 2009
38
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, 2009
39
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, 2009
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
Valor anterior
Ludwig Krippahl, 2009
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
Derivada vezes passo.
Ludwig Krippahl, 2009
42
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
%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, 2009
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
%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, 2009
44
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, 2009
45
Integrar um sistema de
equações diferenciais.
Ludwig Krippahl, 2009
46
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, 2009
47
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, 2009
48
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, 2009
49
Generalizando: cinética com
método de Euler

Exemplo:

A1*B1*kd
C2*ki

• A + B  2C
• er=[1, 1, 0]
• ep=[0, 0, 2]
Ludwig Krippahl, 2009
v directa
v inversa
50
Cinética com método de Euler

Função
function
tconcs=cinetica(er,ep,cis,kd,ki,dt,tmax)
Ludwig Krippahl, 2009
51
Cinética com método de Euler
Integração:
tconcs=[0,cis];
for t=dt: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]; Produto concentrações
elevadas à estequiometria
endfor

Ludwig Krippahl, 2009
52
Cinética com método de Euler
Integração:
tconcs=[0,cis];
for t=dt: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, 2009
53
Cinética com método de Euler
Integração:
tconcs=[0,cis];
for t=dt: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, 2009
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, 2009
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, 2009
56
Cinética com método de Euler
Ludwig Krippahl, 2009
57
Sistema de reacções.

Mesma abordagem, mas a estequiometria é
uma matriz (uma linha por reacção) e cada
constante cinética um vector (um valor por
reacção)
Ludwig Krippahl, 2009
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, 2009
59
Sistema de reacções.

Função
function tconcs=cineticas(ers,eps,cis,kds,kis,dt,tmax)
(para fazer na aula)
Ludwig Krippahl, 2009
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, 2009
61
Sistema de reacções.
ers=[1 1 0 0;0 1 1 0;0 0 1 0];
A + X  2X
eps=[0 2 0 0;0 0 2 0;0 0 0 1];
X + Y  2Y
YQ
kds=[0.05,1,1];
kis=[0,0,0];
cis=[50,1,1,0];
t=cineticas(ers,eps,cis,kds,kis,0.01,50);
plot(t(:,1),t(:,2:columns(t)));
Ludwig Krippahl, 2009
62
Sistema de reacções.
Ludwig Krippahl, 2009
63
Sistema de reacções.
Ludwig Krippahl, 2009
64
Resumo

Integração numérica
•
•
•
Se podemos calcular o y
•
•
Calcular pontos, y, multiplicar pelo dx
Melhor ainda, trapézio
Equações diferenciais
•
Método de Euler, calcular derivada em t, usar
valores em t e extrapolar para t+dt.
Cinética
•
Ludwig Krippahl, 2009
Estequiometria e constantes.
65
Dúvidas
Ludwig Krippahl, 2009
66