clear; clc; clf;
global dados n k
dados.t=[0 2 5 10 20 30 40 60 80 100];
dados.C=[0.2 0.175 0.15 0.123 0.093 0.077 0.065 0.051 0.044 0.039];
n=1;
k=-1/dados.t(2)*log(dados.C(2)/dados.C(1));
options=optimset('Display','iter','PlotFcns',@outfun);
para=fminsearch(@minq,[n k],options);
n=para(1);k=para(2);
%Subrotina de calculo do erro, função a minimizar, neste caso desvios quadrados
function erro=minq(x)
global dados n k
n=x(1);k=x(2);
erro=0;
for i=2:length(dados.t)
[t,C]=ode45(@BM,[0 dados.t(i)],dados.C(1));
erro=erro+(dados.C(i)-C(end))^2;
end
%Subrotina com o modelo, neste caso o modelo é o balanço de massas a um reactor descontínuo com uma
% uma reacção elementar de ordem n
function dydt = BM(t,y)
global n k
dydt = -k*y^n;
end
%Subrotina para fazer os gráficos
function stop = outfun(x,optimValues,state)
global n k dados
stop=false;
switch state
case 'init'
%plot(dados.t,dados.C,'o');
%hold on
%xlabel('t ()');
%ylabel('C (M)');
%[t,C]=ode45(@BM,[0 dados.t(end)],dados.C(1));
%plot(t,C,'g')
case 'iter'
set(gca,'NextPlot','replacechildren')
plot(dados.t,dados.C,'o');
hold on
[t,C]=ode45(@BM,[0 dados.t(end)],dados.C(1));
plot(t,C,'g')
S=strcat('\leftarrow N=', num2str(n),' \wedge ko=',num2str(k));
text(dados.t(int8(length(dados.t)/2+1)),dados.C(int8(length(dados.t)/2)),S,'HorizontalAlignment','left')
case 'done'
hold off
plot(dados.t,dados.C,'o');
hold on
xlabel('t ()');
ylabel('C (M)');
[t,C]=ode45(@BM,[0 dados.t(end)],dados.C(1));
plot(t,C,'r')
strcat('N=',num2str(n),' e ',' ko=',num2str(k))
S=strcat('\leftarrow N=', num2str(n),' \wedge ko=',num2str(k));% escrever...
legend('pontos experimentais','lei cinética obtida');
text(dados.t(int8(length(dados.t)/2+1)),dados.C(int8(length(dados.t)/2)),S,'HorizontalAlignment','left')
hold off
otherwise
end
end
Iteration Func-count min f(x) Procedure
0 1 0.0128089
1 3 0.0102973 initial simplex
2 5 0.00741204 expand
:
:
63 120 2.89856e-06 contract inside
64 122 2.89856e-06 contract inside
65 124 2.89856e-06 contract inside
66 126 2.89838e-06 contract inside
67 128 2.8983e-06 contract inside
68 130 2.89826e-06 contract inside
69 132 2.89826e-06 contract outside
70 134 2.89824e-06 contract inside
71 135 2.89824e-06 reflect
72 137 2.89824e-06 contract inside
73 139 2.89824e-06 contract inside
74 141 2.89823e-06 contract inside
75 142 2.89823e-06 reflect
76 144 2.89823e-06 contract inside
77 146 2.89823e-06 contract outside
Optimization terminated:
the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04
and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04
ans =
N=2.4885 e ko=0.78723
