应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » Matlab拟合动力学参数问题，求指导！
51/1返回列表
查看: 1650  |  回复: 9
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

ZHM10

新虫 (小有名气)

[求助] 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));Matlab拟合动力学参数问题，求指导！
动力学模型1.jpg


Matlab拟合动力学参数问题，求指导！-1
动力学模型2.jpg
回复此楼

» 本帖附件资源列表

  • 欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
    本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com
  • 附件 1 : data_expx.xlsx
    • 2015-11-06 15:59:20, 14.03 K

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2015-11-06 16:05:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

898766282

新虫 (正式写手)
引用回帖:
4楼: Originally posted by hujinqiang at 2015-12-12 16:52:46
你好，我现在有一个复杂的动力学方程，不知道怎么解，你可否指导一下

我最近也在编写求复杂反应动力学的matlabb程序，但是我准备用的是退货算法，不过还没写出来，可以发给我你的参考资料文献吗？谢谢!
赞一下 回复此楼
5楼2015-12-14 10:27:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 10 个回答

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
ZHM10: 金币+10 2015-11-08 16:11:38
这种微分方程拟合问题建议用1stOpt来解决，代码简单十倍不止。
赞一下 回复此楼
2楼2015-11-06 17:04:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ZHM10

新虫 (小有名气)
引用回帖:
2楼: Originally posted by dingd at 2015-11-06 17:04:12
这种微分方程拟合问题建议用1stOpt来解决，代码简单十倍不止。

谢谢！你能帮我看一下matlab程序吗，直接运行就行，我看文章上都是用这种方法，，1stopt这东西一点都没接触过啊
赞一下 回复此楼
3楼2015-11-06 19:26:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hujinqiang

新虫 (初入文坛)
你好，我现在有一个复杂的动力学方程，不知道怎么解，你可否指导一下
赞一下 回复此楼
4楼2015-12-12 16:52:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 10 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 一志愿 西北大学 ，070300化学学硕，总分287，双非一本，求调剂。 +4 晨昏线与星海 2026-03-19 4/200 2026-03-20 22:15 by JourneyLucky
[考研] 一志愿中南化学（0703）总分337求调剂 +8 niko- 2026-03-19 9/450 2026-03-20 21:57 by luoyongfeng
[考研] A区线材料学调剂 +5 周周无极 2026-03-20 5/250 2026-03-20 21:33 by laoshidan
[考研] 材料与化工专硕调剂 +7 heming3743 2026-03-16 7/350 2026-03-20 19:31 by zhukairuo
[考研] 08工学调剂 +5 用户573181 2026-03-20 5/250 2026-03-20 15:47 by xia_2003
[考研] 265求调剂 +8 梁梁校校 2026-03-17 8/400 2026-03-20 14:40 by 27道科特
[考研] 0856调剂，是学校就去 +8 sllhht 2026-03-19 9/450 2026-03-20 14:25 by 无懈可击111
[考博] 申博26年 +3 八6八68 2026-03-19 3/150 2026-03-19 19:43 by nxgogo
[考研] 085600材料与化工求调剂 +6 绪幸与子 2026-03-17 6/300 2026-03-19 13:27 by houyaoxu
[考研] 344求调剂 +6 knight344 2026-03-16 7/350 2026-03-18 20:13 by walc
[考研] 070300化学319求调剂 +6 锦鲤0909 2026-03-17 6/300 2026-03-18 13:22 by Iveryant
[考研] 0703化学调剂 +3 妮妮ninicgb 2026-03-17 3/150 2026-03-18 10:29 by macy2011
[考研] 301求调剂 +9 yy要上岸呀 2026-03-17 9/450 2026-03-18 08:58 by 无际的草原
[考研] 334求调剂 +3 志存高远意在机� 2026-03-16 3/150 2026-03-18 08:34 by lm4875102
[考研] 有没有道铁/土木的想调剂南林，给自己招师弟中～ +3 TqlXswl 2026-03-16 7/350 2026-03-17 15:23 by TqlXswl
[考研] 一志愿苏州大学材料工程（085601）专硕有科研经历三项国奖两个实用型专利一项省级立项 +6 大火山小火山 2026-03-16 8/400 2026-03-17 15:05 by 无懈可击111
[考研] 考研调剂 +3 淇ya_~ 2026-03-17 5/250 2026-03-17 09:25 by Winj1e
[考研] 11408 一志愿西电，277分求调剂 +3 zhouzhen654 2026-03-16 3/150 2026-03-17 07:03 by laoshidan
[考研] 求老师收留调剂 +4 jiang姜66 2026-03-14 5/250 2026-03-15 20:11 by Winj1e
[考研] 一志愿哈工大材料324分求调剂 +5 闫旭东 2026-03-14 5/250 2026-03-14 14:53 by 木瓜膏
信息提示
关闭
请填处理意见
关闭