%**************************************************************************
% Programa para demonstrar o método integral (ordem <=2) por derivação de um
% Conjunto de datos na forma de C(t).
%
% (c)FGF 2008
%**************************************************************************
clear; clc; clf;    %limpar tudo
syms x
%**************************************************************************
%       DADOS
%**************************************************************************
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];
plot(dados.t,dados.C,'o');
xlabel('t ()');
ylabel('C (M)');
hold on;
%**************************************************************************
N=length(dados.t);
%**************************************************************************
%   Ordem 0
%**************************************************************************
R0=corrcoef(dados.t,dados.C)
if R0(1,2)^2<0.99
%**************************************************************************
%   Ordem 1
%**************************************************************************
    R1=corrcoef(dados.t,log(dados.C(1)./dados.C))
    if R1(1,2)^2>R0(1,2)^2
%**************************************************************************
%   Ordem 2
%**************************************************************************
        R2=corrcoef(dados.t,(1./dados.C))
        if R2(1,2)^2<R1(1,2)^2
             r2=R1(1,2)^2;S=strcat('A ordem é 1, R²=',num2str(r2));R=1;
        else r2=R2(1,2)^2;S=strcat('A ordem pode ser maior que 2, R²=',num2str(r2));R=2;
        end
    else
        r2=R0(1,2)^2;S=strcat('A ordem é 0, R²=',num2str(r2));R=0;
    end
else
r2=R0(1,2)^2;S=strcat('A ordem é 0, R²=',num2str(r2));R=0;
end
S
if R==0;y=dados.C;end
if R==1;y=log(dados.C(1)./dados.C);end
if R==2;y=1./dados.C;end
para=polyfit(dados.t,y,1);    % ajuste linear
if R==0;k=-para(1);CAo=para(2);end
if R==1;k=para(1);CAo=dados.C(1);end
if R==2;k=para(1);CAo=1/para(2);end
%**************************************************************************
% 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/(k*x^R),x);                                %integrar -dc/kC^n=dt
y=dados.C(end):dados.C(1)/100:dados.C(1);                % Vector de C's
to=subs(z,CAo);                                   % limite inferior
t=subs(z,y)-to;                                          % t=P-limite inf.
%**************************************************************************
% Escrever ordem,erro e ko no gráfico e na linha de comando
%**************************************************************************
S=strcat('N=', num2str(R), '\newline ko=',num2str(k));% escrever...
%**************************************************************************
% Acabar o gráfico representando o ajuste pela lei cinética obtida
%**************************************************************************
plot(t,y,'r');
S1=strcat('R²=',num2str(r2));
text(dados.t(int8((N+1)/2)),dados.C(int8((N)/2)),S1,'HorizontalAlignment','left')
legend('pontos experimentais',S);

R0 =

    1.0000   -0.8698
   -0.8698    1.0000


R1 =

    1.0000    0.9500
    0.9500    1.0000


R2 =

    1.0000    0.9936
    0.9936    1.0000


S =

A ordem pode ser maior que 2, R²=0.98722

Published with MATLAB® 7.5

dados.t=[0 2 4 6 8 10];
dados.C=[3 2.5 2 1.5 1 0.5];

ou

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