Ludwig Krippahl, 2007 Programação para as Ciências Experimentais 2006/7 Teórica 7

Ludwig Krippahl, 2007

Programação para as Ciências Experimentais

2006/7

Teórica 7

Ludwig Krippahl, 2007 2

Aviso: P1 e P5

Aula prática extra na sexta, dia 4-5, das 8:00 às 11:00.

Ludwig Krippahl, 2007 3

Na aula de hoje...

Integração de funções de uma variável. Integração de equações diferenciais. Interpolação linear Ajuste de modelos a dados

experimentais

Ludwig Krippahl, 2007 4

Integração numérica.

y = exp(-x3)

Ludwig Krippahl, 2007 5

Integração numérica

Aproximar a função considerando cada rectângulo

dx*y

Ludwig Krippahl, 2007 6

Integração numérica

Quanto menor dx mais preciso dx*y

Ludwig Krippahl, 2007 7

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, 2007 8

Integração numérica

octave:113> intexpxcubo(0.2,0,2)ans = 0.99296octave:114> intexpxcubo(0.2,0,2)ans = 0.99296octave:115> intexpxcubo(0.02,0,2)ans = 0.90296octave:116> intexpxcubo(0.002,0,2)ans = 0.89395octave:117> intexpxcubo(0.0002,0,2)ans = 0.89305

Ludwig Krippahl, 2007 9

Integração numérica

Aproximar melhor pela regra do trapézio

Ludwig Krippahl, 2007 10

Integração numérica

Àrea: base*(a+b)/2

a

b

base

Ludwig Krippahl, 2007 11

Integração numérica

Implementação:• Método ingénuo: calcular os dois y em cada

iteração.

• Método mais inteligente: calcular o y2 e guardá-lo no y1 para a próxima

y1 y2

Ludwig Krippahl, 2007 12

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, 2007 13

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, 2007 14

Integração numérica

Implementação mais genérica:• Separar a função que calcula o y da função

que integra.

Ludwig Krippahl, 2007 15

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, 2007 16

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, 2007 17

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, 2007 18

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, 2007 19

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, 2007 20

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, 2007 21

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, 2007 22

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, 2007 23

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, 2007 24

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

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, 2007 26

Integrar um sistema de equações diferenciais.

Ludwig Krippahl, 2007 27

Generalizando: cinética com método de Euler

Separar a função que avalia a derivada. Estequiometria Exemplo

• A + B C

• Reacção reversível: kd e ki

• Derivada = [C]ki - [A][B]kd

Ludwig Krippahl, 2007 28

Generalizando: cinética com método de Euler

Separar a função que avalia a derivada. Estequiometria Caso geral

• Derivada = kd*reagentesesteq - ki*produtosesteq

• Alterar as concentrações• cs=cs+derivada*esteq

• Nota: reagentes têm estequiometria negativa

Ludwig Krippahl, 2007 29

Cinética com método de Euler

function tconcs=cinetica(esteq,cis,kd,ki,dt,tmax)rs=find(esteq<0);ps=find(esteq>0);tconcs=[0,cis];for t=0:dt:tmax

dps=prod(cis(ps).^esteq(ps))*ki;drs=prod(cis(rs).^-esteq(rs))*kd;deriv=(drs-dps)*dt;cis=cis+deriv*esteq;

tconcs=[tconcs;t,cis];endforendfunction

Ludwig Krippahl, 2007 30

Cinética com método de Euler

function tconcs=cinetica(esteq,cis,kd,ki,dt,tmax)rs=find(esteq<0);ps=find(esteq>0);tconcs=[0,cis];for t=0:dt:tmax

dps=prod(cis(ps).^esteq(ps))*ki;drs=prod(cis(rs).^-esteq(rs))*kd;deriv=(drs-dps)*dt;cis=cis+deriv*esteq;

tconcs=[tconcs;t,cis];endforendfunction

Encontrar os índices dos reagentes e produtos na estequiometria

Ludwig Krippahl, 2007 31

Cinética com método de Euler

esteq>0 Devolve vector 0s e 1soctave:110> vector=[1,2,-3,0,-4]vector = 1 2 -3 0 -4octave:111> vector>0ans = 1 1 0 0 0

Ludwig Krippahl, 2007 32

Cinética com método de Euler

find(esteq>0)Índices dos não 0

octave:113> find([0,1,0,-2,0,0.5])

ans =

2 4 6

Ludwig Krippahl, 2007 33

Cinética com método de Euler

function tconcs=cinetica(esteq,cis,kd,ki,dt,tmax)rs=find(esteq<0);ps=find(esteq>0);tconcs=[0,cis];for t=0:dt:tmax

