| 查看: 3594 | 回复: 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 ] |
回复此楼» 收录本帖的淘帖专辑推荐
 动力学拟合 |
» 猜你喜欢
求国际会议网站
已经有1人回复
-
求取一些关于纳米材料和纳米技术相关的英文PPT。
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有166人回复
-
【复旦大学】二维材料方向招收2026年博士研究生1名
已经有0人回复
-
北京纳米能源与系统研究所 王中林院士/曹南颖研究员课题组2026级硕/博/博后招生
已经有10人回复
-
荷兰Utrecht University超快太赫兹光谱王海教授课题招收2026 CSC博士生
已经有19人回复
-
反铁磁体中的磁性切换：两种不同的机制已成功可视化
已经有0人回复
-
26申博推荐：南京航空航天大学国际前沿院光学方向招收博士生！
已经有0人回复
-
求标准粉末衍射卡号 ICDD 01-076-1802
已经有0人回复
» 本主题相关价值贴推荐，对您同样有帮助:
Kinetic Monte carlo动力学蒙特卡罗模拟化学反应。。
已经有17人回复- 对反应动力学参数的拟合及表观反应速率常数的计算
已经有5人回复
- 动力学参数拟合
已经有3人回复
- 求高手相助！matlab求反应动力学模型参数
已经有6人回复
- 这样的动力学方程能够有1stopt拟合出参数来么
已经有4人回复
- 使用MATLAB拟合动力学参数问题
已经有18人回复
- Matlab同时拟合活化能，反应级数，速率常数
已经有26人回复
- 使用matlab最优化方法拟合获得多个动力学参数中的问题
已经有4人回复
- 求大神帮忙matlab拟合函数求参数
已经有9人回复
- 求助用matlab拟合动力学方程
已经有9人回复
- 求助用matlab拟合动力学方程
已经有13人回复
- 关于matlab拟合
已经有4人回复
- 反应动力学的拟合有什么意义？
已经有12人回复
- 谁有成熟的求动力学参数的MATLAB的程序。
已经有11人回复
- 动力学拟合参数
已经有6人回复
- 非线性拟合估计方程参数（matlab，spss，origin都可以）
已经有4人回复
- 【求助】用matlab最优化方法进行参数拟合
已经有17人回复
- 【求助】催化反应动力学matlab计算各基元反应的速率常数时，该如何避免较小量被忽略？
已经有3人回复
feng198800
金虫 (正式写手)
- 应助: 15 (小学生)
- 金币: 860.2
- 散金: 234
- 红花: 1
- 帖子: 415
- 在线: 242.1小时
- 虫号: 700353
- 注册: 2009-02-12
- 性别: GG
- 专业: 化学反应工程
10楼2012-05-29 11:50:35
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.5小时
- 虫号: 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
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.5小时
- 虫号: 291104
- 注册: 2006-10-28
5楼2012-05-02 10:57:30
zt_chem
木虫 (正式写手)
- 应助: 12 (小学生)
- 金币: 2147.8
- 红花: 6
- 帖子: 872
- 在线: 981.6小时
- 虫号: 286920
- 注册: 2006-10-15
- 性别: GG
- 专业: 催化化学
6楼2012-05-02 21:49:24
7楼2012-05-03 00:10:49
8楼2012-05-11 10:40:39
feng198800
金虫 (正式写手)
- 应助: 15 (小学生)
- 金币: 860.2
- 散金: 234
- 红花: 1
- 帖子: 415
- 在线: 242.1小时
- 虫号: 700353
- 注册: 2009-02-12
- 性别: GG
- 专业: 化学反应工程
【答案】应助回帖
|
请问楼主活度系数：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 ] 是怎么回事？程序语句中没有体现你输入实验数据！ |
9楼2012-05-29 10:46:57