24小时热门版块排行榜    

查看: 2323  |  回复: 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 的主题更新
信息提示
请填处理意见