| 查看: 206 | 回复: 0 | |||
| 当前主题已经存档。 | |||
| 【有奖交流】积极回复本帖子，参与交流，就有机会分得作者 只做陌生人 的 40 个金币 | |||
[交流]
【求助】新手求帮助解答 动力学方程的参数估计方面的问题
|
|||
|
小弟最近在做动力学方面的模拟，即用经验型的方程式，里面有8个参数需要估计，前阵子抢了有好多数据了，然后最近在急着写程序，由于是新手，做得比较慢，希望高手帮忙，话正题： 该动力学是两个反应，A+B=C,C+A=D，现在要得到A和B的反应速率，我在程序中是用f1和f2表示的，嗯，我附上我的程序吧，希望高手帮忙指导调试一下，不胜感谢啊。 function mykinetics clc global T,obj,g,w,x g0=[300,30000,0.9,0.8,400,15000,0.6,0.86]; [g,obj]=fminsearch(wlq_obj,g0); %Cc=0; T=[438.1500 443.1500 448.1500 453.1500 458.1500 463.1500]; C_e1=[ 8.2209 0.4252 8.2327 0.4370 7.9779 0.1650 8.1535 0.3579 8.1660 0.3760 8.1419 0.3528 ]; C_e2=[ 8.4981 0.2077 8.5034 0.2095 8.5025 0.2060 8.4996 0.2148 8.5019 0.2275 8.4710 0.1872 ]; C_e3=[ 8.8011 0.1389 8.7987 0.1389 8.7816 0.1233 8.7770 0.1158 8.7783 0.1165 8.7570 0.0959 ]; C_e4=[ 9.0110 0.0577 8.9832 0.0319 8.9843 0.0336 8.9669 0.0189 8.9732 0.0258 8.9644 0.0186 ]; plot(T,Cc1,'*',T,C_e1,'+') ru=1-(nansum((Cc1-C_e1).^2)+nansum((Cc2-C_e2).^2)+nansum((Cc3-C_e3).^2)+nansum((Cc4-C_e4).^2)+nansum((Cc5-C_e5).^2))/(nansum(C_e1.^2)+nansum(C_e2.^2)+nansum(C_e3.^2)); error1=nansum((C_e1-Cc1).^2); s1=nansum(C_e1.^2)-error1; error2=nansum((C_e2-Cc2).^2); s2=nansum(C_e2.^2)-error2; error3=nansum((C_e3-Cc3).^2); s3=nansum(C_e3.^2)-error3; error4=nansum((C_e4-Cc4).^2); s4=nansum(C_e4.^2)-error4; s=s1+s2+s3+s4; error=error1+error2+error3+error4; f=s*(24-8-1)/(error*8); result Cc1 C_e1 ru s error f %-------------------------------------------------------------------------- function obj=wlq_obj(g) global T,obj,g T=[438.1500 443.1500 448.1500 453.1500 458.1500 463.1500]; C_e1=[ 8.2209 0.4252 8.2327 0.4370 7.9779 0.1650 8.1535 0.3579 8.1660 0.3760 8.1419 0.3528 ]; C_e2=[ 8.4981 0.2077 8.5034 0.2095 8.5025 0.2060 8.4996 0.2148 8.5019 0.2275 8.4710 0.1872 ]; C_e3=[ 8.8011 0.1389 8.7987 0.1389 8.7816 0.1233 8.7770 0.1158 8.7783 0.1165 8.7570 0.0959 ]; C_e4=[ 9.0110 0.0577 8.9832 0.0319 8.9843 0.0336 8.9669 0.0189 8.9732 0.0258 8.9644 0.0186 ]; obj1=sum((C_e1-wlq_dynamic_equation1(T,g)).^2); obj2=sum((C_e2-wlq_dynamic_equation2(T,g)).^2); obj3=sum((C_e3-wlq_dynamic_equation3(T,g)).^2); obj4=sum((C_e4-wlq_dynamic_equation4(T,g)).^2); obj=obj1+obj2+obj3+obj4 %-------------------------------------------------------------------------- function C=wlq_dynamic_equation1(T,g) global k01,k02,E1,E2,m1,n1,m1,n2,w,x,Cc1 k01=g(1); k02=g(5); E1=g(2); E2=g(6); m1=g(3); n1=g(4); m2=g(7); n2=g(8); [w,x]=ode45(@wlq_f1,[0,0.222],[9.5074 1.9007]); for i=1:length(T) Cc1(i,  =wlq_dynamic_equation1(T(i),result);end %-------------------------------------------------------------------------- function C=wlq_dynamic_equation2(T,g) global k01,k02,E1,E2,m1,n1,m1,n2,w,x,Cc2 k01=g(1); k02=g(5); E1=g(2); E2=g(6); m1=g(3); n1=g(4); m2=g(7); n2=g(8); [w,x]=ode45(@wlq_f2,[0,0.222],[9.7582 1.6262]); for i=1:length(T) Cc2(i, =wlq_dynamic_equation1(T(i),result);end %-------------------------------------------------------------------------- function C=wlq_dynamic_equation3(T,g) global k01,k02,E1,E2,m1,n1,m1,n2,w,x,Cc3 k01=g(1); k02=g(5); E1=g(2); E2=g(6); m1=g(3); n1=g(4); m2=g(7); n2=g(8); [w,x]=ode45(@wlq_f3,[0,0.222],[9.9463 1.4203]); for i=1:length(T) Cc3(i, =wlq_dynamic_equation1(T(i),result);end %-------------------------------------------------------------------------- function C=wlq_dynamic_equation4(T,g) global k01,k02,E1,E2,m1,n1,m1,n2,w,x,Cc4 k01=g(1); k02=g(5); E1=g(2); E2=g(6); m1=g(3); n1=g(4); m2=g(7); n2=g(8); [w,x]=ode45(@wlq_f4,[0,0.222],[10.0904 1.2625]); for i=1:length(T) Cc4(i, =wlq_dynamic_equation1(T(i),result);end %-------------------------------------------------------------------------- function f=wlq_f1(w,x) global k01,k02,E1,E2,m1,n1,m1,n2,k1,k2,x01,x02,R,f1,f2,x1,x2 %Cb=x1,Cpr=x2 R=8.314; % （单位J?mol-1?K-1） x01=9.5074; % （单位mol?L-1） x02=1.9007;% （单位mol?L-1） k1=k01.*exp(-(E1./(R.*T))); k2=k02.*exp(-(E2./(R.*T))); f1=k1.*x1.^m1.*x2.^n1+k2.*(2*(x01-x1)-(x02-x2)).^m2.*x2.^n2; f2=k1.*x1.^m1.*x2.^n1+2*k2.*(2*(x01-x1)-(x02-x2)).^m2.*x2.^n2; %-------------------------------------------------------------------------- function f=wlq_f2(w,x) global k01,k02,E1,E2,m1,n1,m1,n2,k1,k2,x01,x02,R,f1,f2,x1,x2 %Cb=x1,Cpr=x2 R=8.314; % （单位J?mol-1?K-1） x01=9.7582; % （单位mol?L-1） x02=1.6262;% （单位mol?L-1） k1=k01.*exp(-(E1./(R.*T))); k2=k02.*exp(-(E2./(R.*T))); f1=k1.*x1.^m1.*x2.^n1+k2.*(2*(x01-x1)-(x02-x2)).^m2.*x2.^n2; f2=k1.*x1.^m1.*x2.^n1+2*k2.*(2*(x01-x1)-(x02-x2)).^m2.*x2.^n2; %-------------------------------------------------------------------------- function f=wlq_f3(w,x) global k01,k02,E1,E2,m1,n1,m1,n2,k1,k2,x01,x02,R,f1,f2,x1,x2 %Cb=x1,Cpr=x2 R=8.314; % （单位J?mol-1?K-1） x01=9.9463; % （单位mol?L-1） x02=1.4203;% （单位mol?L-1） k1=k01.*exp(-(E1./(R.*T))); k2=k02.*exp(-(E2./(R.*T))); f1=k1.*x1.^m1.*x2.^n1+k2.*(2*(x01-x1)-(x02-x2)).^m2.*x2.^n2; f2=k1.*x1.^m1.*x2.^n1+2*k2.*(2*(x01-x1)-(x02-x2)).^m2.*x2.^n2; %-------------------------------------------------------------------------- function f=wlq_f4(w,x) global k01,k02,E1,E2,m1,n1,m1,n2,k1,k2,x01,x02,R,f1,f2,x1,x2 %Cb=x1,Cpr=x2 R=8.314; % （单位J/mol-1/K-1） x01=10.0904; % （单位mol/L-1） x02=1.2625;% （单位mol/L-1） k1=k01.*exp(-(E1./(R.*T))); k2=k02.*exp(-(E2./(R.*T))); f1=k1.*x1.^m1.*x2.^n1+k2.*(2*(x01-x1)-(x02-x2)).^m2.*x2.^n2; f2=k1.*x1.^m1.*x2.^n1+2*k2.*(2*(x01-x1)-(x02-x2)).^m2.*x2.^n2; 我运行的时候，总提示我变量没有定义，但是前面明明已经定义过了的啊，不明白～～～ [ Last edited by kuhailangyu on 2010-3-24 at 21:48 ] |
回复此楼» 猜你喜欢
售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急
已经有3人回复
-
国自然面上复盘~欢迎讨论
已经有16人回复
-
售SCI一区文章，我:8.O.551.O.5.4,科目全,可伽急
已经有3人回复
-
售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急
已经有4人回复
-
有没有仍没收到信息的
已经有6人回复
-
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急
已经有4人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有7人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有6人回复
-
售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急
已经有8人回复
-
面上意见出来了
已经有12人回复
30