| 查看: 1840 | 回复: 13 | ||||
| 当前主题已经存档。 | ||||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||||
jlzeng木虫 (正式写手)
|
[交流]
【求助】反应动力学求解
|
|||
|
我在求反应动力学时遇到一个非常棘手的问题，请大家帮忙解决一下，先对每一个关注本贴的人表示感谢！ 问题描述： 有三个可逆反应组成的反应体系如下： TG+M→←DG+ME (正逆反应速率常数分别为k1、k2) DG+M→←MG+ME (正逆反应速率常数分别为k3、k4) MG+M→←GL+ME (正逆反应速率常数分别为k5、k6) 建立反应以上方程的数学模型如下： dTG/dt=-k1[TG][M]+ k2[DG][ME] dDG/dt=k1[TG][M]- k2[DG][ME] –k3[DG][M]+ k4[MG][ME] dMG/dt= k3[DG][M]- k4[MG][ME] –k5[MG][M]+ k6[GL][ME] dME/dt= k1[TG][M]- k2[DG][ME]+k3[DG][M]- k4[MG][ME] +k5[MG][M]- k6[GL][ME] 其中TG、DG、MG、ME的浓度是可测定的，各物质浓度间存在如下关系： [GL]=[TG]0-[TG]-[DG]-[MG]= 0.6380-[TG]-[DG]-[MG] ; [M]=[M]0-[ME]=4.0626-[ME]; 方程初始值为 [TG]0=0.6380 [DG]0=0 [MG]0=0 [GL]0=0 [ME]0=0 [M]0=4.0626 我用MATLAB求解某一条件下的动力学常数k，程序如下： clear all k0 = [1 1 3 1 5 1]; % 随意给定的参数初值 lb = [0 0 0 0 0 0]; % 随意给定的参数下限 ub = [1000 1000 1000 1000 1000 1000]; % 随意给定的参数上限 X0=[0.638;0;0;0]; %t=0时刻四种物质TG、DG、MG、ME的初始浓度 Xexp=... [0.6380 0 0 0 0.1512 0.1088 0.1244 1.0757 0.0955 0.0816 0.1454 1.2740 0.0719 0.0632 0.1431 1.3858 0.0500 0.0493 0.1400 1.4838 0.0451 0.0436 0.1283 1.5247 0.0361 0.0351 0.1226 1.5762 0.0346 0.0301 0.1028 1.6149 0.0261 0.0282 0.1016 1.6460 0.0237 0.0257 0.0977 1.6630 ]; % 25度0.8%催化剂条件下不同时间（0 0.5 1.5 2 3 4 5 7 10min）各浓度数据 % 第一种计算方法，使用函数fmincon()进行参数估计 [k,fval,flag] = fmincon(@jlzengObjFmincon,k0,[],[],[],[],lb,ub,[],[],X0,Xexp); fprintf('\n使用函数fmincon()估计得到的参数值为:\n') fprintf('\tk1 = %.4f\n',k(1)) fprintf('\tk2 = %.4f\n',k(2)) fprintf('\tk3 = %.4f\n',k(3)) fprintf('\tk4 = %.4f\n',k(4)) fprintf('\tk5 = %.4f\n',k(5)) fprintf('\tk6 = %.4f\n',k(6)) fprintf(' The sum of the squares is: %.1e\n\n',fval) k_fmincon = k; %第二种计算方法，使用函数Lsqnonlin()进行参数估计 [k,resnorm,residual,exitflag,Output,lambda,jacobian]=... lsqnonlin(@jlzengObjLsqnonlin,k0,lb,ub,[],X0,Xexp); ci=nlparci(k,residual,jacobian); fprintf('\n\n使用lsqnonlin()函数估计得到的参数值为：\n'),Output fprintf('\n使用函数lsqqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.4f\n',k(1)) fprintf('\tk2 = %.4f\n',k(2)) fprintf('\tk3 = %.4f\n',k(3)) fprintf('\tk4 = %.4f\n',k(4)) fprintf('\tk5 = %.4f\n',k(5)) fprintf('\tk6 = %.4f\n',k(6)) fprintf(' The sum of the squares is: %.1e\n\n',fval) k_fmincon = k; 三个m文件如下： function dXdt=jlzengkinetic(t,X,k) %建立微分方程组，供oed23s调用 dXdt=[ (-k(1)*X(1)*(4.0626-X(4))+k(2)*X(2)*X(4)) (k(1)*X(1)*(4.0626-X(4))-k(2)*X(2)*X(4)-k(3)*X(2)*(4.0626-X(4))+k(4)*X(3)*X(4)) (k(3)*X(2)*(4.0626-X(4))-k(4)*X(3)*X(4)-k(5)*X(3)*(4.0626-X(4))+k(6)*(0.638-X(1)-X(2)-X(3))*X(4)) (k(1)*X(1)*(4.0626-X(4))-k(2)*X(3)*(4.0626-X(4))+k(3)*X(2)*(4.0626-X(4))-k(4)*X(3)*X(4)+k(5)*X(3)*(4.0626-X(4))-k(6)*(0.638-X(1)-X(2)-X(3))*X(4))]; %-------------------------------------------------- function f = jlzengObjFmincon(k,X0,Xexp) %建立使用fmincon()进行参数估计的函数 tspan = [0 0.5 1 1.5 2 3 4 5 7 10]; [t X] = ode23s(@jlzengkinetic,tspan,X0,[],k); f = sum((X(:,1)-Xexp(:,1)).^2) + sum((X(:,2)-Xexp(:,2)).^2) ... + sum((X(:,3)-Xexp(:,3)).^2) + sum((X(:,4)-Xexp(:,4)).^2); %给物质浓度的计算值 %------------------------------------------------------ function f = jlzengObjLsqnonlin(k,X0,Xexp) %在前面已经用fmincon求得参数的基础上使用Lsqnonlin对参数进行更精确的计算 tspan =[0 0.5 1 1.5 2 3 4 5 7 10]; [t X] = ode23s(@jlzengkinetic,tspan,X0,[],k); f1 = X(:,1) - Xexp(:,1); f2 = X(:,2) - Xexp(:,2); f3 = X(:,3) - Xexp(:,3); f4 = X(:,4) - Xexp(:,4); f = [f1; f2; f3; f4]; %------------------------------------------------------- 运行结果发现1)两种方法求出的k值差别极大；2）不管将那组k值带入[t X] = ode23s(@jlzengkinetic,tspan,X0,[],k)求出的浓度与实际浓度差别很大；3）任意给定的k值的起始值对结果影响极大。 请高手帮忙！ 如果还有没表述清楚的，请提出来. [ Last edited by sunxiao on 2009-3-8 at 11:17 ] |
回复此楼» 收录本帖的淘帖专辑推荐
 matlab |
» 猜你喜欢
不自信的我
已经有11人回复
-
北核录用
已经有3人回复
-
要不要辞职读博？
已经有6人回复
-
实验室接单子
已经有3人回复
-
磺酰氟产物，毕不了业了！
已经有8人回复
-
求助:我三月中下旬出站，青基依托单位怎么办？
已经有10人回复
-
26申博（荧光探针方向，有机合成）
已经有4人回复
-
论文终于录用啦！满足毕业条件了
已经有26人回复
-
2026年机械制造与材料应用国际会议 （ICMMMA 2026）
已经有4人回复
-
Cas 72-43-5需要30g，定制合成，能接单的留言
已经有8人回复
luckyyjjun
铁虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 103
- 散金: 20
- 帖子: 547
- 在线: 233小时
- 虫号: 545943
- 注册: 2008-04-15
- 专业: 理论和计算化学
13楼2009-04-14 16:54:20