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.


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


Integração numérica.

y = exp(-x3)


Integração numérica

Aproximar a função considerando cada rectângulo

dx*y



Quanto menor dx mais preciso dx*y



function int=intexpxcubo(dx,x0,x1);

int=0;

for x=x0:dx:x1

int=int+dx*exp(-x^3);

endfor

endfunction



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



Aproximar melhor pela regra do trapézio



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

a

b

base



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



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



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



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

que integra.



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


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



Implementação mais genérica:

octave:185> trapezio("expxcubo",0.2,0,2)

ans = 0.89293

Nota: “expxcubo” em vez de expxcubo!


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.


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


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



Inicio, t0

• [A]0 [B]0 [C]0

Passo 1



• [A]1 = [A]0 + d[A]/dt * passo

Próximo valor



Inicio, t0

• [A]0 [B]0 [C]0

Passo 1



• [A]1 = [A]0 + d[A]/dt * passo

Valor anterior



Inicio, t0

• [A]0 [B]0 [C]0

Passo 1



• [A]1 = [A]0 + d[A]/dt * passo

Derivada vezes passo.


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



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




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


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


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



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



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

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

ans =

2 4 6






Encontrar os índices dos reagentes e produtos na estequiometria






Guardar valores para t0






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






Calcular derivada e actualizar concentrações






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



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);



Reacção do tipo• A + B C

hold off

axiseixo automático

plot(pontos(:,1),pontos(:,2:columns(pontos)))




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)


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.


Ajuste de um modelo

Dados Experimentais Simulação

Discrepância

Minimizar


Ajuste de um modelo

Dados Experimentais Simulação

Discrepância

Minimizarminfn

cinetica


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


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.



xix1 x2

y1

y2



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

xix1 x2

y1

y2

yi


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





endifendfor

endfor

Cria vector yi, dos valores interpolados





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.





endifendfor

endfor

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



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;;");




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



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



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


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


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



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;;");




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


Dúvidas

Documents

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