Programação para as
Ciências Experimentais
2007/8
Teórica 10
Ludwig Krippahl, 2008
Na aula de hoje...






Ajustar um modelo a dados
experimentais.
Interpolação linear
Minimização de funções
Cálculo de erros
Estimar uma constante cinética
ajustando o modelo aos dados.
Conceitos básicos de Excel
Ludwig Krippahl, 2008
2
Ajuste de um modelo
Dados
Experimentais
Simulação
Discrepância
Minimizar
Ludwig Krippahl, 2008
3
Ajuste de um modelo

Exemplo: reacção química
Dados
Experimentais
cinetica
Simulação
Discrepância
minfn
Ludwig Krippahl, 2008
Minimizar
4
Ajuste de um modelo




Dados: matriz com tempo na primeira coluna e
concentração (ou concentrações) na segunda
(ou outras).
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, 2008
5
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, 2008
6
Interpolação linear
y2
y1
x1
Ludwig Krippahl, 2008
xi
x2
7
Interpolação linear

yi = (y1*(x2-xi) + y2*(xi-x1)) / (x2 – x1)
y2
yi
y1
x1
Ludwig Krippahl, 2008
xi
x2
8
Interpolação linear
function 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
endif
endfor
endfor
Ludwig Krippahl, 2008
9
Interpolação linear
function yi=interpol(matxy,xi)
yi=0*xi;
for f=1:length(xi)
for g=2:rows(matxy)
if matxy(g,1)>=xi(f);
Cria
vector yi, dos valores interpolados
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
endif
endfor
endfor
Ludwig Krippahl, 2008
10
Interpolação linear
function 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);
Para
cada xi onde interpolar percorre
y2 = matxy(g,2);
os
da matriz até encontrar o primeiro
d =xx2-x1;
yi(f) =ultrapassa
(y1*(x2-xi(f))+y2*(xi(f)-x1))/d;
que
xi. Começa do 2º
break
elemento
porque precisa do anterior
endif
para interpolar.
endfor
endfor
Ludwig Krippahl, 2008
11
Interpolação linear
function yi=interpol(matxy,xi)
yi=0*xi;
for f=1:length(xi)
for g=2:rows(matxy)
if matxy(g,1)>=xi(f);
Calcula a interpolação
x1 = matxy(g-1,1);
e termina o ciclo
x2 = matxy(g,1);
y1 = matxy(g-1,2);
interno (g).
y2 = matxy(g,2);
d = x2-x1;
yi(f) = (y1*(x2-xi(f))+y2*(xi(f)-x1))/d;
break
endif
endfor
endfor
Ludwig Krippahl, 2008
12
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, 2008
13
Interpolação linear
Ludwig Krippahl, 2008
14
Medir a discrepância (erro)

Reacção

