function KineticsEst1_Diff  
% 動力學參數辨識: 用微分法進行反應速率分析得到速率常數k和反應級數n
% Analysis of kinetic rate data by using the differential method
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/07/27 $
%
% Reaction of the type -- rate = kCA^order
% order - reaction order
% rate -- reaction rate vector
% CA -- concentration vector for reactant A
% T -- vector of reaction time
% N -- number of data points
% k- reacion rate constant

clear all
clc

% 動力學數據
t = [0 20  40  60  120  180  300];
CAm = [10  8  6  5  3  2  1];

% 用最小二乘樣條擬合法計算微分dCA/dt--使用不經過實驗點的B樣條插值函數
knots = 3;
K = 3;                  % 三次B樣條
sp = spap2(knots,K,t,CAm);
pp = fnder(sp);         % 計算B樣條函數的導函數
dCAdt = fnval(pp,t);    % 計算t處的導函數值
rAm = dCAdt;

% 繪制濃度擬合曲線
ti = linspace(t(1),t(end),200);
CAi = fnval(sp,ti);
plot(t,CAm,'ro',ti,CAi,'b-')
xlabel('t')
ylabel('C_A')
legend('實驗值','B樣條擬合')

% 非線性擬合
beta0 = [0.0053 1.39];
[beta,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@OptObjFunc,beta0,[],[],[],rAm,CAm);          
ci = nlparci(beta,residual,jacobian);

% 參數辨識結果
fprintf('Estimated Parameters:\n')
fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))
fprintf('  The sum of the squares is: %.1e\n\n',sum(residual.^2))

% 繪制反應速率擬合曲線
figure
plot(t,rAm,'ro',t,Rate(CAm,beta),'b*')
xlabel('t')
ylabel('dC_Adt')
legend('Experiment','Kinetic Model')

% ------------------------------------------------------------------
function f = OptObjFunc(beta,rAm,CAm)
rAc = Rate(CAm,beta);
f = rAc - rAm;

% ------------------------------------------------------------------
function rA = Rate(CA,beta)
rA = -beta(1)*CA.^beta(2);   % -rA = -dCAdt = k*CA^n, 其中k=beta(1), n=beta(2)