版块导航: 正在加载中...

客户端APP下载

登录注册

应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜

heisekaishui

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 215
帖子: 52
在线: 59.2小时
虫号: 2445260

[交流] 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可以让我学习一下（我没搜到啊

）。谢谢

» 猜你喜欢

求国际会议网站已经有1人回复
求取一些关于纳米材料和纳米技术相关的英文PPT。已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有97人回复
【复旦大学】二维材料方向招收2026年博士研究生1名已经有0人回复
北京纳米能源与系统研究所王中林院士/曹南颖研究员课题组2026级硕/博/博后招生已经有10人回复
荷兰Utrecht University超快太赫兹光谱王海教授课题招收2026 CSC博士生已经有22人回复
反铁磁体中的磁性切换：两种不同的机制已成功可视化已经有0人回复
26申博推荐：南京航空航天大学国际前沿科学研究院光学方向招收博士生！已经有0人回复
求标准粉末衍射卡号 ICDD 01-076-1802 已经有0人回复

» 本主题相关价值贴推荐，对您同样有帮助:

动力学反应级数已经有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楼2014-08-25 17:30:13

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

回帖支持 ( 显示支持度最高的前 50 名 )

heisekaishui

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 215
帖子: 52
在线: 59.2小时
虫号: 2445260

纯粹的初学者，如果叙述有什么问题希望大家指出。

赞一下(1人)

2楼2014-08-25 17:31:36

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

whqs8426212

铜虫 (正式写手)

应助: 5 (幼儿园)
金币: 13.2
帖子: 782
在线: 685.7小时
虫号: 2572504

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

引用回帖:

16楼: Originally posted by heisekaishui at 2014-08-27 11:00:36
非常感谢，对，数据格式是这样的。
主要是数据还没有做出来，如果可以的话可以让我看一下程序吗，我想看一下大概的思路。...

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];

赞一下(1人)

» 本帖已获得的红花（最新10朵）

18楼2014-08-27 11:33:16

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

tangjiguo

木虫 (职业作家)

应助: 0 (幼儿园)
金币: 1379
帖子: 3705
在线: 370.2小时
虫号: 2353662

★
heisekaishui(金币+5): 谢谢参与

得金币了，说声谢谢！

5楼2014-08-25 17:54:40

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 215
帖子: 52
在线: 59.2小时
虫号: 2445260

同学们金币诚可贵，回答更重要啊。

9楼2014-08-25 18:17:25

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

计算强帖: 4
应助: 1641 (讲师)
金币: 15037.3
帖子: 3410
在线: 1223.5小时
虫号: 291104

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

用1stOpt应该可以很好解决啊。

10楼2014-08-25 21:13:01

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

whqs8426212

铜虫 (正式写手)

应助: 5 (幼儿园)
金币: 13.2
帖子: 782
在线: 685.7小时
虫号: 2572504

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

程序大致写好，楼主可以将试验数据按楼下格式贴出，我们运行看看

11楼2014-08-25 23:33:17

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

whqs8426212

铜虫 (正式写手)

应助: 5 (幼儿园)
金币: 13.2
帖子: 782
在线: 685.7小时
虫号: 2572504

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

楼主方不方便给我一些数据，如：

某T，P（应该是多个吧？）
t（/s） C(mol)
0          XX
1          XX
2          XX

12楼2014-08-26 00:14:30

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

程源洪

铁虫 (初入文坛)

应助: 0 (幼儿园)
金币: 32.4
帖子: 21
在线: 5.5小时
虫号: 1282324

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

直接利用origin 的非线性拟合

14楼2014-08-26 21:40:00

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 215
帖子: 52
在线: 59.2小时
虫号: 2445260

引用回帖:

10楼: Originally posted by dingd at 2014-08-25 21:13:01
用1stOpt应该可以很好解决啊。

主要是实验室没有买1stOpt，所以想用Matlab做一下。

15楼2014-08-27 10:58:07

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 215
帖子: 52
在线: 59.2小时
虫号: 2445260

引用回帖:

12楼: Originally posted by whqs8426212 at 2014-08-26 00:14:30
楼主方不方便给我一些数据，如：

某T，P（应该是多个吧？）
t（/s） C(mol)
0          XX
1          XX
2          XX

非常感谢，对，数据格式是这样的。
主要是数据还没有做出来，如果可以的话可以让我看一下程序吗，我想看一下大概的思路。

16楼2014-08-27 11:00:36

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 215
帖子: 52
在线: 59.2小时
虫号: 2445260

引用回帖:

14楼: Originally posted by 程源洪 at 2014-08-26 21:40:00
直接利用origin 的非线性拟合

谢谢回复，如果用Origin的话是分开拟合的吧？

17楼2014-08-27 11:02:13

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

whqs8426212

铜虫 (正式写手)

