| 查看: 4565 | 回复: 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可以让我学习一下（我没搜到啊  ）。谢谢 |
回复此楼» 猜你喜欢
Error in routine lambda (100): wrong # or too many modes
已经有0人回复
-
祈福国自然面上
已经有61人回复
-
物理学I论文润色/翻译怎么收费?
已经有200人回复
-
请问castep计算出的能带图横轴数据缩放因子怎么计算
已经有1人回复
-
Result in Physics 从见刊开始 一般多久可以被SCI收录
已经有26人回复
-
面青地会评时间？？？
已经有8人回复
-
科研求助
已经有1人回复
-
深圳大学高等研究院暑期夏令营
已经有10人回复
-
深圳大学高等研究院暑期夏令营
已经有8人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 用MATLAB求解化学反应速率
已经有5人回复
- 动力学反应级数
已经有4人回复
- 求助用matlab三次样条法算一个反应的反应速率
已经有13人回复
- 求助用matlab拟合动力学方程
已经有9人回复
- 速率方程的问题求教！
已经有4人回复
- 有关反应速率常数的计算
已经有9人回复
- SCR反应动力学方程求解的问题
已经有4人回复
- 催化剂增大了速率常数但没有改变活化能-来自一篇经典文献，如何理解？
已经有3人回复
- 动力学参数拟合
已经有26人回复
- matlab 拟合反应动力学参数结果很差。大家帮忙看一下
已经有14人回复
- 化学动力学---------反应常数
已经有9人回复
- 根据活化能怎么判断一级反应还是零级反应
已经有12人回复
- 焓变 活化能 键离解能 能垒 反应速率 的巩固 求教
已经有14人回复
- 帮我看看聚合反应三个速率常数怎么算？
已经有9人回复
- 【求助】拟合动力学方程求助
已经有13人回复
- 【求助】催化反应动力学matlab计算各基元反应的速率常数时，该如何避免较小量被忽略？
已经有3人回复
- 【求助】反应动力学拟合
已经有16人回复
- 【求助】反应速率常数怎么算
已经有5人回复
- 【求助】催化反应动力学
已经有9人回复
- 【求助】关于催化反应中的活化能
已经有24人回复
- 【求助】使用Matlab拟合反应动力学方程问题
已经有7人回复
- 【求助】使用Matlab预估动力学方程问题
已经有13人回复
» 抢金币啦！回帖就可以得到:
-
年薪最高36万，南京大学诚聘天文、行星、地球和大气科学领域研究系列岗位（含博士后）
+1/473
-
招保研硕士-化学材料专业-北京交通大学2024年全国优秀大学生暑期夏令营通知
+1/185
-
计算机土博工作建议交流
+1/90
-
内蒙古大学-武利民教授-赵岩研究员-张江威研究员联合团队招聘高层次青年人才及优秀博�
+1/89
-
上海交大储能方向高薪招聘博士后研究员
+1/50
-
诚征女友坐标郑州非诚勿扰
+1/49
-
中科院大连化物所先进生物制造研究中心诚聘生物/化学科研助理博士后
+2/48
-
招聘制备纯化研究员
+1/44
-
分析化学专业教师招聘（硕士）
+1/18
-
复旦大学化学系凡勇教授团队招聘博士后
+1/18
-
【招聘】四川大学高分子科学与工程学院-仿生医用高分子课题组招收科研助理/硕士/博士
+1/16
-
西北工业大学 力学与土木建筑学院 断裂力学计算仿真方向 招收硕士研究生
+1/9
-
强磁场科学中心马文彬教授课题组招聘博士后
+1/6
-
一家专门做高纯度化学品提纯的公司，有前途吗？
+1/6
-
机器学习方向博士后（三位）招聘-深圳大学及KAUST深圳创新中心
+1/4
-
Sustainable Structures 第二届青年编委招募启事（正式通知）
+1/4
-
北欧-东芬兰大学招博士生一名-月工资2600欧以上-7月12号截止
+1/3
-
[18/6更新] 港教大招生 | 无线通信 | 博后Postdoc | 博士生PhD (全奖,可直博)
+1/2
-
深圳大学高等研究院暑期夏令营
+1/2
-
聚氨酯水凝胶
+1/1
+1/473
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
回复
