2楼: Originally posted by dingd at 2014-03-12 09:56:13
k1, k2, k3要求都大于0的话效果会很差的。

3楼: Originally posted by 805592836_01 at 2014-03-12 10:02:43
嗯，都要大于0的……...

4楼: Originally posted by dingd at 2014-03-12 10:35:53
不可能有好结果，检查下模型吧。...

function test_123
clear all;clc

format long

tspan=[0,180,300,420,600,900,1200,1800,2400];

yexp=[1.414,1.237,1.065,0.873,0.629,0.440,0.226,0.124]';




y0=1.845;

k0=[-1 1 1];
lb=-[inf 0 0];
ub=+[inf inf inf];        

yy=[y0 yexp'];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1=%.4f\n',k(1)),  %本征动力学参数
fprintf('\tk2=%.4f\n',k(2)),  % A的吸附平衡常数KA
fprintf('\tk3=%.4f\n',k(3)),  % H2的吸附平衡常数KB
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,y0,[],k);
[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k);
y=XXsim(2:end);
xexp=yexp;
R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2);
fprintf('\n\t相关系数之平方R^2 = %.6f',R2);
plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best'),


%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ;
ysim = Xsim(2:end);
f=ysim-yexp;
%----------------------------------------------------------

function dxdt = KineticsEqs(t,y,k)
y0=1.845;
CB=y;
CB0=0.4;
CA0=y0;
k1=k(1);
KA=k(2);
KB=k(3);
dxdt=k(1)*k(2)*k(3)*CA0*(CB.^3)./(1+k(2)*CA0+k(3)^(1/6)*CB.^(1/2)).^7;