应助: 5 (幼儿园)
金币: 13.2
帖子: 782
在线: 685.7小时
虫号: 2572504

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

引用回帖:

16楼: Originally posted by heisekaishui at 2014-08-27 11:00:36
非常感谢，对，数据格式是这样的。
主要是数据还没有做出来，如果可以的话可以让我看一下程序吗，我想看一下大概的思路。...

数据造了三组，对于看四个参数背景不了解，所以初值也选的不合理，以至于运行结果比较糟糕

» 本帖已获得的红花（最新10朵）

19楼2014-08-27 11:35:31

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 215
帖子: 52
在线: 59.2小时
虫号: 2445260

送红花一朵

引用回帖:

19楼: Originally posted by whqs8426212 at 2014-08-27 11:35:31
数据造了三组，对于看四个参数背景不了解，所以初值也选的不合理，以至于运行结果比较糟糕...

非常感谢，其实我理解没错的话，其实Matlab解此类问题的方法就是同时拟合不同条件的下的多组数据，得到最优解？
以及我还请问程序里lsqnonlin函数中的数据dataTP1等，我看其它程序中有的会设置yexp = dataTP1(:,2);然后lsqnonlin中的数据使用不包含时间t的数据 lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,C0,yexp);这样的，我想请问这和您直接使用dataTP1有什么区别，谢谢。

20楼2014-08-28 19:50:16

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 215
帖子: 52
在线: 59.2小时
虫号: 2445260

送红花一朵

引用回帖:

18楼: Originally posted by whqs8426212 at 2014-08-27 11:33:16
function k1k2k3k4
format long
clear all
clc
tspan = ;
C0 = ; %初值只有一个C0
k0 = ; %k0 E M n 参数的初值
lb = ; %上届
%ub = ; %下界

dataTP1=...
; %实验数据1
dataTP2=.. ...

非常感谢！

21楼2014-08-28 19:50:58

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

guanwei

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 1075.1
帖子: 247
在线: 53.7小时
虫号: 1638561

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

引用回帖:

12楼: Originally posted by whqs8426212 at 2014-08-26 00:14:30
楼主方不方便给我一些数据，如：

某T，P（应该是多个吧？）
t（/s） C(mol)
0          XX
1          XX
2          XX

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

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

guanwei

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 1075.1
帖子: 247
在线: 53.7小时
虫号: 1638561

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

引用回帖:

18楼: Originally posted by whqs8426212 at 2014-08-27 11:33:16
function k1k2k3k4
format long
clear all
clc
tspan = ;
C0 = ; %初值只有一个C0
k0 = ; %k0 E M n 参数的初值
lb = ; %上届
%ub = ; %下界

dataTP1=...
; %实验数据1
dataTP2=.. ...

能否帮我指导下，乳酸与正丁醇酯化反应的反应级数和反应速率常数
数据如代码

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

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

whqs8426212

铜虫 (正式写手)

应助: 5 (幼儿园)
金币: 13.2
帖子: 782
在线: 685.7小时
虫号: 2572504

引用回帖:

20楼: Originally posted by heisekaishui at 2014-08-28 19:50:16
非常感谢，其实我理解没错的话，其实Matlab解此类问题的方法就是同时拟合不同条件的下的多组数据，得到最优解？
以及我还请问程序里lsqnonlin函数中的数据dataTP1等，我看其它程序中有的会设置yexp = dataTP1(:,2 ...

没区别，一样的

[ 发自小木虫客户端 ]

24楼2014-09-09 12:31:40

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 215
帖子: 52
在线: 59.2小时
虫号: 2445260

引用回帖:

24楼: Originally posted by whqs8426212 at 2014-09-09 12:31:40
没区别，一样的
...

谢谢，非常感谢。

25楼2014-09-10 16:26:23

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

fusheng1990

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 1
帖子: 3
在线: 4.6小时
虫号: 3347068

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

还有人关注这个帖吗？我的数据里产物种类比较多，每个温度下都有好几个数据
比如
t c1 c2 c3
x xx xx xx
x xx xx xx
x xx xx xx
这个怎么模拟啊？

26楼2014-12-24 16:30:24

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

fusheng1990

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 1
帖子: 3
在线: 4.6小时
虫号: 3347068

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

引用回帖:

25楼: Originally posted by heisekaishui at 2014-09-10 16:26:23
谢谢，非常感谢。...

我的数据里产物种类比较多
比如
t时间 c1 c2 c3表示各产物的浓度，每个温度得到类似的数据
t c1 c2 c3
x xx xx xx
x xx xx xx
x xx xx xx
这个怎么模拟啊？

27楼2014-12-24 16:35:05

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

简单回复

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

2014-08-25 18:12 回复

heisekaishui(金币+5): 谢谢参与

2014-08-26 15:56 回复

相关版块跳转我要订阅楼主 heisekaishui 的主题更新