Função erro mede o erro quadrático
médio, que é a média dos quadrados
das diferenças entre os vectores
• 2A  B
• Só kd
Ludwig Krippahl, 2008
15
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 a simulação.
A função codifica a concentração inicial e
reacção, recebe como argumentos o kd e os
valores para comparar.
Ludwig Krippahl, 2008
16
Medir a discrepância (2A  B)
function r=erro2AB(vals,k)
er=[2,0];
define a reacção
ep=[0,1];
cis=[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. Para resolver na prática...
endfunction
Ludwig Krippahl, 2008
17
Medir a discrepância (2A  B)


Para simular a reacção podemos usar a
função cinetica da aula anterior.
Para comparar com os dados
experimentais precisamos interpolar
para os valores de t experimentais (que
podem não coincidir com os da
simulação)
Ludwig Krippahl, 2008
18
Medir a discrepância (2A  B)

O erro é o erro quadrático:
r=sum((vals(:,2)-int).^2);
vals é a matriz com as concentrações
de A na segunda coluna
int é o vector das concentrações de A
obtido interpolando a simulação para os
valores na 1ª coluna de vals.
Ludwig Krippahl, 2008
19
O mínimo de uma função

Método da razão dourada
Ludwig Krippahl, 2008
20
O mínimo de uma função

Tal como “encurralámos” a raiz num
intervalo, vamos fazer o mesmo com o
mínimo, mas precisamos de 3 pontos:
a
c
b
Ludwig Krippahl, 2008
21
O mínimo de uma função

Se x1<x2<x3 e y2<y1 e y2<y3 então tem
que haver um mínimo local entre x1 e x3
x1
Ludwig Krippahl, 2008
x2
x3
22
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2008
x2
x3
23
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2008
x2
x3
24
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2008
x2
x3
25
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2008
x2
x3
26
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2008
x2
x3
27
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2008
x2
x3
28
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2008
x2
x3
29
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2008
x2
x3
30
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2008
x2
x3
31
O mínimo de uma função

Guardar sempre os 3 pontos
consecutivos em que o y do meio é
menor que os extremos.
x1
Ludwig Krippahl, 2008
x2
x3
32
O mínimo de uma função

Como dividir o intervalo:
• O ideal é manter as proporções. Dividir ao
meio não é ideal.
x1
Ludwig Krippahl, 2008
x3
x2
33
O mínimo de uma função

Como dividir o intervalo:
• O ideal é manter as proporções. Dividir ao
meio não é ideal.
x1
Ludwig Krippahl, 2008
x4
x5
x3
x2
34
O mínimo de uma função

Como dividir o intervalo:
• Escolher o ponto novo no intervalo maior e
• Partir pela razão dourada:
(a+b)/a = a / b
a= 0.618 (a+b)
b= (1-0.618) (a+b)
Ludwig Krippahl, 2008
35
O mínimo de uma função
function xm=minfn(func,params,x1,xm,x2,prec)
c=1-0.618;
Nome
da
função,
ym=feval(func,params,xm);
parâmetros (como no
zerpol), os 3 pontos
iniciais e precisão
Ludwig Krippahl, 2008
36
O mínimo de uma função
function xm=minfn(func,params,x1,xm,x2,prec)
c=1-0.618;
ym=feval(func,params,xm);
Constante c para os
intervalos (razão dourada)
Ludwig Krippahl, 2008
37
O mínimo de uma função
function xm=minfn(func,params,x1,xm,x2,prec)
c=1-0.618;
ym=feval(func,params,xm);
Avalia a função no ponto
do meio. Nota: assume-se
que y é maior em x1 e x2.
Ludwig Krippahl, 2008
38
O mínimo de uma função
while abs(x2-x1)>prec
if abs(x1-xm)>abs(x2-xm)
intervalo maiorEnquanto
é x1 a xmo intervalo é
else
maior que a precisão
intervalo maior é xm a x2
endif
endwhile
Ludwig Krippahl, 2008
39
O mínimo de uma função
while abs(x2-x1)>prec
if abs(x1-xm)>abs(x2-xm)
intervalo maior é x1 a xm
else
intervalo maior é xm a x2
Encontra o sub-intervalo
endif
maior, (x1 a xm ou xm a
endwhile
x2)
Ludwig Krippahl, 2008
40
O mínimo de uma função
x1
Ludwig Krippahl, 2008
xm
x2
41
O mínimo de uma função
Se o intervalo maior é de x1 a xm
o novo x será entre x1 e xm, próximo de xm
xn=xm-c*(xm-x1)
o novo y será feval(func,params,xn)
Se o novo y for menor que o anterior (em xm)
passar o x2 para onde está xm, xm para o
novo x, e ym será o novo y.
Ludwig Krippahl, 2008
42
O mínimo de uma função
ym
yn
x1
Ludwig Krippahl, 2008
xn
xm
x2
43
O mínimo de uma função
ym
x1
Ludwig Krippahl, 2008
xm
x2
x2
44
O mínimo de uma função
Se o intervalo maior é de xm a x2
o novo x será entre xm e x2, mais próximo
de xm.
xn=xm+c*(x2-xm);
Se o novo y for menor que o anterior (em xm)
passar o x1 para onde está xm, xm para o
novo x, e ym será o novo y.
Ludwig Krippahl, 2008
45
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)
•
Ludwig Krippahl, 2008
k = 0.97843
46
Ajustar o modelo (2A  B)

Comparar o modelo com os dados
er=[2,0]
ep=[0,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), "x");
Ludwig Krippahl, 2008
47
Ajustar o modelo (2A  B)
Ludwig Krippahl, 2008
48
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, 2008
49
Conceitos básicos de Excel



Célula: A5
Grupo de células: A5:B12
Referência relativa ou absoluta:
• O cifrão marca uma referência absoluta.
• A$5, $B$5 Nestes casos o 5 e o B estão
•
fixos.
Sem cifrão a referência é relativa, e muda
com copy/paste ou fill down/right
Ludwig Krippahl, 2008
50
Conceitos básicos de Excel