dps=prod(cis(ps).^esteq(ps))*ki;drs=prod(cis(rs).^-esteq(rs))*kd;deriv=(drs-dps)*dt;cis=cis+deriv*esteq;

tconcs=[tconcs;t,cis];endforendfunction

Encontrar os índices dos reagentes e produtos na estequiometria

Ludwig Krippahl, 2007 34

Cinética com método de Euler

function tconcs=cinetica(esteq,cis,kd,ki,dt,tmax)rs=find(esteq<0);ps=find(esteq>0);tconcs=[0,cis];for t=0:dt:tmax

dps=prod(cis(ps).^esteq(ps))*ki;drs=prod(cis(rs).^-esteq(rs))*kd;deriv=(drs-dps)*dt;cis=cis+deriv*esteq;

tconcs=[tconcs;t,cis];endforendfunction

Guardar valores para t0

Ludwig Krippahl, 2007 35

Cinética com método de Euler

function tconcs=cinetica(esteq,cis,kd,ki,dt,tmax)rs=find(esteq<0);ps=find(esteq>0);tconcs=[0,cis];for t=0:dt:tmax

dps=prod(cis(ps).^esteq(ps))*ki;drs=prod(cis(rs).^-esteq(rs))*kd;deriv=(drs-dps)*dt;cis=cis+deriv*esteq;

tconcs=[tconcs;t,cis];endforendfunction

Contribuição dos produtos e reagentes para a derivada (reacção inversa e directa)

Ludwig Krippahl, 2007 36

Cinética com método de Euler

function tconcs=cinetica(esteq,cis,kd,ki,dt,tmax)rs=find(esteq<0);ps=find(esteq>0);tconcs=[0,cis];for t=0:dt:tmax

dps=prod(cis(ps).^esteq(ps))*ki;drs=prod(cis(rs).^-esteq(rs))*kd;deriv=(drs-dps)*dt;cis=cis+deriv*esteq;

tconcs=[tconcs;t,cis];endforendfunction

Calcular derivada e actualizar concentrações

Ludwig Krippahl, 2007 37

Cinética com método de Euler

function tconcs=cinetica(esteq,cis,kd,ki,dt,tmax)rs=find(esteq<0);ps=find(esteq>0);tconcs=[0,cis];for t=0:dt:tmax

dps=prod(cis(ps).^esteq(ps))*ki;drs=prod(cis(rs).^-esteq(rs))*kd;deriv=(drs-dps)*dt;cis=cis+deriv*esteq;

tconcs=[tconcs;t,cis];endforendfunction

Acrescentar à matriz o tempo t e as concentrações numa nova linha.

Ludwig Krippahl, 2007 38

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

esteq=[-1,-1,1] reacção

pontos=cinetica(esteq,cis,kd,ki,0.1,5);

Ludwig Krippahl, 2007 39

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, 2007 40

Cinética com método de Euler

Ludwig Krippahl, 2007 41

Sistema de reacções.

Mesma abordagem, mas estequiometria é matriz (uma linha por reacção) e cada constante cinética um vector (um valor por reacção)

Ludwig Krippahl, 2007 42

Sistema de reacções.

Alterações:• Dentro do ciclo for com os passos é preciso dois

ciclos, cada um percorrendo todas as reacções.

• O primeiro calcula os índices dos reagentes e produtos e o valor da derivada para cada reacção.

• O segundo ciclo actualiza as concentrações para todas as reacções considerando todas as derivadas calculadas.

É preciso dois ciclos internos para processar todas as reacções em simultâneo.

Ludwig Krippahl, 2007 43

Ajuste de um modelo

Dados Experimentais Simulação

Discrepância

Minimizar

Ludwig Krippahl, 2007 44

Ajuste de um modelo

Dados Experimentais Simulação

Discrepância

Minimizarminfn

cinetica

Ludwig Krippahl, 2007 45

Ajuste de um modelo

Dados: matriz com tempo na primeira coluna e concentração (ou concentrações) na segunda.

Função erro compara cada vector com o correspondente na simulação.

Mas os valores de t podem ser diferentes. É preciso interpolar.

Primeiro, função interpol

Ludwig Krippahl, 2007 46

Interpolação linear

Função interpol Recebe: uma matriz x, y, em colunas, e um

vector x1 com os pontos a interpolar. Devolve: vector y1 com os valores em x1

interpolados de x, y.

Ludwig Krippahl, 2007 47

Interpolação linear

xix1 x2

y1

y2

Ludwig Krippahl, 2007 48

Interpolação linear

yi = (y1*(x2-xi) + y2*(xi-x1)) / (x2 – x1)

xix1 x2

y1

y2

yi

