| 查看: 1918 | 回复: 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 |
» 猜你喜欢
一志愿西电085401数一英一299求调剂 六级521
已经有3人回复
-
生物考研337分求调剂
已经有5人回复
-
求调剂：085600材料与化工，考材科基，总分319
已经有14人回复
-
285求调剂
已经有10人回复
-
化学工程085602 305分求调剂
已经有29人回复
-
理学07化学 303求调剂
已经有5人回复
-
081200-11408-276学硕求调剂
已经有3人回复
-
土木304求调剂
已经有5人回复
-
一志愿中国科学院大学265求调剂
已经有7人回复
-
286求调剂
已经有8人回复
hitzhang
木虫 (正式写手)
- 仿真EPI: 1
- 应助: 0 (幼儿园)
- 贵宾: 2.15
- 金币: 1376.7
- 散金: 969
- 红花: 8
- 帖子: 862
- 在线: 226.4小时
- 虫号: 390575
- 注册: 2007-06-02
- 性别: GG
- 专业: 无机非金属类电介质与电解
★ ★ ★ ★ ★ ★ ★ ★
sunxiao(金币+8,VIP+0):重奖，欢迎来做我们的专家顾问 3-6 23:58
kuhailangyu(金币+0,VIP+0):你参与的主题已有新的进展，敬请关注，谢谢！ 3-14 15:57
sunxiao(金币+8,VIP+0):重奖，欢迎来做我们的专家顾问 3-6 23:58
kuhailangyu(金币+0,VIP+0):你参与的主题已有新的进展，敬请关注，谢谢！ 3-14 15:57
|
微分方程组： 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] 有4个方程，可简写为：c'=a*k。你测定了10个时刻4个物种的浓度，可计算出每个时刻的矩阵a，另外这10个时刻各个物种的反应速率c'已估计出。在每个时刻上满足c'=a*k，有4个方程，那么十个时刻共有40个方程。可把这40个方程简写为：C'=A*K。数学上称这类方程为超定方程，即方程组中方程的个数多于未知量的个数，一般情况超定方程无解，只能找到一个最近似的解，主要应用最小二乘法，使方程左右两边尽量相等。具体算法可参阅相关文献。 在matlab中解这个方程显得异常简单，K=A\C'。 一个严重的问题：结果代入微分方程模型进行述职求解时和实验偏差较大。我想一个很重要的原因是我把矩阵A算错了，因为繁琐避免不了马虎，如果矩阵A没有错误，我想最可能的原因是微分方程模型不够合理。仅供参考。 |
7楼2009-03-06 19:36:58
jlzeng
木虫 (正式写手)
- 应助: 10 (幼儿园)
- 金币: 2413.7
- 散金: 5
- 红花: 4
- 帖子: 471
- 在线: 214小时
- 虫号: 427388
- 注册: 2007-07-29
- 性别: GG
- 专业: 能源化工
2楼2009-03-04 10:45:37
jlzeng
木虫 (正式写手)
- 应助: 10 (幼儿园)
- 金币: 2413.7
- 散金: 5
- 红花: 4
- 帖子: 471
- 在线: 214小时
- 虫号: 427388
- 注册: 2007-07-29
- 性别: GG
- 专业: 能源化工
~~~~ 3楼2009-03-05 09:10:15
hitzhang
木虫 (正式写手)
- 仿真EPI: 1
- 应助: 0 (幼儿园)
- 贵宾: 2.15
- 金币: 1376.7
- 散金: 969
- 红花: 8
- 帖子: 862
- 在线: 226.4小时
- 虫号: 390575
- 注册: 2007-06-02
- 性别: GG
- 专业: 无机非金属类电介质与电解
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
sunxiao(金币+3,VIP+0):欢迎交流，得到楼主认可后，将继续追加奖励 3-5 23:33
jlzeng(金币+10,VIP+0):很有用，但还没有完全解决，非常感谢！ 3-6 15:48
sunxiao(金币+3,VIP+0):欢迎交流，得到楼主认可后，将继续追加奖励 3-5 23:33
jlzeng(金币+10,VIP+0):很有用，但还没有完全解决，非常感谢！ 3-6 15:48
|
我看了一下你的方法，感觉很新颖，是一种纯粹的数值解法，想必你的MATLAB功底不浅。你主要是想通过两个优化工具箱函数利用最小二乘法解微分方程组进而求出速率常数。但是效果不好，我想你下面的分析是有道理的，就算法有效性而言fmincon这个函数口碑不是很好，多数都是用LINGO软件求优化问题，更别说求微分方程组了，而且透明度也不高好像在解决一个黑箱问题，求解对初始条件很敏感。 我这里有个想法仅供你参考： 这个问题的最终目的是求最佳速率常数以尽量满足微分方程组，现在，你有了各个物种随时间变化的浓度，如果知道浓度随时间的变化速率，就可以直接代入微分方程组应用最小而乘求出速率常数k。 浓度变化速率的最佳估计：可以根据你的实验数据构造出在实验范围内符合较好的各个物种的浓度方程。以下是拟和结果： /////////////////////////////////////////////////////////////////////////////// TG: f(x) = a*exp(b*x) + c*exp(d*x) Coefficients (with 95% confidence bounds): a = 0.5502 (0.5149, 0.5856) b = -3.983 (-4.797, -3.17) c = 0.08751 (0.05922, 0.1158) d = -0.1857 (-0.2901, -0.08124) Goodness of fit: SSE: 0.0004627 R-square: 0.9985 Adjusted R-square: 0.9978 RMSE: 0.008782 DG: f(x) = (p1*x^2 + p2*x + p3) / (x^2 + q1*x + q2) Coefficients (with 95% confidence bounds): p1 = 0.01908 (0.01341, 0.02475) p2 = 0.05717 (0.01411, 0.1002) p3 = -1.43e-007 (-0.0008759, 0.0008756) q1 = -0.2444 (-0.8719, 0.3831) q2 = 0.1788 (0.04132, 0.3162) Goodness of fit: SSE: 1.816e-005 R-square: 0.9979 Adjusted R-square: 0.9963 RMSE: 0.001906 MG: f(x) = (p1*x^2 + p2*x + p3) / (x^2 + q1*x + q2) Coefficients (with 95% confidence bounds): p1 = 0.07364 (0.04582, 0.1015) p2 = 0.3095 (-0.1137, 0.7328) p3 = 1.802e-005 (-0.00794, 0.007976) q1 = 0.9642 (-1.498, 3.426) q2 = 0.6633 (0.1897, 1.137) Goodness of fit: SSE: 0.0001089 R-square: 0.9933 Adjusted R-square: 0.988 RMSE: 0.004666 ME: f(x) = (p1*x + p2) / (x^2 + q1*x + q2) Coefficients (with 95% confidence bounds): p1 = 2169 (-1.654e+004, 2.088e+004) p2 = 0.9198 (-23.22, 25.06) q1 = 1265 (-9708, 1.224e+004) q2 = 406.4 (-3050, 3862) Goodness of fit: SSE: 0.003249 R-square: 0.9986 Adjusted R-square: 0.9978 RMSE: 0.02327 /////////////////////////////////////////////////////////////////////// 这几个浓度方程的拟和结果如如附件中所示。 然后就可以计算在实验时间点上TG, DG, MG, ME的浓度对时间的导数了，结果如下：//////////////////////////////////////////////////////////////////////// >> dTGdt=[-2.20812330989885;-0.313903194394831;-0.0543072642775238;-0.0178664548644437;-0.0119669262873315;-0.00932216291313879;-0.00773110042346866;-0.00642090783831632;-0.00442930197092527;-0.00253772213004114]' dTGdt = -2.2081 -0.3139 -0.0543 -0.0179 -0.0120 -0.0093 -0.0077 -0.0064 -0.0044 -0.0025 >> dDGdt=[0.319813865326457;-0.0194568056055164;-0.0513088546928754;-0.0279068843505723;-0.0164124937537988;-0.00737171342598024;-0.00411902140843388;-0.00261612659035434;-0.00131890646666286;-0.000639095343223927;]' dDGdt = 0.3198 -0.0195 -0.0513 -0.0279 -0.0164 -0.0074 -0.0041 -0.0026 -0.0013 -0.0006 >> dMGdt=[0.466367050839715;0.0999072143408668;0.00933499024234685;-0.00972884650823934;-0.0127379996363272;-0.0104546020606972;-0.00764618068387975;-0.00566858666344666;-0.00339187374122912;-0.00186232961806277;]' dMGdt = 0.4664 0.0999 0.0093 -0.0097 -0.0127 -0.0105 -0.0076 -0.0057 -0.0034 -0.0019 >> dMEdt=[5.28500993498458;0.818198974835815;0.315950868952596;0.165592835535058;0.101261692838514;0.0484953510802743;0.0278523652553689;0.0177021334686043;0.00842980322026890;0.00327667396220621;]' dMEdt = 5.2850 0.8182 0.3160 0.1656 0.1013 0.0485 0.0279 0.0177 0.0084 0.0033 >> //////////////////////////////////////////////////////////////////////////// 之后的工作就是把这四组数据和实验数据代入微分方程组，求出最佳速率常数。 |
4楼2009-03-05 21:16:14
