| 查看: 1678 | 回复: 9 | ||
[求助]
Matlab拟合动力学参数问题，求指导！ 已有1人参与
|
|
小弟现在要用matlab拟合出一个复杂反应体系的动力学参数，先用ode45求解微分方程组再用最小二乘拟合，参照文献编了一个程序，但运行老出一些问题，快急死了… 求大神们帮忙看看，有什么问题…还有个问题，就是我看文献上需要先给动力学参数一个初值，那这个初值该怎么选择啊？ 下面是程序、动力学模型以及实验数据，先谢谢各位了！ clear all; %%清理内存 clc; %%清屏 global C_A_EXP C_B_EXP C_C1_EXP C_C2_EXP C_C3_EXP C_D_EXP C_E_EXP %%各个物质浓度实验值 global C_A_CAL C_B_CAL C_C1_CAL C_C2_CAL C_C3_CAL C_D_CAL C_E_CAL %%各个物质浓度计算值 global R tspan C00 %%理想气体常数、积分区间、初值 R=8.314; tspan={[0 5];[0 10];[0 20];[0 30];[0 45];[0 60];[0 90];[0 120];[0 180]}; %% 读取数据 display '读取50 下数据...';data_expx_50=xlsread('data_expx.xlsx','50','B5:H13'); display '读取60 下数据...';data_expx_60=xlsread('data_expx.xlsx','60','B5: H13'); display '读取70 下数据...';data_expx_70=xlsread('data_expx.xlsx','70','B5: H13'); display '读取80 下数据...';data_expx_80=xlsread('data_expx.xlsx','80','B5: H13'); display '读取初值...' display '读取50 下初值...';data_expx_50_0=xlsread('data_expx.xlsx','50','B4:H4'); display '读取60 下初值...';data_expx_60_0=xlsread('data_expx.xlsx','60','B4:H4'); display '读取70 下初值...';data_expx_70_0=xlsread('data_expx.xlsx','70','B4:H4'); display '读取80 下初值...';data_expx_80_0=xlsread('data_expx.xlsx','80','B4:H4'); display '读取数据结束...' %% 浓度实验值数据 C_A_EXP=[ data_expx_50(:,1)...data_expx_60(:,1)...data_expx_70(:,1)... data_expx_80(:,1)]; C_B_EXP=[ data_expx_50(:,2)...data_expx_60(:,2)...data_expx_70(:,2)... data_expx_80(:,2)]; C_C1_EXP=[ data_expx_50(:,3)...data_expx_60(:,3)...data_expx_70(:,3)... data_expx_80(:,3)]; C_C2_EXP=[ data_expx_50(:,4)...data_expx_60(:,4)...data_expx_70(:,4)... data_expx_80(:,4)]; C_C3_EXP=[ data_expx_50(:,5)...data_expx_60(:,5)...data_expx_70(:,5)... data_expx_80(:,5)]; C_D_EXP=[ data_expx_50(:,6)...data_expx_60(:,6)...data_expx_70(:,6)... data_expx_80(:,6)]; C_E_EXP=[ data_expx_50(:,7)...data_expx_60(:,7)...data_expx_70(:,7)... data_expx_80(:,7)]; %% 初值 C_A_EXP_0=[data_expx_50_0(:,1)...data_expx_60_0(:,1)...data_expx_70_0(:,1)... data_expx_80_0(:,1)]; C_B_EXP_0=[data_expx_50_0(:,2)...data_expx_60_0(:,2)...data_expx_70_0(:,2)... data_expx_80_0(:,2)]; C_C1_EXP_0=[data_expx_50_0(:,3)...data_expx_60_0(:,3)...data_expx_70_0(:,3)... data_expx_80_0(:,3)]; C_C2_EXP_0=[data_expx_50_0(:,4)...data_expx_60_0(:,4)...data_expx_70_0(:,4)... data_expx_80_0(:,4)]; C_C3_EXP_0=[data_expx_50_0(:,5)...data_expx_60_0(:,5)...data_expx_70_0(:,5)... data_expx_80_0(:,5)]; C_D_EXP_0=[data_expx_50_0(:,6)...data_expx_60_0(:,6)...data_expx_70_0(:,6)... data_expx_80_0(:,6)]; C_E_EXP_0=[data_expx_50_0(:,7)...data_expx_60_0(:,7)...data_expx_70_0(:,7)... data_expx_80_0(:,7)]; C00=[C_A_EXP_0 C_B_EXP_0 C_C1_EXP_0 C_C2_EXP_0 C_C3_EXP_0 C_D_EXP_0 C_E_EXP_0]; p0=[2.0e+5 1.0e+004 1.0... 2.0e+5 1.0e+004 1.0... 2.0e+5 1.0e+004 1.0... 2.0e+5 1.0e+004 1.0... 2.0e+5 1.0e+004 1.0... 2.0e+5 1.0e+004 1.0... 2.0e+5 1.0e+004 1.0... 2.0e+5 1.0e+004 1.0]; options=optimset('tolfun',0.1,'TolX',1e4,'MaxFunEvals',9e25,'display','iter'); [p,exitflag]=fminsearch(@minfun_,p0,options); %% 显示结果 disp 'p1=';p(1) disp 'p2=';p(2) disp 'p3=';p(3) disp 'p4=';p(4) disp 'p5=';p(5) disp 'p6=';p(6) disp 'p7=';p(7) disp 'p8=';p(8) disp 'p9=';p(9) disp 'p10=';p(10) disp 'p11=';p(11) disp 'p12=';p(12) disp 'p13=';p(13) disp 'p14=';p(14) disp 'p15=';p(15) disp 'p16=';p(16) disp 'p17=';p(17) disp 'p18=';p(18) disp 'p19=';p(19) disp 'p20=';p(20) disp 'p21=';p(21) disp 'p22=';p(22) disp 'p23=';p(23) disp 'p24=';p(24) function dCdt=myfun_(t,C,p,m) global R %% %% 计算各个物质反应速率 %% if m<=1 TC=323.15; elseif m>1 && m<=2 TC=433.15; elseif m>2 && m<=3 TC=443.15; elseif m>3 TC=453.15; end R1=p(1).*exp(-p(2)./(R*TC)).*(C(1).^p(3)); R2=p(4).*exp(-p(5)./(R*TC)).*(C(2).^p(6)); R3=p(7).*exp(-p(8)./(R*TC)).*(C(2).^p(9)); R4=p(10).*exp(-p(11)./(R*TC)).*(C(2).^p(12)); R5=p(13).*exp(-p(14)./(R*TC)).*(C(3).^p(15)); R6=p(16).*exp(-p(17)./(R*TC)).*(C(4).^p(18)); R7=p(19).*exp(-p(20)./(R*TC)).*(C(5).^p(21)); R8=p(22).*exp(-p(23)./(R*TC)).*(C(6).^p(24)); dCdt=[-R1;R1-R2-R3-R4;R2-R5;R3-R6;R4-R7;R5+R6+R7-R8;R8]; function min=minfun_(p) global R tspan C00 for i=1:1:9 for m=1:1:7 C0=[C00(m) C00(m+7) C00(m+14) C00(m+21)]; [t,C]=ode45(@myfun,tspan{i},C0,[],p,m); C_A_CAL(i,m)=C(end,1); C_B_CAL(i,m)=C(end,2); C_C1_CAL(i,m)=C(end,3); C_C2_CAL(i,m)=C(end,4); C_C3_CAL(i,m)=C(end,5); C_D_CAL(i,m)=C(end,6); C_E_CAL(i,m)=C(end,7); end end min = sum(sum((C_A_EXP-C_A_CAL).^2+(C_B_EXP-C_B_CAL).^2+… (C_C1_EXP-C_C1_CAL).^2+(C_C2_ EXP-C_C2_CAL).^2+… (C_C3_EXP-C_C3_CAL).^2+(C_D_EXP-C_D_CAL).^2+… (C_E_ EXP-C_E_CAL).^2));  动力学模型1.jpg  动力学模型2.jpg |
回复此楼» 本帖附件资源列表
-
欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com - 附件 1 : data_expx.xlsx
2015-11-06 15:59:20, 14.03 K
» 猜你喜欢
070200求调剂，一志愿某211，288分
已经有20人回复
-
法国博士后职位
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有67人回复
-
重庆交大26年硕士生招生拟调剂通知已出！欢迎加入机器视觉与3D光学成像课题组。
已经有0人回复
-
**
已经有1人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有3人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有2人回复
-
B区学生调剂-兰州交通大学材料科学与工程学院
已经有8人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有9人回复
-
桂林理工大学物理学专业招收调剂
已经有6人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 向大神求教一个用matlab拟合动力学参数的问题，求指导
已经有5人回复
- 给我大神求教一个用matlab拟合动力学参数的问题，求指导
已经有0人回复
- matlab 拟合多参数动力学方程 求程序
已经有4人回复
- MATLAB拟合动力学问题
已经有6人回复
- MATLAB拟合动力学问题
已经有11人回复
- matlab拟合动力学参数，非常谢谢~
已经有0人回复
- 求助matlab拟合动力学K值
已经有4人回复
- matlab拟合吸附动力学
已经有0人回复
- 使用MATLAB拟合动力学参数问题
已经有18人回复
- 拜求！matlab拟合动力学参数
已经有14人回复
- 求助用matlab拟合动力学方程
已经有9人回复
- 求助用matlab拟合动力学方程
已经有13人回复
- matlab 拟合反应动力学参数
已经有8人回复
- matlab拟合酶促反应机制
已经有7人回复
- 求大家帮帮忙看看matlab拟合动力学方程的问题，已经走投无路了，求大家帮帮忙啊！
已经有19人回复
- matlab 拟合反应动力学参数结果很差。大家帮忙看一下
已经有14人回复
- 用matlab 拟合动力学级数
已经有2人回复
- matlab拟合拟均相动力学参数
已经有6人回复
- matlab拟合动力学参数
已经有4人回复
- 【求助】使用Matlab拟合反应动力学方程问题
已经有7人回复
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
2楼2015-11-06 17:04:12
3楼2015-11-06 19:26:17
hujinqiang
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 5.8
- 散金: 5
- 帖子: 15
- 在线: 29.7小时
- 虫号: 3961731
- 注册: 2015-07-09
- 专业: 控制理论与方法
4楼2015-12-12 16:52:46
5楼2015-12-14 10:27:29
6楼2015-12-14 10:49:23
7楼2015-12-14 10:55:36
8楼2015-12-14 11:03:39
9楼2015-12-14 11:05:51
10楼2015-12-16 00:11:16
