|
|
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
lxyy: 金币+10, ★★★很有帮助 2012-03-30 12:30:33
function KineticsEst
clear all
clc
x0 = [4.639975 5.289842 0 0.039974 0.063806]; % A B C R W 组成的初始浓
k0=[1 1 1]; % k1, k2, k3 初始向值
lb = [0 0 0];
ub = [1 +inf +inf]; % 上下限
KineticsData3;
yexp = Kinetics(:,2:6);
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian]=...
lsqnonlin(@ObjFunc,k0,lb,ub,[],x0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f ± %.4f\n',k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2))
fprintf('\tk3 = %.4f ± %.4f\n',k(3))
% ------------------------------------------------------------------
function f = ObjFunc(k,x0,yexp) % 目标函数
tspan = [0;5;10;15;30;46;60;90;120;160;200;240;300;360];
[t x] = ode45(@Euqations,tspan,x0,[],k);
y(:,1) = x(:,1);
y(:,2:5) = x(:,2:5);
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f5 = y(:,5) - yexp(:,5);
f = [f1; f2; f3; f4;f5];
% ------------------------------------------------------------------
function dcdt = Euqations(t,x,k)
r1 = k(1)*x(1)*x(2)-k(1)*x(3)/3.94;
r2 = k(2)*x(1)*x(5);
r3 = k(3)*x(2)*x(4)-k(3)*x(3)*x(5)/1.08;
dcAdt = -r1-r2;
dcBdt = -r1-r3;
dcCdt = r1+r3;
dcRdt = r2-r3;
dcWdt = -r2+r3;
dcdt = [dcAdt;dcBdt; dcCdt;dcRdt;dcWdt;]; |
|
回复此楼