| 查看: 4386 | 回复: 14 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[求助]
matlab 拟合反应动力学参数结果很差。大家帮忙看一下
|
|||
|
已知动力学模型为： dA1dt = r1 + r2; dA2dt = - r1-r2; dA3dt = -r1； dA4dt = r1 - r2; dA5dt = r2; r1 = k(1)*((a(3)*a(2) - (1/2.20)*a(4)*a(1))); r2 = k(2)*(a(4)*a(2) - (1/0.40)*a(5)*a(1))); 已知的参数： 时间：tspan = [0 60 120 180 240 300 360 420 600 1200 6900 7500], 活度系数：KineticsData1 % 动力学数据 % a1 a2 a3 a4 a5 ExpData=... [0.0000 0.8122 0.3894 0.0000 0.0000 0.1827 0.6353 0.1599 0.2163 0.0444 0.2547 0.5668 0.0979 0.2485 0.0805 0.2854 0.5390 0.0788 0.2415 0.1061 0.2978 0.5258 0.0738 0.2359 0.1192 0.3019 0.5210 0.0710 0.2309 0.1288 0.3092 0.5143 0.0691 0.2264 0.1349 0.3178 0.5062 0.0680 0.2223 0.1349 0.3178 0.5062 0.0680 0.2186 0.1433 0.3211 0.5030 0.0674 0.2154 0.1478 0.3209 0.5036 0.0673 0.2140 0.1482 0.3205 0.5045 0.0674 0.2136 0.1476 ] 下面是写的matlab程序： function Kinetic2 % 动力学ODE方程模型的参数估计 clear all clc k0 = [900000 50000]; % 参数初值 lb = [0 0]; % 参数下限 ub = [+inf +inf ]; % 参数上限 a0 = [0.0000 0.8122 0.3894 0.0000 0.0000]; KineticsData1; yexp = ExpData(:,1:5); % yexp: 实验数据[a1 a2 a3 a4 a5] % 使用函数fmincon()进行参数估计 [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],a0,yexp); fprintf('\n使用函数fmincon()估计得到的参数值为:\n') fprintf('\tk1 = %.4f\n',k(1)) fprintf('\tk2 = %.4f\n',k(2)) fprintf(' The sum of the squares is: %.1e\n\n',fval) k_fmincon = k; % 使用函数lsqnonlin()进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],a0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.4f\n',k(1)) fprintf('\tk2 = %.4f\n',k(2)) fprintf(' The sum of the squares is: %.1e\n\n',resnorm) exitflag % 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计 k0 = k_fmincon; [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[1e-15],a0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.4f\n',k(1)) fprintf('\tk2 = %.4f\n',k(2)) fprintf(' The sum of the squares is: %.1e\n\n',resnorm) exitflag function f = ObjFunc4Fmincon(k,a0,yexp) tspan = [0 60 120 180 240 300 360 420 600 1200 6900 7500]; [t a] = ode23s(@KineticEqs,tspan,a0,[],k); y(:,1:5) = a(:,1:5); f = sum(((y(:,1)-yexp(:,1))/0.01).^2) + sum(((y(:,2)-yexp(:,2))/0.01).^2) ... + sum(((y(:,3)-yexp(:,3))/0.005).^2) + sum(((y(:,4)-yexp(:,4))/0.01).^2) ... + sum(((y(:,5)-yexp(:,5))/0.01).^2); function f = ObjFunc4LNL(k,a0,yexp) tspan = [0 60 120 180 240 300 360 420 600 1200 6900 7500]; [t a] = ode23s(@KineticEqs,tspan,a0,[],k); y(:,1:5) = a(:,1:5); f1 = (y(:,1) - yexp(:,1))/0.01; f2 = (y(:,2) - yexp(:,2))/0.01; f3 = (y(:,3) - yexp(:,3))/0.005; f4 = (y(:,4) - yexp(:,4))/0.01; f5 = (y(:,5) - yexp(:,5))/0.01; f = [f1; f2; f3; f4; f5]; plot(tspan,yexp(:,1), 'b^',tspan,yexp(:,2), 'ro',tspan,yexp(:,3),'y*',... tspan,yexp(:,4), 'go',tspan,yexp(:,5), 'cs') hold on plot(tspan,y(:,1), 'b-',tspan,yexp(:,2), 'r-',tspan,yexp(:,3),'y-', ... tspan,yexp(:,4), 'g-',tspan,yexp(:,5), 'c-') axis([0 800 0 1]) %动力学模型----------------------------------------------------------------------- function dadt = KineticEqs(t,a ,k) Keq=[2.20 0.40]; dadt=... [((k(1)*(a(3)*a(2)-(1/Keq(1))*a(4)*a(1))) + (k(2)*(a(4)*a(2)-(1/Keq(2))*a(5)*a(1)))); -((k(1)*(a(3)*a(2)-(1/Keq(1))*a(4)*a(1))) + (k(2)*(a(4)*a(2)-(1/Keq(2))*a(5)*a(1)))); -(k(1)*(a(3)*a(2)-(1/Keq(1))*a(4)*a(1))); ((k(1)*(a(3)*a(2)-(1/Keq(1))*a(4)*a(1))) - (k(2)*(a(4)*a(2)-(1/Keq(2))*a(5)*a(1)))); (k(2)*(a(4)*a(2)-(1/Keq(2))*a(5)*a(1)))]; 大家帮忙看一看，哪里有问题，从做出的图可以看出，处理的结果很差啊。 无论K初值给出怎么样，都是拟合出来的a1和实验值a1相差很远，而其他值都拟合很好。  [ Last edited by 楚天湘水 on 2012-5-1 at 00:44 ] |
回复此楼» 猜你喜欢
散金币祈福
已经有186人回复
-
求助，如何提取ELK的rt-TDDFT在某一时刻的自旋密度分布和电子密度分布
已经有1人回复
-
物理学I论文润色/翻译怎么收费?
已经有140人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
深圳大学应用技术学院招聘凝聚态物理博士后
已经有114人回复
-
深圳大学2027级光电信息工程/物理学研究生推免
已经有70人回复
-
求助中医院临床中医中标的青年标书学习，谢谢
已经有4人回复
-
浙大物院表面物理课题组长期招收普博及直博学生
已经有257人回复
-
【征稿】Springer Nature热电专题 Discover Materials / Discover Chemistry
已经有39人回复
|
本帖内容被屏蔽 |
14楼2021-04-24 21:39:10
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
2楼2012-05-02 10:22:42
3楼2012-05-02 10:42:50
hitzhengwei
木虫 (著名写手)
- 应助: 34 (小学生)
- 金币: 1530.2
- 散金: 100
- 红花: 1
- 帖子: 1921
- 在线: 76.5小时
- 虫号: 1730549
- 注册: 2012-04-01
- 性别: GG
- 专业: 爆炸与冲击动力学
4楼2012-05-02 10:57:15
