resid matlab,小弟初学matlab，有大神解读下下这个例子么？？ - 计算模拟 - 小木虫 - 学术 科研 互动社区...

CODE:

function KineticsEst1_int1

% 动力学参数辨识: 用积分法进行反应速率分析得到速率常数k和反应级数n

% Analysis of kinetic rate data by using the integral 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

global CAm

t = [2.5 5 10 15 20 30 40 ];

CAm = [22.7063 13.4063 4.9063 1.6063 0.9063 0.9063 0.3063 ];

% 非线性拟合

beta0 = [0.1656  0.6832];

tspan = [2.5 5 10 15 20 30 40 ];

CA0 = 22.7063;

[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...

lsqnonlin(@OptObjFunc,beta0,[],[],[],tspan,CA0,CAm)

ci = nlparci(beta,resid,jacobian)

% 拟合效果图(实验与拟合的比较)

[t4plot CA4plot] = ode45(@KineticsEqs,[tspan(1)  tspan(end)],CA0,[],beta);

plot(t,CAm,'bo',t4plot,CA4plot,'k-')

legend('Exp','Model')

xlabel('时间t, s')

ylabel('浓度C_A, mol/L')

% 残差关于拟合值的残差图

[t CAc] = ode45(@KineticsEqs,tspan,CA0,[],beta);

figure

plot(CAc,resid,'*')

xlabel('浓度拟合值(mol/L)')

ylabel('残差R (mol/L)')

refline(0,0)

% 参数辨识结果

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))

% ------------------------------------------------------------------

function f = OptObjFunc(beta,tspan,CA0,CAm)

[t CAc] = ode45(@KineticsEqs,tspan,CA0,[],beta);

f = CAc - CAm;

% ------------------------------------------------------------------

function dCAdt = KineticsEqs(t,CA,beta)

dCAdt = -beta(1)*CA^beta(2);            % k = beta(1), n = beta(2)