Referência relativa:
• Nota: fórmulas começam sempre por =
Ludwig Krippahl, 2008
51
Conceitos básicos de Excel

Referência relativa:
• O B passou a C e o C a D copiando para a
direita
Ludwig Krippahl, 2008
52
Conceitos básicos de Excel

Referência relativa:
• O 2 passou a 3 copiando para baixo
Ludwig Krippahl, 2008
53
Conceitos básicos de Excel

Referência absoluta
Ludwig Krippahl, 2008
54
Conceitos básicos de Excel

Referência absoluta
• Fill down (seleccionar, ctrl+d)
Ludwig Krippahl, 2008
55
Conceitos básicos de Excel

Referência absoluta
• Multiplicar pelo C1, mas sem mudar o 1...
Ludwig Krippahl, 2008
56
Conceitos básicos de Excel

Referência absoluta
• Marcar o 1 como ref. absoluta
Ludwig Krippahl, 2008
57
Conceitos básicos de Excel

Referência absoluta
• Marcar o 1 como ref. absoluta
Ludwig Krippahl, 2008
58
Conceitos básicos de Excel

Dar nomes às células.
• Exemplo: 2A  B
• Parâmetros
• Constante
• DeltaT
Ludwig Krippahl, 2008
59
Conceitos básicos de Excel

Dar nomes às células.
• Exemplo: 2A  B
• Parâmetros
• Constante
• DeltaT
Ludwig Krippahl, 2008
60
Conceitos básicos de Excel

Dar nomes às células.
• Exemplo: 2A  B
• Parâmetros
• Constante
• DeltaT
Ludwig Krippahl, 2008
61
Conceitos básicos de Excel

Dar nomes às células.
• Exemplo: 2A  B
• Parâmetros
• Constante
• DeltaT
Ludwig Krippahl, 2008
62
Conceitos básicos de Excel

Dar nomes às células.
• Exemplo: 2A  B
• Parâmetros
• Constante
• DeltaT
Ludwig Krippahl, 2008
63
Conceitos básicos de Excel

Dar nomes às células.
• Exemplo: 2A  B
• Parâmetros
• Constante
• DeltaT
Ludwig Krippahl, 2008
64
Conceitos básicos de Excel

Dar nomes às células.
• Exemplo: 2A  B
• Parâmetros
• Constante
• DeltaT
Ludwig Krippahl, 2008
65
Conceitos básicos de Excel

Dar nomes às células.
• Exemplo: 2A  B
• Parâmetros
• Constante
• DeltaT
Ludwig Krippahl, 2008
66
Conceitos básicos de Excel

Dar nomes às células.
• Exemplo: 2A  B
• Fill down...
• Mas falta o tempo.
Ludwig Krippahl, 2008
67
Conceitos básicos de Excel

Seleccionar a primeira coluna (click no
topo da coluna, no A).
Ludwig Krippahl, 2008
68
Conceitos básicos de Excel

Insert, Columns
Ludwig Krippahl, 2008
69
Conceitos básicos de Excel

Insert, Columns
Ludwig Krippahl, 2008
70
Conceitos básicos de Excel

Definir a fórmula, e fill down.
Ludwig Krippahl, 2008
71
Conceitos básicos de Excel


IF(condição; valor se verdade; valor se
falso)
Ex:
Ludwig Krippahl, 2008
72
Conceitos básicos de Excel


IF(condição; valor se verdade; valor se
falso)
Ex:
Ludwig Krippahl, 2008
73
Conceitos básicos de Excel

Exemplo: raiz do polinómio x3+2
Ludwig Krippahl, 2008
74
Conceitos básicos de Excel

Exemplo: raiz do polinómio x3+2
Ludwig Krippahl, 2008
75
Conceitos básicos de Excel

Exemplo: raiz do polinómio x3+2
Fill right, fill down
Ludwig Krippahl, 2008
76
Conceitos básicos de Excel

Exemplo: raiz do polinómio x3+2
Ludwig Krippahl, 2008
77
Conceitos básicos de Excel

Exemplo: raiz do polinómio x3+2
Fill down
Ludwig Krippahl, 2008
78
Dúvidas
Ludwig Krippahl, 2008
79