版块导航
正在加载中...
客户端APP下载
论文辅导
登录
注册
帖子
帖子
用户
本版
应《网络安全法》要求,自2017年10月1日起,未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用,请尽快对帐号进行手机号验证,感谢您的理解与支持!
24小时热门版块排行榜
>
论坛更新日志
(5604)
>
虫友互识
(1061)
>
休闲灌水
(226)
>
文献求助
(190)
>
招聘信息布告栏
(184)
>
论文投稿
(128)
>
基金申请
(119)
>
导师招生
(95)
>
找工作
(95)
>
硕博家园
(67)
>
职场人生
(64)
>
考博
(62)
>
博后之家
(44)
>
绿色求助(高悬赏)
(36)
>
育儿交流
(22)
>
留学DIY
(15)
小木虫论坛-学术科研互动平台
»
计算模拟区
»
计算模拟
»
1stOpt编写常微分方程组计算反应速率问题
5
1/1
返回列表
查看: 1112 | 回复: 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论文润色/翻译怎么收费?
已经有253人回复
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
金虫
(初入文坛)
应助: 0
(幼儿园)
金币: 967
帖子: 36
在线: 19.8小时
虫号: 2054876
注册: 2012-10-11
专业: 化学反应工程
引用回帖:
8楼
:
Originally posted by
CelestialCYJ
at 2015-01-19 17:13:07
你P3~10,Q2~10在0时刻的值是不是都等于0???
是的是的。
回复此楼
高级回复
9楼
2015-01-20 13:06:19
已阅
回复此楼
关注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
金虫
(初入文坛)
ben_ladeng: 屏蔽内容
2015-01-17 13:11:02
本帖内容被屏蔽
4楼
2015-01-16 20:58:17
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
查看全部 12 个回答
不应助
确定回帖应助
(注意:应助才可能被奖励,但不允许灌水,必须填写15个字符以上)
如果回帖内容含有宣传信息,请如实选中。否则帐号将被全论坛禁言
普通表情
龙
兔
虎
猫
论文翻译-EnPapers美助理教授团队-SCI英文润色
百度网盘
|
360云盘
|
千易网盘
|
华为网盘
在新窗口页面中打开自己喜欢的网盘网站,将文件上传后,然后将下载链接复制到帖子内容中就可以了。
最具人气热帖推荐
[查看全部]
作者
回/看
最后发表
[
基金申请
]
F口信息学部拿面上,大概需要什么样的成果
+3
_奋黎_
2024-06-16
3/150
2024-06-16 15:38
by
我是王小帅
[
教师之家
]
每次骚扰女学生的都是院系领导,而不是普通教师,小编们要注意措辞正确
+9
zju2000
2024-06-15
11/550
2024-06-16 14:49
by
appleapple2
[
基金申请
]
希望今年自己国自然面上项目和老婆青年项目能中!
+7
恐龙爸爸
2024-06-14
7/350
2024-06-16 14:48
by
redfish105
[
找工作
]
成都产品质量检测研究院
200
+3
鲸鱼663
2024-06-11
9/450
2024-06-16 10:08
by
SNaiL1995
[
考博
]
34岁读博士晚吗
+26
emitdne
2024-06-13
26/1300
2024-06-16 07:16
by
liyeqik
[
文学芳草园
]
累并快乐着
+13
MYHLD521
2024-06-14
13/650
2024-06-15 22:59
by
zeolitess
[
基金申请
]
为什么我的博后基金还在流动站审核中?不会是学院给我卡了吧?
+14
王凯12
2024-06-13
26/1300
2024-06-15 15:22
by
好人与坏人
[
教师之家
]
饶议:什么制度能保障大学普通教师不用为领导拎包,不用看领导脸色
+8
zju2000
2024-06-12
14/700
2024-06-15 13:59
by
chemhua
[
论文投稿
]
求机械类四区sci推荐
5
+3
迷茫小旷
2024-06-14
4/200
2024-06-15 11:25
by
bobvan
[
食品
]
食品博士导师
+6
小李醒yy
2024-06-11
9/450
2024-06-14 23:37
by
小李醒yy
[
论文投稿
]
审稿问题:为什么荧光激发波长和紫外吸收波长差的大?
10
+4
sdawege
2024-06-14
8/400
2024-06-14 22:39
by
东北读书人
[
基金申请
]
博士后基金需要结题吗?
+8
zhouchuck
2024-06-13
8/400
2024-06-14 17:27
by
liuyupu132
[
基金申请
]
E12面上申请
+4
汉风之遗
2024-06-13
4/200
2024-06-14 15:28
by
天外飞去来
[
硕博家园
]
关于硕博连读的一些疑问?
+4
Lwenter
2024-06-14
4/200
2024-06-14 14:32
by
ou0551
[
论文投稿
]
ACS Nano投稿后分配到副编辑手里12天了,能确定送审了吗?
+5
潇洒怡惜
2024-06-12
10/500
2024-06-14 09:56
by
潇洒怡惜
[
论文投稿
]
文章proof要求使用机构的邮箱
5
+3
不可不信缘
2024-06-11
11/550
2024-06-14 07:00
by
3001160025
[
基金申请
]
连续两年医学口青年项目初审体会
+11
进击的荣耀
2024-06-09
18/900
2024-06-13 17:27
by
进击的荣耀
[
考博
]
博导选择
+3
bing85977
2024-06-12
3/150
2024-06-13 15:34
by
我是邱尧
[
论文投稿
]
with editor日期变更
+3
慎独的小花卷
2024-06-12
8/400
2024-06-13 11:00
by
慎独的小花卷
[
基金申请
]
博后特助这周出结果吗?往年都是啥时候啊?
+13
jsqy
2024-06-12
17/850
2024-06-12 19:55
by
Lynn212
信息提示
关闭
请填处理意见
关闭
确定