24小时热门版块排行榜    

查看: 2645  |  回复: 5

lvst

新虫 (初入文坛)

[求助] matlab拟合反应速率常数 已有2人参与

小弟刚刚接触反应动力学,需要用matlab拟合动力学参数,很困惑,希望有高手指点,我的问题是下面的微分方程怎么解出来反应速率常数,每种物质的随时间变化的浓度都已经由实验测出。
A,B两种物质相互转换,k1为A向B转变的速率常数,k2为B向A转变的速率常数。x1,x2分别为A、B的浓度。
dx1/dt=-k1*x1+k2*x2
dx2/dt=-dx1/dt
附件为原始模型

matlab拟合反应速率常数
QQ截图.png
回复此楼

» 猜你喜欢

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

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

lsjc77

铜虫 (初入文坛)

如果你有一系列C和t的数据,而且足够表示浓度随时间的变化,可以直接简单的先做一个曲线拟合,这样你就得到了C和t的关系,然后回归k,不需要解微分方程,因为解微分方程还是要得到C和t的关系,但是你通过实验数据已经有了。
如果还想继续通过微分方程算:
计算这两个常微分方程,或者一个(x2 = 1-x1)代入可以弄掉第二个方程。
然后做参数的拟合,s = E(i)E(k){w(ki)*(c(ki) - ce(ki)^2} -> minimum, E在这里表示连加号,你可以用不一样的方程,反正意思就是让实验浓度和计算浓度的差达到最小。
但是在解微分方程之前你得给k0 和 E一个初值,这个初值最好接近终值,不要瞎猜,获得一个不错的初值的方式就是一开始说的方法,通过数据得到c和t的关系然后计算k.
不过动力学的话,从你这个体系来看,k1和k2是正逆,所以一般满足热力学约束,就是k1 = k2* Keq,Keq是平衡常数,如果你的时间足够,观测到体系浓度不变,那么可以通过那个平衡浓度计算keq, 对于微分方程,就是给一个足够的t,发现体系的浓度不变了,也得到了keq,然后检测你的k1和k2是否符合热力学关系。
这是一般情况,具体实现你自己琢磨吧
2楼2015-11-09 07:36:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lvst

新虫 (初入文坛)

引用回帖:
2楼: Originally posted by lsjc77 at 2015-11-09 07:36:02
如果你有一系列C和t的数据,而且足够表示浓度随时间的变化,可以直接简单的先做一个曲线拟合,这样你就得到了C和t的关系,然后回归k,不需要解微分方程,因为解微分方程还是要得到C和t的关系,但是你通过实验数据已 ...

非常感谢您的解答,请问第一步做曲线拟合的时候一般用什么模型拟合,多项式吗?回归k,具体的实现途径是什么,用matlab可以做吗、?由于本人是这方面的菜鸟,刚刚接触,可能一些基础的还不是很清楚,希望多解释一下
3楼2015-11-09 13:28:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

whyjackeyson

金虫 (小有名气)

【答案】应助回帖

感谢参与,应助指数 +1
引用回帖:
3楼: Originally posted by lvst at 2015-11-09 13:28:17
非常感谢您的解答,请问第一步做曲线拟合的时候一般用什么模型拟合,多项式吗?回归k,具体的实现途径是什么,用matlab可以做吗、?由于本人是这方面的菜鸟,刚刚接触,可能一些基础的还不是很清楚,希望多解释一下 ...

首先,要看你的三点分布如何,你的问题是有关反应动力学的,是不是还要考虑反应级数的因素。可以尝试几种非线性拟合模型,比较拟合优度。matlab完成这些工作没有问题,spss,origin也都可以。有数据的话,你可以再联系我

发自小木虫IOS客户端
4楼2015-11-09 18:21:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
微分方程拟合问题吧,1stOpt可以很容易求解的
5楼2015-11-09 20:53:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lsjc77

铜虫 (初入文坛)

引用回帖:
3楼: Originally posted by lvst at 2015-11-09 06:28:17
非常感谢您的解答,请问第一步做曲线拟合的时候一般用什么模型拟合,多项式吗?回归k,具体的实现途径是什么,用matlab可以做吗、?由于本人是这方面的菜鸟,刚刚接触,可能一些基础的还不是很清楚,希望多解释一下 ...

因为你的反应很简单,是一个一般的一级可逆反应,所以你应该可以找到一个解析解,
这个解析解是:
- ln { (C1 - C1e) / (C10 - C1e)} = { (R+1) / (R + X1e) }*k1 * t                 (1)
上边等式的左边变为:
-ln { 1 - X1/X1e}
X1 - 反应物摩尔分数
X1e - 反应物 平衡摩尔分率
R - R = X20 / X10 t=0 时的产物和反应物的比,如果X20 = 0 那 R= 0
所以实际上,你根据方程(1),画一条 -ln { 1 - X1/X1e}  和 t 的曲线,这个曲线应该是一条直线,你就得到k1了,
k2 = k1/Keq 得到k2
但也就到此为止了,如果想要得到 k = k0 exp( -e/RT) 中的 k0 和 e, 那你需要的是不同温度下的实验,重复以上步骤,得到不同温度下的 k1,k2
然后做一个图,ln{k} = ln{k0} - e/R * (1/T) 也就是ln{k}和 {1/T}的图,也是条直线,截距和斜率就是你想要的东西了。

所以事实上根据你的这个问题,完全没有数值计算的必要,当你遇到更加复杂的反应体系的时候再考虑数值方法,或者你需要首先把一般问题的解析解弄明白,当你遇到没法用解析解的时候再考虑数值方法。通常对于这种batch模型,不可逆反应一般都能找到一个解析解,可逆稍稍复杂一些,当反应级数比较奇怪,不是1或2,那基本上就只能用数值方法了。
由于基本上不接触有解析解的体系,我也没有立刻想到这个问题的适用方法,不过看到这个反应体系比较简单,就翻了翻书,你随便拿一本反应工程看看,就能找到了。所以对于这个问题没必要用matlab,用origin,或者excel就可以了。不过你可以试着学学怎么用数值方法解决这个问题,因为可以和解析解做对比,以后遇到复杂问题就好办了。
6楼2015-11-10 18:35:45
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 lvst 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 投票:  有多少人是今天查系统知道结果的? +16 爱看书的可乐 2026-08-26 18/900 2026-08-28 09:33 by winsaint
[基金申请] 哪位高人中了,把查询到的截图贴出来让我看看,让我长长见识 +5 yuleib84 2026-08-26 6/300 2026-08-28 00:02 by yudaoqian88
[教师之家] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +3 XLGfUIdOmEAO 2026-08-27 4/200 2026-08-27 22:26 by 3RAFyYXYSwGS
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +7 苏东坡二世 2026-08-23 7/350 2026-08-27 18:14 by 瞬息宇宙
[基金申请] 基金不中,共勉 +10 eulota 2026-08-26 10/500 2026-08-27 17:23 by lqllinqiaoli
[基金申请] 怎么看青基中了没有啊 +5 叶九微 2026-08-26 5/250 2026-08-27 10:35 by l_zh2008
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[硕博家园] 售SCI一区T0P文章,我:8.O.55.1.O.54,科目齐全,可+急 +3 LIbGuocjEEYw 2026-08-26 4/200 2026-08-27 02:01 by Ie9AyIAvGbvs
[基金申请] 范进中举一文的中心思想 +9 炎黄贵胄 2026-08-22 10/500 2026-08-26 15:40 by semaglutide
[基金申请] 2026国自然函评费到账 +22 羊腰板 2026-08-21 25/1250 2026-08-26 15:00 by zuocuiping
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件,有关系吗? +25 木水思豆 2026-08-25 28/1400 2026-08-26 14:53 by draco1987
[基金申请] 为什么国自然不能直接公布 +4 bjdxyxy 2026-08-26 4/200 2026-08-26 13:12 by qingmu1201
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
[基金申请] 明天应该可查了!? +6 chengyan1220 2026-08-23 6/300 2026-08-25 19:45 by zfd97
[基金申请] 如果此刻你正在为国基感到焦虑,不妨来听听这首《基金之外》 +8 scalable 2026-08-24 8/400 2026-08-25 12:52 by jnhyjjm
[基金申请] 建议基金发布提前给出明确的时间点 +13 kulium 2026-08-21 16/800 2026-08-24 16:27 by superceng
[基金申请] 看来今天不会放榜了? +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
信息提示
请填处理意见