| 查看: 5530 | 回复: 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可以让我学习一下（我没搜到啊  ）。谢谢 |
回复此楼» 猜你喜欢
北京纳米能源与系统研究所王中林院士/曹南颖研究员课题组招收2026硕士研究生1名
已经有0人回复
-
山西大同大学物理学硕士研究生招收调剂生
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有270人回复
-
[调剂信息]211智能人工感知方向国家青年特聘专家课题组招收调剂研究生
已经有0人回复
-
［调剂信息］211智能人工感知方向国家青年特聘专家课题组招收调剂研究生
已经有0人回复
-
[调剂信息]211智能人工感知方向国家青年特聘专家课题组招收调剂研究生
已经有0人回复
-
一志愿物理所凝聚态物理352求调剂
已经有0人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有2人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有0人回复
-
070200求调剂，一志愿某211，288分
已经有17人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 动力学反应级数
已经有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人回复
» 抢金币啦！回帖就可以得到:
-
江西科技师范大学 材料与能源招收调剂研究生
+1/92
-
一志愿前4位代码0703均可报名-河北大学分析化学专业招收调剂考生
+2/92
-
西安石油大学新能源学院接收材料类、能源动力类、机械类、计算机类等专业专硕调剂生！
+1/89
-
2026 硕士研究生调剂｜河南工业大学化学化工学院名额开放
+2/58
-
2026年宁波工程学院招收硕士研究生（调剂）
+1/55
-
材料与化工专业有大量调剂指标
+1/52
-
湖北师范大学化学化工学院招收调剂
+1/36
-
河南工业大学化学化工学院-调剂
+1/34
-
吉大生科博一博士联培自荐-已发表+在投一作SCI3篇-踏实肯干能加班
+1/24
-
【轻松准入】【生物医药】2026青岛大学招博士生
+1/16
-
同济大学环境学院 肖倩特聘研究员课题组 招聘硕士/博士（长期有效）
+1/13
-
齐鲁工业大学 绿色造纸与资源循环全国重点实验室 招收调剂研究生
+1/11
-
青年长江学者团队 2026年 硕士研究生 招收调剂考生
+1/10
-
渤海大学储能材料课题组招收调剂研究生
+1/8
-
苏州大学医学部纳米生物医学方向招收调剂生1名
+1/6
-
中国科学院城市环境研究所与福建农林大学联合培养研究生调剂公告（第一批）
+1/6
-
滁州学院0860生物与医药专业调剂57人
+1/5
-
广东工业大学生态环境与海洋学院 杨盼教授团队招收硕士调剂生
+1/4
-
苏州科技大学胡老师课题组招2026年硕士生若干名
+1/2
-
Top-88悉尼科技大学数据科学/AI 招收2027年入学 校奖 博士生1到2名(国际和本地学生）
+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
回复