Ludwig Krippahl, 2007 49

Interpolação linearfunction yi=interpol(matxy,xi)yi=0*xi;for f=1:length(xi)

for g=2:rows(matxy)if matxy(g,1)>=xi(f);

x1 = matxy(g-1,1);x2 = matxy(g,1);y1 = matxy(g-1,2);y2 = matxy(g,2);d = x2-x1;yi(f) = (y1*(x2-xi(f))+y2*(xi(f)-x1))/d;break

endifendfor

endfor

Ludwig Krippahl, 2007 50

Interpolação linearfunction yi=interpol(matxy,xi)yi=0*xi;for f=1:length(xi)

for g=2:rows(matxy)if matxy(g,1)>=xi(f);

x1 = matxy(g-1,1);x2 = matxy(g,1);y1 = matxy(g-1,2);y2 = matxy(g,2);d = x2-x1;yi(f) = (y1*(x2-xi(f))+y2*(xi(f)-x1))/d;break

endifendfor

endfor

Cria vector yi, dos valores interpolados

Ludwig Krippahl, 2007 51

Interpolação linearfunction yi=interpol(matxy,xi)yi=0*xi;for f=1:length(xi)

for g=2:rows(matxy)if matxy(g,1)>=xi(f);

x1 = matxy(g-1,1);x2 = matxy(g,1);y1 = matxy(g-1,2);y2 = matxy(g,2);d = x2-x1;yi(f) = (y1*(x2-xi(f))+y2*(xi(f)-x1))/d;break

endifendfor

endfor

Para cada xi onde interpolar percorre os x da matriz até encontrar o primeiro que ultrapassa xi. Começa do 2º elemento porque precisa do anterior para interpolar.

Ludwig Krippahl, 2007 52

Interpolação linearfunction yi=interpol(matxy,xi)yi=0*xi;for f=1:length(xi)

for g=2:rows(matxy)if matxy(g,1)>=xi(f);

x1 = matxy(g-1,1);x2 = matxy(g,1);y1 = matxy(g-1,2);y2 = matxy(g,2);d = x2-x1;yi(f) = (y1*(x2-xi(f))+y2*(xi(f)-x1))/d;break

endifendfor

endfor

Calcula a interpolação e termina o ciclo interno (g).

Ludwig Krippahl, 2007 53

Interpolação linear

xy=[[1:10]',[2:2:20]'];

xi=[2.5:2:8];

yi=interpol(xy,xi)

hold off

plot(xy(:,1), xy(:,2))

hold on

plot(xi,yi,"ob;;");

Ludwig Krippahl, 2007 54

Interpolação linear

Ludwig Krippahl, 2007 55

Medir a discrepância (erro)

Reacção• 2A B

• Só kd

Função erro mede o erro quadrático médio, que é a média dos quadrados das diferenças entre os vectores

Ludwig Krippahl, 2007 56

Medir a discrepância (erro)

Função erro mede o erro quadrático: soma dos quadrados das diferenças.

Ludwig Krippahl, 2007 57

Medir a discrepância (erro)

Exemplo:• 2A B

• Só kd (irreversível) Função erro2AB mede o erro quadrático

entre os dados experimentais e o cálculo para esta reacção.

A função inclui a concentração inicial e reacção

Ludwig Krippahl, 2007 58

Medir a discrepância (2A B)

function r=erro2AB(vals,k)esteq=[-2,1]; define a reacçãocis=[1,0]; e as concentrações

aqui falta calcular os valores previstos pelo modelo para este k e comparar com o vector vals para calcular o erro, interpolando os valores

endfunction

Ludwig Krippahl, 2007 59

Ajustar o modelo (2A B)

Basta usar a minfn para calcular o k que minimiza o erro

Exemplo:• vals=[0.5,0.5;2,0.2;6,0.07;9,0.055];

• k=minfn("erro2AB",vals,0,1,2,0.001)• k = 0.97843

Ludwig Krippahl, 2007 60

Ajustar o modelo (2A B)

Comparar o modelo com os dadosesteq=[-2,1];

cis=[1,0];

xy=cinetica(esteq,cis,k,0,0.01,10);

hold off

plot(xy(:,1),xy(:,2))

hold on

plot(vals(:,1),vals(:,2),"ob;;");

Ludwig Krippahl, 2007 61

Ajustar o modelo (2A B)

Ludwig Krippahl, 2007 62

Ajustar um modelo

Abordagem genérica• Simular dados previstos para um conjunto

de parâmetros

• Minimizar a discrepância entre os valores previstos e observados alterando os parâmetros.

• Na prática pode ser difícil...

Ludwig Krippahl, 2007 63

Dúvidas