%************************************************************************** % Programa para demonstrar o método diferencial por derivação de um % polinómio de ajuste a um conjunto de datos C(t). % (c)FGF 2008 %************************************************************************** clear; clc; clf; %limpar tudo syms x z; %precisamos de variáveis simbólicas %************************************************************************** % DADOS %************************************************************************** dados.t=[0 2.5 5 7.5 10 15 20];
dados.C=[0.12 0.0964 0.0802 0.065 0.054 0.037 0.0248]; subplot(2,1,2);
plot(dados.t,dados.C,'o');
hold on; %************************************************************************** N=length(dados.t); % Nº de pontos experimentais
ordem=int8(N/2.5); if ordem>4; ordem=4;end %ordem do poliómio de ajuste
pC=polyfit(dados.t,dados.C,ordem); % ajuste polinomial aos dados
t=dados.t(1):dados.t(end)/(30):dados.t(end);% vector de t's C=polyval(pC,t); % vector de C's ajustados
plot(t,C);
xlabel('t ()'); ylabel('C (M)');
%************************************************************************** % Obtenção de dC/dt
%**************************************************************************
dcdt=-polyval(polyder(pC),t); %derivação do polinómio
%************************************************************************** % Método diferencial
% obtenção dos xlabel('t ()');
ylabel('C (M)'); %parâmetros de ajuste linear de ln(dC/dt)=ln(ko)+n ln(C)
%************************************************************************** l_dcdt=log(dcdt);
l_C=log(C);
para=polyfit(l_C,l_dcdt,1); % ajuste linear
M=int8(para(1)); e=double(int8((para(1)-double(M))*10))/10; % qual o erro disso %************************************************************************** % Vamos ver agora se a ordem obtida é razoável % uma das formas é representar a lei obtida no gráfico que é C(t)
% para isso é preciso integrar o balanço de massa % que separando variáveis fica dt=-dC/(-r)
z=int(-1/(exp(para(2))*x^para(1)),x); %integrar -dc/kC^n=dt
y=dados.C(end):dados.C(1)/100:dados.C(1); % Vector de C's
to=subs(z,dados.C(1)); % limite inferior
t=subs(z,y)-to; % t=P-limite inf.
%************************************************************************** % Acabar o gráfico representando o ajuste pela lei cinética obtida
%************************************************************************** S=strcat('ajuste polinomial de ordem = ',num2str(ordem));
plot(t,y,'r'); legend('pontos experimentais',S,'lei cinética obtida');
subplot(2,1,1);
for i=1:N-1
dados.log_C(i)=log((dados.C(i)+dados.C(i+1))/2);
dados.log_r(i)=log((dados.C(i+1)-dados.C(i))/(dados.t(i)-dados.t(i+1)));
end
plot(dados.log_C,dados.log_r,'o', l_C,polyval(para,l_C));
legend('pontos experimentais', 'lei cinética obtida','Location', 'NorthWest');
xlabel('ln(C)'); ylabel('ln(- \Delta C/ \Delta t)'); %************************************************************************** % Escrever ordem,erro e ko no gráfico e na linha de comando %************************************************************************** S=strcat('\leftarrow N=', num2str(M),'\pm', num2str(e), ' \wedge ko=',num2str(exp(para(2))));% escrever... text(dados.log_C(int8((N+1)/2)-1),dados.log_r(int8((N)/2)),S,'HorizontalAlignment','left')
strcat('N=',num2str(para(1)), ' e ',' ln(ko)=',num2str(para(2)),' => ko=',num2str(exp(para(2))))
ans =
N=0.9954 e ln(ko)=-2.5513 => ko=0.077982
