| 查看: 5559 | 回复: 26 | |||
[交流]
Matlab同时拟合活化能，反应级数，速率常数
|
|||
|
各位，您好! 本人研究的方向是催化化学，最近在评估催化剂效果时想引入动力学的内容，对其中的一些地方有一些疑问。 受月之蓝版主的这个帖子http://muchong.com/bbs/viewthread.php?tid=6425538&target=self&page=1以及其它前辈的启发，对于反应级数以及反应速率常数的拟合有一定的了解。 但是如果我的动力学模型是这种形式：-dC/dt=k0*exp(-E/R*T)*C^m*P^n（只有一个反应物和氢气），其中P与T分别为压力与时间，就是说是一个加氢反应，试验中可以得到不同压力、不同温度下浓度随时间的变化值，应该如何拟合k0,E,m,n。 我的一个想法是先拟合不同压力下的k1（k1=k0*exp(-E/R*T)*P^n)（例如以上述帖子的代码方法），再根据k1=k2*P^n（k2=k0*exp(-E/R*T)）拟合出k2和n，此后拟合出不同温度下的k2，根据k2=k0*exp(-E/R*T)）拟合出E和k0 但是我觉得这种方法较为复杂，一个matlab里有好几个拟合程序，不知道是否有更简洁的CODE。 PS:不知道版里有没有相似的帖子有类似的CODE可以让我学习一下（我没搜到啊  ）。谢谢 |
回复此楼» 猜你喜欢
B区学生调剂-兰州交通大学材料科学与工程学院
已经有12人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有13人回复
-
物理学I论文润色/翻译怎么收费?
已经有83人回复
-
桂林理工大学物理学专业招收调剂
已经有18人回复
-
VASP 的一组 GPU / CPU 基准测试记录
已经有0人回复
-
津理工大学晶体材料全国重点实验室刘红军教授课题组招收博士生一名
已经有0人回复
-
【原创讨论】从电子约束到物质编辑：一套可迭代的环形磁场科技树
已经有0人回复
-
【方案分享】单环磁场+轴心控制+偏转导出电子束约束系统（可行性实验）
已经有6人回复
-
【修正版】单环用磁约束低速电子实验方案（简化版）
已经有0人回复
-
桂林理工大学物理学专业招收调剂，还有三个名额！！！
已经有19人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 动力学反应级数
已经有4人回复
- 求助用matlab三次样条法算一个反应的反应速率
已经有13人回复
- 求助用matlab拟合动力学方程
已经有9人回复
- 速率方程的问题求教！
已经有4人回复
- 有关反应速率常数的计算
已经有9人回复
- SCR反应动力学方程求解的问题
已经有4人回复
- 催化剂增大了速率常数但没有改变活化能-来自一篇经典文献，如何理解？
已经有3人回复
- 动力学参数拟合
已经有26人回复
- matlab 拟合反应动力学参数结果很差。大家帮忙看一下
已经有14人回复
- 根据活化能怎么判断一级反应还是零级反应
已经有12人回复
- 焓变 活化能 键离解能 能垒 反应速率 的巩固 求教
已经有14人回复
- 帮我看看聚合反应三个速率常数怎么算？
已经有9人回复
- 【求助】拟合动力学方程求助
已经有13人回复
- 【求助】催化反应动力学matlab计算各基元反应的速率常数时，该如何避免较小量被忽略？
已经有3人回复
- 【求助】反应速率常数怎么算
已经有5人回复
- 【求助】关于催化反应中的活化能
已经有24人回复
- 【求助】使用Matlab拟合反应动力学方程问题
已经有7人回复
- 【求助】使用Matlab预估动力学方程问题
已经有13人回复
» 抢金币啦！回帖就可以得到:
-
苏州国家实验室---南京大学联合培养博士招生（2026秋季入学）
+5/175
-
2026年手套箱应该怎么选？
+1/94
-
燕山大学亚稳材料全国重点实验室2026年博士研究生招生信息
+1/81
-
香港城市大学（东莞）雷进程博士课题组招聘研究助理
+1/77
-
万人领军
+3/47
-
南京航空航天大学 太阳能蒸发方向招收博士/科研助理
+1/39
-
厦门大学机械工程学科2026年招收博士后和科研助理
+1/38
-
郑州大学物理学院倪佩楠教授课题组诚邀博士后、青年教师及副教授加盟
+1/32
-
西工大国家级青年人才招聘微尺度传热方向博士后/研究助理
+1/29
-
北京航空航天大学-仿生界面材料科学全国重点实验室 康建新 教授招收博士一名
+2/18
-
广东公办省属高校招聘博士后2名，特别优秀可申请留校。年薪28-40万，科研奖励另计
+1/8
-
河南省医学科学院王宁利院士/苗彦彦研究员科研团队、伏羲慧眼云平台团队2026年博士、�
+1/7
-
贵州师范大学山地环境重点实验室调剂生态学5名，时间4月18日0点——4月19日10点
+1/7
-
全合成博士想找香港药学相关博后
+1/6
-
青岛科技大学高分子学院硕士招生
+1/5
-
湖北文理学院物理与电子工程学院材料与化工专业研究生调剂
+1/5
-
为什么很多发了优秀论文的博士却不再走科研之路？
+1/3
-
小身材大能量，易度质量流量计在雾化器测量中的精准控制
+1/2
-
文献求助
+1/1
-
中国药科大学招聘药物化学方向博士后（长期有效）
+1/1
2楼2014-08-25 17:31:36
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
小木虫: 金币+0.5, 给个红包，谢谢回帖
|
function k1k2k3k4 format long clear all clc tspan = [ 0.6 1.2 2 3 ]; C0 = [2 20 200]; %初值只有一个C0 k0 = [0 0 0 0]; %k0 E M n 参数的初值 lb = [0 0 0 0]; %上届 %ub = [100 10000 1000 1000 ]; %下界 dataTP1=... [ %t C XX XX XXX XX XXX X XX X ]; %实验数据1 dataTP2=... [ %t C XX XX XXX XX XXX X XX X ]; %实验数据2 dataTP3=... [ %t C XX XX XXX XX XXX X XX X ]; %实验数据3 [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,k0,lb,[],options,tspan,C0,dataTP1,dataTP2,dataTP3); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk0 = %.9f \n',k(1)) fprintf('\tE = %.9f \n',k(2)) fprintf('\tm = %.9f \n',K(3)) fprintf('\tn = %.9f \n',k(4)) fprintf(' The sum of the squares is: %.9e\n\n',resnorm) %----------------------------------------------------- function f = ObjFunc(k,tspan,C0,dataTP1,dataTP2,dataTP3) % 目标函数 T1=?;P1=?; [t XsimTP1] = ode23s(@KineticsEqs,tspan,C0(1),[],k,T1,P1); T2=? ;P2=?; [t XsimTP2] = ode23s(@KineticsEqs,tspan,C0(2),[],k,T2,P2); T3=?;P3=?; [t XsimTP3] = ode23s(@KineticsEqs,tspan,C0(3),[],k,T3,P3); f = [(XsimTP1(:,1)-dataTP1(:,2)) (XsimTP2(:,1)-dataTP2(:,2)) (XsimTP3(:,1)-dataTP3(:,2))]; %---------------------------------------------------------- function dCdt = KineticsEqs(t,C,k,T,P) % ODE模型方程 R=8; dC=-k(1)*exp(-k(2)/R*T)*(C^k(3))*(P^k(4)); %k(1)=k0,k(2)=E,k(3)=m;K(4)=n dCdt = [dC]; |
» 本帖已获得的红花（最新10朵）
heisekaishui18楼2014-08-27 11:33:16
得金币了，说声谢谢！5楼2014-08-25 17:54:40
9楼2014-08-25 18:17:25
10楼2014-08-25 21:13:01
11楼2014-08-25 23:33:17
12楼2014-08-26 00:14:30
14楼2014-08-26 21:40:00
15楼2014-08-27 10:58:07
16楼2014-08-27 11:00:36
17楼2014-08-27 11:02:13
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
小木虫: 金币+0.5, 给个红包，谢谢回帖
|
数据造了三组，对于看四个参数背景不了解，所以初值也选的不合理，以至于运行结果比较糟糕 |
» 本帖已获得的红花（最新10朵）
19楼2014-08-27 11:35:31
20楼2014-08-28 19:50:16
21楼2014-08-28 19:50:58
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
小木虫: 金币+0.5, 给个红包，谢谢回帖
|
function KineticsEst1_int11 clear all clc t=[0 20 40 100 270 450 750 1410 1590 1810];%时间 rA=[0.0004765 0.0004507 0.0004250 0.0003477 0.001549 0.0001117 0.0000396 0.0000018 0.0000038 0.0000062];%反应速率 Pa=[0.042729739 0.041609737 0.040801805 0.038526838 0.033412812 0.030519156 0.028695168 0.02792077 0.027918742 0.0278887981]; %乳酸的摩尔分率 Pb=[0.123670897 0.122550895 0.121742963 0.119467996 0.11435397 0.111460314 0.109636326 0.108861928 0.1088599 0.108829956];%正T醇的摩尔分率 Pc=[0.833599364 0.834719366 0.835527299 0.837802265 0.842916291 0.845809947 0.847633933 0.848408333 0.848410361 0.848440305];%水的摩尔分率 Pd=[0.001120002 0.001927934 0.004202901 0.009316927 0.012210583 0.014034571 0.014808969 0.014810996 0.014840941 0.014850941];%乳酸正丁酯的摩尔分率 %线性拟合 P=2.337132745*Pa.*Pb-0.564655417*Pc.*Pd;y=rA';X=[ones(size(y)) P']; b=X\y;k=b(2); %非线性拟合 beta0=[k] [beta,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,beta0,[],[],[],rA,Pa,Pb,Pc,Pd) %拟合效果图(实验与拟合的比较) figure(1);plot(t,rA,'.') r_poly=beta(1)*(2.337132745*Pa.*Pb-0.564655417*Pc.*Pd); hold on;plot(t, r_poly,'g') figure(2);plot(Pb,rA,'.'); hold on;plot(Pb, r_poly,'g') %-----------------------------------------------------------------------------一 Function f=ObjFunc(beta,rA,Pa,Pb,Pc,Pd) f=rA-beta(1)*(2.337132745*Pa.*Pb-0.564655417*Pc.*Pd); 我用的这个程序拟合反应速率常数为什么有问题的呢 |
22楼2014-09-05 17:20:40
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
小木虫: 金币+0.5, 给个红包，谢谢回帖
|
能否帮我指导下，乳酸与正丁醇酯化反应的反应级数和反应速率常数 数据如代码 function KineticsEst1_int11 % 动力学参数辨识: 用积分法进行反应速率分析得到速率常数k和反应级数n % Analysis of kinetic rate data by using the integral method % % Author: HUANG Huajiang % Copyright 2003 UNILAB Research Center, % East China University of Science and Technology, Shanghai, PRC % $Revision: 1.0 $ $Date: 2003/07/27 $ % % Reaction of the type -- rate = kCA^order % order - reaction order % rate -- reaction rate vector % CA -- concentration vector for reactant A % T -- vector of reaction time % N -- number of data points % k- reacion rate constant clear all clc t=[0 20 40 100 270 450 750 1410 1590 1810];%时间 rA=[0.0004765 0.0004507 0.0004250 0.0003477 0.001549 0.0001117 0.0000396 0.0000018 0.0000038 0.0000062];%反应速率 Pa=[0.042729739 0.041609737 0.040801805 0.038526838 0.033412812 0.030519156 0.028695168 0.02792077 0.027918742 0.0278887981]; %乳酸的摩尔分率 Pb=[0.123670897 0.122550895 0.121742963 0.119467996 0.11435397 0.111460314 0.109636326 0.108861928 0.1088599 0.108829956];%正T醇的摩尔分率 Pc=[0.833599364 0.834719366 0.835527299 0.837802265 0.842916291 0.845809947 0.847633933 0.848408333 0.848410361 0.848440305];%水的摩尔分率 Pd=[0.001120002 0.001927934 0.004202901 0.009316927 0.012210583 0.014034571 0.014808969 0.014810996 0.014840941 0.014850941];%乳酸正丁酯的摩尔分率 %线性拟合 P=2.337132745*Pa.*Pb-0.564655417*Pc.*Pd;y=rA';X=[ones(size(y)) P']; b=X\y;k=b(2); %非线性拟合 beta0=[k] [beta,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,beta0,[],[],[],rA,Pa,Pb,Pc,Pd) %拟合效果图(实验与拟合的比较) figure(1);plot(t,rA,'.') r_poly=beta(1)*(2.337132745*Pa.*Pb-0.564655417*Pc.*Pd); hold on;plot(t, r_poly,'g') figure(2);plot(Pb,rA,'.'); hold on;plot(Pb, r_poly,'g') %-----------------------------------------------------------------------------一 Function f=ObjFunc(beta,rA,Pa,Pb,Pc,Pd) f=rA-beta(1)*(2.337132745*Pa.*Pb-0.564655417*Pc.*Pd); 这个代码我照着书上弄的有错误 |
23楼2014-09-05 17:23:36
24楼2014-09-09 12:31:40
25楼2014-09-10 16:26:23
26楼2014-12-24 16:30:24
27楼2014-12-24 16:35:05
简单回复
2014-08-25 17:50
回复
heisekaishui(金币+5): 谢谢参与
2014-08-25 17:53
回复
heisekaishui(金币+5): 谢谢参与
dmbb6楼
2014-08-25 17:57
回复
heisekaishui(金币+5): 谢谢参与
xiejf7楼
2014-08-25 18:07
回复
heisekaishui(金币+5): 谢谢参与
1
@gemini@8楼
2014-08-25 18:12
回复
heisekaishui(金币+5): 谢谢参与
coal-xust13楼
2014-08-26 15:56
回复
