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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » Matlab同时拟合活化能，反应级数，速率常数
271/1返回列表
查看: 6099  |  回复: 26
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

heisekaishui

银虫 (小有名气)

[交流] 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楼 2014-08-25 17:30:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

heisekaishui

银虫 (小有名气)
纯粹的初学者，如果叙述有什么问题希望大家指出。
赞一下(1人) 回复此楼
2楼2014-08-25 17:31:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

whqs8426212

铜虫 (正式写手)
★
小木虫: 金币+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朵）

heisekaishui
18楼2014-08-27 11:33:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

tangjiguo

木虫 (职业作家)
★
heisekaishui(金币+5): 谢谢参与
得金币了，说声谢谢！
赞一下 回复此楼
5楼2014-08-25 17:54:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)
同学们金币诚可贵，回答更重要啊。
赞一下 回复此楼
9楼2014-08-25 18:17:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
用1stOpt应该可以很好解决啊。
赞一下 回复此楼
10楼2014-08-25 21:13:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

whqs8426212

铜虫 (正式写手)
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
程序大致写好，楼主可以将试验数据按楼下格式贴出，我们运行看看
赞一下 回复此楼
11楼2014-08-25 23:33:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

whqs8426212

铜虫 (正式写手)
★
小木虫: 金币+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.5, 给个红包，谢谢回帖
直接利用origin 的非线性拟合
赞一下 回复此楼
14楼2014-08-26 21:40:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)
引用回帖:
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

银虫 (小有名气)
引用回帖:
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

银虫 (小有名气)
引用回帖:
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

铜虫 (正式写手)
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
引用回帖:
16楼: Originally posted by heisekaishui at 2014-08-27 11:00:36
非常感谢，对，数据格式是这样的。
主要是数据还没有做出来，如果可以的话可以让我看一下程序吗，我想看一下大概的思路。...

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

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

heisekaishui
19楼2014-08-27 11:35:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

heisekaishui

银虫 (小有名气)
送红花一朵
引用回帖:
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

银虫 (小有名气)
送红花一朵
引用回帖:
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.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.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

铜虫 (正式写手)
引用回帖:
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

银虫 (小有名气)
引用回帖:
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.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.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的回帖
简单回复
4993972323楼
2014-08-25 17:50   回复  
heisekaishui(金币+5): 谢谢参与
aixingang4楼
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   回复  
相关版块跳转 我要订阅楼主 heisekaishui 的主题更新
271/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