24小时热门版块排行榜

返回列表

飞鸿印雪jay

银虫 (小有名气)

应助: 1 (幼儿园)
金币: 340.5
散金: 60
帖子: 69
在线: 60.3小时
虫号: 2338094
注册: 2013-03-11
性别: GG
专业: 生物化工与食品化工

[求助] 我想求这个方程里的k，用这个代码但是有错误，求大神改改，万分感谢！！！

function k1k2k3
format long
clear all
clc
tspan = [0  30  50  80 140 170 200 230 260 290 320 360 400 460 520 580 600];
x0 = [9;0;0;0;0];
k0 = [1  1  1  1  1  1  1  1  1  1];
lb = [0  0  0  0  0  0  0  0  0  0];
ub = [1  1  1  1  1  1  1  1  1  1];

data=[30       7.939       1.458 7.939       1.458    7.939
   50       7.687       1.535 7.687       1.535    7.687
   80       7.289       1.602 7.289       1.602    7.289
   140       6.658       1.717 6.658       1.717    6.658
   170       6.531       1.722 6.531       1.722    6.531
   200       6.218       1.671 6.218       1.671    6.218
   230       5.979       1.620 5.979       1.620    5.979
   260       5.591       1.550 5.591       1.550    5.591
   290       5.414       1.488 5.414       1.488    5.414
   320       4.968       1.433 4.968       1.433    4.968
   360       4.692       1.350 4.692       1.350    4.692
   400       4.438       1.319 4.438       1.319    4.438
   460       4.144       1.294 4.144       1.294    4.144
   520       4.041       1.294 4.041       1.294    4.041
   580       4.052       1.287 4.052       1.287    4.052
   600       4.052       1.287 4.052       1.287    4.052];
yexp = data(:,2:6);

[k,resnorm,residual,exitflag,output,lambda,jacobian] =...
lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f ± %.9f\n',ci(1,2)-k(1),k(1))
fprintf('\tk2 = %.9f ± %.9f\n',ci(1,2)-k(2),k(2))
fprintf('\tk3 = %.9f ± %.9f\n',ci(1,2)-k(3),k(3))
fprintf('\tk3 = %.9f ± %.9f\n',ci(1,2)-k(4),k(4))
fprintf('\tk3 = %.9f ± %.9f\n',ci(1,2)-k(5),k(5))
fprintf('\tk3 = %.9f ± %.9f\n',ci(1,2)-k(6),k(6))
fprintf('\tk3 = %.9f ± %.9f\n',ci(1,2)-k(7),k(7))
fprintf('\tk3 = %.9f ± %.9f\n',ci(1,2)-k(8),k(8))
fprintf('\tk3 = %.9f ± %.9f\n',ci(1,2)-k(9),k(9))
fprintf('\tk3 = %.9f ± %.9f\n',ci(1,2)-k(10),k(10))
fprintf('  The sum of the squares is: %.9e\n\n',resnorm)

function f = ObjFunc(k,tspan,x0,yexp)          % 目标函数
[t, Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,3);
Xsim4=Xsim(:,4);
Xsim5=Xsim(:,5);
% Xsim6=Xsim(:,6);
ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);
ysim(:,4) = Xsim4(2:end);
ysim(:,5) = Xsim5(2:end);
% ysim(:,6) = Xsim6(2:end);
size(ysim(:,1));
size(ysim(:,2));
size(ysim(:,3));
size(ysim(:,4));
size(ysim(:,5));
% size(ysim(:,6));
size(yexp(:,1));
size(yexp(:,2));
size(yexp(:,3));
size(yexp(:,4));
size(yexp(:,5));
% size(yexp(:,6));
f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3)) (ysim(:,4)-yexp(:,4)) (ysim(:,5)-yexp(:,5))];

function dCdt = KineticsEqs(t,C,k)             % ODE模型方程
dCAdt = (k(1)+k(5)+k(6)+k(7))*C(1);
dCBdt = k(1)*C(1)-(k(2)+k(8)+k(9))*C(2);
dCCdt = k(5)*C(1)+k(2)*C(2)-(k(3)+k(10))*C(3);
dCDdt = k(6)*C(1)+k(9)*C(2)+k(3)*C(3)-k(4)*C(4);
dCEdt = k(7)*C(1)+k(8)*C(2)+k(10)*C(3)+k(4)*C(4);
dCdt = [dCAdt; dCBdt;dCCdt;dCDdt;dCEdt];

回复此楼

» 猜你喜欢

1楼 2013-11-11 19:40:49

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

相关版块跳转我要订阅楼主飞鸿印雪jay 的主题更新

返回列表