Upload
internet
View
106
Download
2
Embed Size (px)
Citation preview
Ludwig Krippahl, 2009
Programação para as Ciências Experimentais
2008/9
Teórica 10
Ludwig Krippahl, 2009 2
Ficha 7
Não é para entregar• Mas é para fazer...
Aula de compensação?• Turnos P3-6
Ludwig Krippahl, 2009 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.
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 5
Objectivo
Acertar reacções químicas
?H2 + ?O2 ?H2O
Ludwig Krippahl, 2009 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.
Ludwig Krippahl, 2009 7
Objectivo
Acertar reacções químicas
?H2 + ?O2 ?H2O
2H2 + O2 2H2O
Ludwig Krippahl, 2009 8
Partir o problema
Ler o ficheiro Estruturar cada reacção Procurar coeficientes
• Testar se acertou
Gravar resultados
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 10
Tirar os espaços
Menos intuitivo, mas mais sucinto:
function t=semespacos(s)
t=s(s>" ");
endfunction
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 12
Testar coeficientes
Estrutura com
.esteq
.el
A estequiometria para o teste
O elemento para podermos organizar a estequiometria
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 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: ...
Ludwig Krippahl, 2009 15
Procurar coeficientes
testatodos.m
Ludwig Krippahl, 2009 16
Criar os coficientes
Lista de termos• .termo
• .produto Cada termo, decompor Cada elemento, acrescentar a lista de
elementos• .el
• .esteq
Ludwig Krippahl, 2009 17
Criar os coficientes
acrescentar.m• Acrescentar cada elemento à lista
crialista.m• Percorre lista de termos
• Decompõe e acrescenta
Ludwig Krippahl, 2009 18
Processar cada reacção
Partir em termos Tirar os coeficientes Gerar o vector com
• .termo
• .produto
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 22
Trabalhos práticos
Começar com tempo Aproveitar as aulas
• práticas e teóricas
Ludwig Krippahl, 2009 23
Integração numérica.
y = exp(-x3)
Ludwig Krippahl, 2009 24
Integração numérica
Aproximar a função considerando cada rectângulo
dx*y
Ludwig Krippahl, 2009 25
Integração numérica
Quanto menor dx mais aproximado dx*y
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 28
Integração numérica
Aproximar melhor pela regra do trapézio
Ludwig Krippahl, 2009 29
Integração numérica
Àrea: base*(a+b)/2
a
b
base
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 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;
endforendfunction
Ludwig Krippahl, 2009 32
Integração numérica
octave:180> intexcubot(0.2,0,2)ans = 0.89293octave:181> intexpxcubo(0.2,0,2)ans = 0.99296octave:182> intexcubot(0.002,0,2)ans = 0.89295octave:183> intexpxcubo(0.002,0,2)ans = 0.89395
Ludwig Krippahl, 2009 33
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 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
Ludwig Krippahl, 2009 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;
endforendfunction
Ludwig Krippahl, 2009 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!
Ludwig Krippahl, 2009 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.
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 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
Ludwig Krippahl, 2009 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.
Ludwig Krippahl, 2009 43
Integrar um sistema de equações diferenciais.function tcs=reacabc(cis,dt,tmax,k)tcs=[0,cis]; %valores em t0for t=dt:dt:tmax
abk=cis(1)*cis(2)*k*dt; calcular a derivada% A e Bcis(1)=cis(1)-abk;cis(2)=cis(2)-abk;% Ccis(3)=cis(3)+abk; actualizar concentrações% guarda o novo pontotcs=[tcs;t,cis];
endforendfunction
Ludwig Krippahl, 2009 44
Integrar um sistema de equações diferenciais.function tcs=reacabc(cis,dt,tmax,k)tcs=[0,cis]; %valores em t0for t=dt:dt:tmax
abk=cis(1)*cis(2)*k*dt; calcular a derivada% A e Bcis(1)=cis(1)-abk;cis(2)=cis(2)-abk;% Ccis(3)=cis(3)+abk; actualizar concentrações% guarda o novo pontotcs=[tcs;t,cis];
endforendfunction
Ludwig Krippahl, 2009 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, 2009 46
Integrar um sistema de equações diferenciais.
Ludwig Krippahl, 2009 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, 2009 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
Ludwig Krippahl, 2009 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, 2009 50
Generalizando: cinética com método de Euler
Exemplo:• A + B 2C
• er=[1, 1, 0]
• ep=[0, 0, 2]
A1*B1*kdv directa
C2*ki v inversa
Ludwig Krippahl, 2009 51
Cinética com método de Euler
Função
function tconcs=cinetica(er,ep,cis,kd,ki,dt,tmax)
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; cis=cis+deriv*dt; tconcs=[tconcs;t,cis];endfor
Velocidade da reacção
Produto concentrações elevadas à estequiometria
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; cis=cis+deriv*dt; tconcs=[tconcs;t,cis];endfor
Calcula derivadas e actualiza concentrações
Ludwig Krippahl, 2009 54
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; cis=cis+deriv*dt; tconcs=[tconcs;t,cis];endfor
Guarda valores na matriz
Ludwig Krippahl, 2009 55
Cinética com método de Euler
Reacção do tipo• A + B C
kd=1; constante directaki=0.5; constante inversa
cis=[1,2,0]; concentrações t0
er=[1,1,0] esteq. reagentesep=[0,0,1] esteq. produtospontos=cinetica(er,ep,cis,kd,ki,0.1,5);
Ludwig Krippahl, 2009 56
Cinética com método de Euler
Reacção do tipo• A + B C
hold off
axiseixo automático
plot(pontos(:,1),pontos(:,2:columns(pontos)))
Ludwig Krippahl, 2009 57
Cinética com método de Euler
Ludwig Krippahl, 2009 58
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 59
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 60
Sistema de reacções.
Funçãofunction tconcs=cineticas(ers,eps,cis,kds,kis,dt,tmax)
(para fazer na aula)
Ludwig Krippahl, 2009 61
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 62
Sistema de reacções.
ers=[1 1 0 0;0 1 1 0;0 0 1 0];
eps=[0 2 0 0;0 0 2 0;0 0 0 1];
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)));
A + X 2XX + Y 2YY Q
Ludwig Krippahl, 2009 63
Sistema de reacções.
Ludwig Krippahl, 2009 64
Sistema de reacções.
Ludwig Krippahl, 2009 65
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• Estequiometria e constantes.
Ludwig Krippahl, 2009 66
Dúvidas