版块导航
正在加载中...
客户端APP下载
论文辅导
登录
注册
帖子
帖子
用户
本版
应《网络安全法》要求,自2017年10月1日起,未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用,请尽快对帐号进行手机号验证,感谢您的理解与支持!
24小时热门版块排行榜
>
论坛更新日志
(5402)
>
虫友互识
(938)
>
休闲灌水
(224)
>
招聘信息布告栏
(174)
>
文献求助
(160)
>
论文投稿
(127)
>
基金申请
(102)
>
找工作
(94)
>
导师招生
(93)
>
硕博家园
(65)
>
职场人生
(64)
>
考博
(61)
>
博后之家
(36)
>
绿色求助(高悬赏)
(27)
>
育儿交流
(21)
>
留学DIY
(14)
小木虫论坛-学术科研互动平台
»
计算模拟区
»
计算模拟
»
1stOpt编写常微分方程组计算反应速率问题
5
1/1
返回列表
查看: 1109 | 回复: 11
只看楼主
@他人
存档
新回复提醒
(忽略)
收藏
在APP中查看
【悬赏金币】回答本帖问题,作者ykai_1991将赠送您 100 个金币
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖
ykai_1991
金虫
(初入文坛)
应助: 0
(幼儿园)
金币: 967
帖子: 36
在线: 19.8小时
虫号: 2054876
注册: 2012-10-11
专业: 化学反应工程
[
求助
]
1stOpt编写常微分方程组计算反应速率问题
已有2人参与
还请各位大神半忙查看下,我也是来学习的哈。另外烦请dingd帮忙查看下代码是否合理并计算下结果。
反应过程和微分方程见附图。
在编写中,数字代表不同的链长,物质是随着反应增长的。
比较麻烦的是微分方程中的多组分的加和,并且每个物质都有个方程,都写出来比较费时间,不知道有没有简便点的方法。
CODE:
本人这么写代码:([b]需要求每步反应速率k值(即为ka和kc),参数很多,不知能否计算出来;并且提供三组数据[/b])
Parameters ka(10), kc1(10), kc2(10), kc3(10), kc4(10), kc5(10), kc6(10), kc7(10), kc8(10), kc9(10), kc10(10) ; [b]这里的kc有方法简化下么?不然就要写10个了[/b]
Variable t, P(10), Q(10), F; [b]P1~P10也可以这么定义变量么,直接用P(10)?[/b]
Plot t[x], P1, P2, Q1, Q2, F; [b]不同组分的含量不一样,在一张图中会画出来看不出来?如果不同物质分开作图,怎么写?[/b]
Plot t[x], P3, P4, Q3, Q4;
ODEFunction
P1'=-ka1*P1*F-P1*Sum(j=1:10)(kc1[j]*Q[j]);
P2'=-ka2*P2*F-P2*Sum(j=1:10)(kc2[j]*Q[j])+Sum(j=1:(2-1))(kc2[j]*P2*Q[2-j]); [b]式子中的数字2都和P2的2有关系,如果P2变成P3,那么其他的2也都为3,那么针对P3'到P10',有没有更为方便的方法?下面再写一个P3'的,要一直写到P10'的
P3'=-ka3*P3*F-P3*Sum(j=1:10)(kc3[j]*Q[j])+Sum(j=1:(3-1))(kc3[j]*P3*Q[3-j]);[/b]
同样类似递推编写的方法也适用于Q1'到Q10',因为都只是角标不一样,加和的物质产生变化而已,在matlab中就可以定义出来比较方便,针对Q1'如下:
Q1'=-Q1*Sum(j=1:10)(kc1[j]*P[j])+ka1*P1*F; [b]同理,Q2'=-Q2*Sum(j=1:10)(kc2[j]*P[j])+ka2*P2*F; 也要一直写到Q10'[/b]
F'=-F*Sum(j=1:10)(ka[j]*P[j]);
Data t, P1, Q1, P2;
10, 30, 40, 50, 60, 75, 110, 150, 230;
10, 30, 40, 50, 60, 75, 110, 150, 230;
6.288097876, 4.91765691, 3.878352807, 3.298605506, 2.850889397, 2.241683203, 1.781706784, 0.910037534, 0.693987349;
0.257670552, 0.312407281, 0.244993608, 0.189709656, 0.133296065, 0.0803651, 0.033400217, 0.010654615, 0.005606342;
0.236493005, 0.607274759, 0.615323542, 0.623870854, 0.61234489, 0.562467845, 0.538393919, 0.347783178, 0.305460229;
Data t, P1, Q1, P2;
10, 30, 50, 60, 70, 85, 120, 170, 230, 305;
5.882896195, 4.593250328, 3.701714247, 2.910582239, 2.571710069, 2.416218035, 1.920578469, 1.68785903, 1.416030469, 1.287574203;
0.255794223, 0.297297319, 0.219253251, 0.161355927, 0.121348782, 0.081854014, 0.039094897, 0.016390052, 0.008676848, 0.005710866;
0.203988176, 0.534247028, 0.637609417, 0.634725031, 0.638037303, 0.569498973, 0.479824153, 0.477902906, 0.424218287, 0.411113563;
Data t, P1, Q1, P2;
10, 30, 50, 60, 70, 85, 130, 190, 237;
5.879605129, 4.718556629, 3.450891407, 3.045516702, 2.681262558, 2.428241655, 1.569209161, 1.234288027, 1.237674719;
0.22471112, 0.268468006, 0.177155403, 0.1332699, 0.094942111, 0.071844393, 0.019781466, 0.008520798, 0.00503865;
0.226510416, 0.569687824, 0.596106833, 0.599362637, 0.580304305, 0.666507432, 0.464066174, 0.389844442, 0.423282218;
最后再附上编写过的微分方程组部分的matlab程序,便于您对照。
function dydt=KineticsEqs(t, y, ka, kc)
m = 10;
n = 10;
for i=1:m %-----------------------将P赋值给y-----------------------------
P(i) = y(i);
end
for j=1:n %-----------------------将Q赋值给y-----------------------------
Q(j) = y(m+j);
end
F = y(m+n+1);%-----------------------将F赋值给y-----------------------------
%----------------------------以下为动力学微分方程----------------------------
for j=1:n
polynomial_1(j)=kc(1,j)*Q(j);
end
r(1) = -( ka(1)*P(1)*F+P(1)*sum(polynomial_1) ); %-------微分方程1--------
for i=1:m
for j=1:n
polynomial_2(j)=kc(i,j)*Q(j);
end
polysum_2(i)=sum(polynomial_2);
end
for i=2:m
for j=1:i-1
polynomial_3(j)=kc(i,j)*P(j)*Q(i-j);
end
polysum_3(i)=sum(polynomial_3);
end
r(i) = -( ka(i)*P(i)*F+P(i)*polysum_2(i)-polysum_3(i) );
for i=1:m
for j=1:n
polynomial_4(j)=kc(i,j)*P(j);
end
r(i+n) = -( Q(i)*sum(polynomial_4)-ka(i)*P(i)*F );%-----微分方程3-----
end
for j=1:n
polynomial_5(j)=ka(j)*P(j);
end
r(m+n+1) = -( F*sum(polynomial_5) ); %-------微分方程4--------
dydt = r';
反应方程
微分方程组
[
Last edited by dbb627 on 2015-1-15 at 23:45
]
回复此楼
» 猜你喜欢
求助将cif文件转化为jade
已经有1人回复
金属氧化物薄膜空气退火和真空退火对表面粗糙度的影响
已经有0人回复
物理学I论文润色/翻译怎么收费?
已经有119人回复
Apex2单晶测试
已经有0人回复
结构建模
已经有0人回复
声子谱计算问题
已经有2人回复
已结束
已经有23人回复
Error in routine lambda (100): wrong # or too many modes
已经有0人回复
祈福国自然面上
已经有55人回复
请问castep计算出的能带图横轴数据缩放因子怎么计算
已经有1人回复
» 本主题相关价值贴推荐,对您同样有帮助:
关于1stOpt程序问题
已经有5人回复
1stopt-- 求给定范围初始值,求满足方程组的初始值的最优解
已经有3人回复
请教高手,如何解非线性方程组!!!
已经有13人回复
1stopt/MATLAB 解非线性方程组问题
已经有4人回复
求大神相助,请问这个方程组存不存在解呀?
已经有16人回复
求助:用matlab或其它软件进行如下计算模拟,多谢
已经有3人回复
用MATLAB解方程作图,请大神帮忙
已经有7人回复
求助用matlab三次样条法算一个反应的反应速率
已经有13人回复
1stopt解方程组
已经有12人回复
用1stopt软件做微分方程组的参数估计
已经有13人回复
哪位有1stopt软件把我拟合一下数据好吗?我用matlab非线性曲线拟合lsqnonli误差太大
已经有18人回复
1stopt拟合问题求助
已经有7人回复
这个编号后为什么1stopt没有反应,不能执行
已经有5人回复
一个1stOpt的初级问题:复数方程组计算问题(有代码)
已经有5人回复
(急):用1stopt跑一段程序,解决非线性方程组的解(已有程序,需要较高版本的1stopt
已经有4人回复
求助,用1stopt解方程组,两次得到的解不唯一?什么原因,大神帮忙
已经有6人回复
1stopt中使用sharedmodel拟合曲线出现问题,求指示!!
已经有6人回复
关于matlab微分、及非线性拟合问题
已经有5人回复
1stopt4阶非线性常微分方程,帮小弟看一下,跪求。
已经有7人回复
请教matlab反应动力学参数估计遇到的问题,谢谢
已经有15人回复
1stopt问题请教(附代码)
已经有5人回复
1stOpt求解常微分方程边值问题
已经有5人回复
1楼
2015-01-15 22:31:26
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
ykai_1991
金虫
(初入文坛)
ben_ladeng: 屏蔽内容
2015-01-17 13:11:02
本帖内容被屏蔽
高级回复
4楼
2015-01-16 20:58:17
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
查看全部 12 个回答
ykai_1991
金虫
(初入文坛)
应助: 0
(幼儿园)
金币: 967
帖子: 36
在线: 19.8小时
虫号: 2054876
注册: 2012-10-11
专业: 化学反应工程
本人重新写下代码,有些错误的地方哈
赞
一下
回复此楼
2楼
2015-01-16 13:27:13
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
ykai_1991
金虫
(初入文坛)
应助: 0
(幼儿园)
金币: 967
帖子: 36
在线: 19.8小时
虫号: 2054876
注册: 2012-10-11
专业: 化学反应工程
ben_ladeng: 屏蔽内容
2015-01-17 13:10:30
ben_ladeng: 取消屏蔽
2015-01-17 13:11:09
本人重新写下代码,有些错误的地方哈
对于两种下角标的,这里的kc有方法简化下么?不然就要写10个了
P1~P10也可以这样直接用P(1:10)定义变量么?
Parameters ka(1:10), kc1(1:10), kc2(1:10), kc3(1:10), kc4(1:10), kc5(1:10), kc6(1:10), kc7(1:10), kc8(1:10), kc9(1:10), kc10(1:10) ;
Variable t, P(1:10), Q(1:10), F;
Plot t[x], P1, P2, Q1, Q2, F;
Plot t[x], P3, P4, Q3, Q4;
Plot t[x], P5, P6, Q5, Q6;
Plot t[x], P7, P8, Q7, Q8;
ODEFunction
P1'=-ka1*P1*F-P1*Sum(j=1:10)(kc1[j]*Q[j]);
P2'=-ka2*P2*F-P2*Sum(j=1:10)(kc2[j]*Q[j])+Sum(j=2-1))(kc2[j]*P2*Q[2-j]);
P3'=-ka3*P3*F-P3*Sum(j=1:10)(kc3[j]*Q[j])+Sum(j=3-1))(kc3[j]*P3*Q[3-j]);
P4'=-ka4*P4*F-P4*Sum(j=1:10)(kc4[j]*Q[j])+Sum(j=4-1))(kc4[j]*P4*Q[4-j]);
P5'=-ka5*P5*F-P5*Sum(j=1:10)(kc5[j]*Q[j])+Sum(j=5-1))(kc5[j]*P5*Q[5-j]);
P6'=-ka6*P6*F-P6*Sum(j=1:10)(kc6[j]*Q[j])+Sum(j=6-1))(kc6[j]*P6*Q[6-j]);
P7'=-ka7*P7*F-P7*Sum(j=1:10)(kc7[j]*Q[j])+Sum(j=7-1))(kc7[j]*P7*Q[7-j]);
P8'=-ka8*P8*F-P8*Sum(j=1:10)(kc8[j]*Q[j])+Sum(j=8-1))(kc8[j]*P8*Q[8-j]);
P9'=-ka9*P9*F-P9*Sum(j=1:10)(kc9[j]*Q[j])+Sum(j=9-1))(kc9[j]*P9*Q[9-j]);
P10'=-ka10*P10*F-P10*Sum(j=1:10)(kc10[j]*Q[j])+Sum(j=10-1))(kc10[j]*P10*Q[10-j]);
Q1'=-Q1*Sum(j=1:10)(kc1[j]*P[j])+ka1*P1*F;
Q2'=-Q2*Sum(j=1:10)(kc2[j]*P[j])+ka2*P2*F;
Q3'=-Q3*Sum(j=1:10)(kc3[j]*P[j])+ka3*P3*F;
Q4'=-Q4*Sum(j=1:10)(kc4[j]*P[j])+ka4*P4*F;
Q5'=-Q5*Sum(j=1:10)(kc5[j]*P[j])+ka5*P5*F;
Q6'=-Q6*Sum(j=1:10)(kc6[j]*P[j])+ka6*P6*F;
Q7'=-Q7*Sum(j=1:10)(kc7[j]*P[j])+ka7*P7*F;
Q8'=-Q8*Sum(j=1:10)(kc8[j]*P[j])+ka8*P8*F;
Q9'=-Q9*Sum(j=1:10)(kc9[j]*P[j])+ka9*P9*F;
Q10'=-Q10*Sum(j=1:10)(kc10[j]*P[j])+ka10*P10*F;
F'=-F*Sum(j=1:10)(ka[j]*P[j]);
Data t, P1, Q1, P2;
赞
一下
回复此楼
3楼
2015-01-16 13:32:37
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
ykai_1991
金虫
(初入文坛)
应助: 0
(幼儿园)
金币: 967
帖子: 36
在线: 19.8小时
虫号: 2054876
注册: 2012-10-11
专业: 化学反应工程
除了时间t是min的单位,其他的物质都是mol/L的单位。
赞
一下
回复此楼
5楼
2015-01-16 20:59:23
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
查看全部 12 个回答
不应助
确定回帖应助
(注意:应助才可能被奖励,但不允许灌水,必须填写15个字符以上)
如果回帖内容含有宣传信息,请如实选中。否则帐号将被全论坛禁言
普通表情
龙
兔
虎
猫
SCI英文润色翻译-EnPapers美团队-先修改后付款
百度网盘
|
360云盘
|
千易网盘
|
华为网盘
在新窗口页面中打开自己喜欢的网盘网站,将文件上传后,然后将下载链接复制到帖子内容中就可以了。
最具人气热帖推荐
[查看全部]
作者
回/看
最后发表
[
基金申请
]
希望今年自己国自然面上项目和老婆青年项目能中!
+7
恐龙爸爸
2024-06-14
7/350
2024-06-16 14:48
by
redfish105
[
找工作
]
初始合伙人来啦!(生物试剂耗材标准品)
+7
欢快的小科研人
2024-06-15
12/600
2024-06-16 14:11
by
水水水中游
[
基金申请
]
博后基金,博管会会提前知道消息吗?
+4
yuyiang
2024-06-13
4/200
2024-06-16 11:40
by
yangyuzhong4
[
博后之家
]
山东大学(青岛)“天然药物生物智造”课题组 招聘“博士后”(年薪20.4-55.6万元)
+5
第二种态度
2024-06-11
8/400
2024-06-16 10:14
by
午睡未进行
[
找工作
]
成都产品质量检测研究院
200
+3
鲸鱼663
2024-06-11
9/450
2024-06-16 10:08
by
SNaiL1995
[
基金申请
]
博士后创新人才支持计划公示
+9
aishida144
2024-06-14
15/750
2024-06-16 09:52
by
msjy
[
硕博家园
]
硕博巨婴,也许才刚刚开始
+30
SNaiL1995
2024-06-12
72/3600
2024-06-16 08:34
by
小龙虾biubiubiu
[
基金申请
]
面青地会评时间
+5
tanjydd
2024-06-15
5/250
2024-06-15 23:08
by
andywei1028
[
考博
]
上海交大招收材料化学方向科研助理/“申请考核”博士(请勿回复帖子或站内投条)
+3
灵梦and紫
2024-06-12
4/200
2024-06-15 20:58
by
1822836277
[
教师之家
]
我们学院常年位居 各学院 倒数第二。专业撤销的话,在编者有什么补偿?
+13
河西夜郎
2024-06-09
14/700
2024-06-15 19:44
by
LittleBush
[
基金申请
]
有没有机械的前辈分享一下评上海优都是什么成果啊
+7
wulala800
2024-06-10
7/350
2024-06-15 09:33
by
晓目崇
[
基金申请
]
博后基金,以往的结果点不开,怎么回事呢?最后一次机会了,两次都没中前面。
+7
kyukitu
2024-06-14
13/650
2024-06-15 06:46
by
我是王小帅
[
食品
]
食品博士导师
+6
小李醒yy
2024-06-11
9/450
2024-06-14 23:37
by
小李醒yy
[
考研
]
物理化学一对一辅导
+3
林大diao
2024-06-12
5/250
2024-06-14 20:57
by
林大diao
[
基金申请
]
工材E口JQ有消息了吗
+4
babyduck
2024-06-11
4/200
2024-06-14 17:23
by
firepick
[
基金申请
]
国自然基金公布的时候基金号有吗
+8
潇洒怡惜
2024-06-13
11/550
2024-06-14 11:24
by
JRfei
[
基金申请
]
化学会年会改成一年开一次了?
+5
babu2015
2024-06-12
13/650
2024-06-14 10:46
by
babu2015
[
基金申请
]
75批博后基金
+10
kyukitu
2024-06-13
13/650
2024-06-14 10:31
by
kyukitu
[
考博
]
博导选择
+3
bing85977
2024-06-12
3/150
2024-06-13 15:34
by
我是邱尧
[
论文投稿
]
摩擦磨损论文投稿
+3
jmysan
2024-06-12
3/150
2024-06-13 08:36
by
莱茵润色
信息提示
关闭
请填处理意见
关闭
确定