[求助]matlab运行了10个小时还不出结果

mathtype · 发表于 2005-12-24 13:52

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

我把几个小程序封装在一起后，运行10个小时却还不出结果，但把ode45改为ode23s结果2分钟就出来了，程序如下： function Simulate
clear all;clc format long
R=8.31434; t=700; % 气体温度，t: Temperatre, Celsius T = t + 273.15; % T: Temperatre, K M=101325/(R*T) % 第三体浓度M xCH4=0.2;xCO2=1-xCH4; % xCH4为甲烷的摩尔百分含量 cCH4=xCH4*M,cCO2=xCO2*M, Y0=[cCH4 cCO2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; % 各物质的初始浓度 [t,Y]=ode23(@f,[0:0.005:0.05],Y0,[],T,M); disp('Results by using ode45()：') disp([t Y]) plot(t,Y),xlabel('t'),ylabel('Y') legend('CH4','CO2','H2','CO','H2O','O2','CH3','CH2','CH','C','H','O','OH','HCO','CH2O','CH3O','HO2') Cbanlance=Y0(1)+Y0(2)-Y(:,1)-Y(:,2)-Y(:,4)-Y(:,7)-Y(:,8)-Y(:,9)-Y(:,10)-Y(:,14)-Y(:,15)-Y(:,16); Hbanlance=4*Y0(1)-4*Y(:,1)-2*Y(:,3)-2*Y(:,5)-3*Y(:,7)-2*Y(:,8)-Y(:,9)-Y(:,11)-Y(:,13)-Y(:,14)-2*Y(:,15)-3*Y(:,16)-Y(:,17); Obanlance=2*Y0(2)-2*Y(:,2)-Y(:,4)-Y(:,5)-2*Y(:,6)-Y(:,12)-Y(:,13)-Y(:,14)-Y(:,15)-Y(:,16)-2*Y(:,17); disp('Cbanlance,Hbanlance,Obanlance:') disp([Cbanlance Hbanlance Obanlance]) xlswrite('ode45forSimulate3',[t Y]); %-------------------------------------------------------------------------- function dYdt=f(t,Y,T,M)
NA=6.02217e+23; electronm=9.10956e-31;kb=1.38062e-23;electronC=1.60219e-19; electrondensity=7e+11;Ne=electrondensity*1e+06/NA; % 电子密度electrondensity Te=2; % 平均电子温度，单位eV
% ------------------------------------------------------------------------- Tch4 = Te*11600; % T: Temperatre, K % 求甲烷与电子碰撞的反应 K01=Krate(5.16e-05,-0.642,169882.6,Tch4);K02=Krate(3.02e-04,-0.781,190073.5,Tch4); K03=Krate(5.75e-04,-0.828,204462.5,Tch4);K04=Krate(6.29e-16,0,0,Tch4); % -------------------------------------------------------------------------
% ------------------------------------------------------------------------- % 求二氧化碳与电子碰撞的反应 Eth = 1.1e-02; % 积分下限,单位keV Eupper = 5.0e-02; % 积分下限，单位keV f3_quadl = quadl(@func,Eth,Eupper,1e-30,[],Eth,Te); K05=((8*electronC/(pi*electronm))^(1/2))*((Te)^(-3/2))*f3_quadl*(1e-4)*NA; % -------------------------------------------------------------------------
% ------------------------------------------------------------------------- % 求自由基反应速率常数 K1 = rate1(7.23e+14,0.0,14970,T,2);K2 = rate1(4.10e+14,0.0,13970,T,2);K3 = rate1(2.83e+13,0.0,4970,T,2); K4 = rate1(2.60e+28,-5.11,-2627,T,2);K5 = rate1(8.00e+12,0.0,0,T,2);K6 = rate1(8.00e+13,0.0,0,T,2); K7 = rate1(3.40e+11,0.0,8940,T,2);K8 = rate1(2.00e+13,0.0,0,T,2);K9 = rate1(1.00e+12,0.0,400,T,2); K10 = rate1(4.00e+13,0.0,0,T,2);K11 = rate1(1.00e+14,0.0,3700,T,2);K12 = rate1(8.64e+10,0.0,-500,T,2); K13 = rate1(1.59e+12,0.0,1000,T,2);K14 = rate2(2.50e+13,0.00,0.0,T,2);K15 = rate1(2.00e+13,0.0,0,T,2); K16 = rate1(1.00e+13,0.0,0,T,2);K17 = rate1(4.00e+13,0.0,0,T,2);K18 = rate1(5.00e+13,0.0,2100,T,2); K19 = rate1(1.00e+12,0.0,6000,T,2);K20 = rate1(3.50e+13,0.0,3510,T,2);K21 = rate1(2.50e+13,0.0,3990,T,2); K22 = rate1(7.53e+12,0.0,170,T,2);K23 = rate1(1.00e+10,0.5,6000,T,2);K24 = rate1(3.55e+14,0.0,16790,T,2); K25 = rate1(5.00e+13,0.0,0,T,2);K26 = rate1(2.00e+14,0.0,0,T,2);K27 = rate1(3.00e+13,0.0,0,T,2); K28 = rate1(3.00e+13,0.0,0,T,2);K29 = rate1(3.00e+12,0.0,0,T,2);K30 = rate1(5.93e+8,1.05,-613,T,2); K31 = rate1(3.00e+11,0.5,0,T,2);K32 = rate1(2.20e+22,-2.0,0,T,3);K33 = rate1(1.00e+16,0.0,0,T,3); K34 = rate1(4.70e+15,-0.28,0,T,3);K35 = rate1(1.83e+18,-1.0,0,T,3);K36 = rate1(2.30e+18,-0.8,0,T,3); K37 = rate1(1.50e+9,1.14,0,T,2);K38 = rate1(1.00e+13,0.0,0,T,2);K39 = rate1(8.30e+9,1.0,6950,T,2); K40 = rate1(1.80e+10,1.0,8910,T,2);K41 = rate1(4.40e+6,1.5,-740,T,2);K42 = rate1(1.20e+9,1.3,3630,T,2); K43 = rate1(1.50e+14,0.84,1000,T,2);K44 = rate1(2.50e+13,0.0,700,T,2);K45 = rate1(2.00e+13,0.0,0,T,2); K46 = rate1(1.50e+13,0.0,0,T,2);K47 = rate3(6.3e+11,0.5,0,T,2);K48 = rate3(6.3e+11,0.5,0,T,2); K49 = rate2(4.80e+01,3.22,-13.5,T,2); %-------------------------------------------------------------------------- % 解常微分方程组 f1=K4*Y(7)*Y(11)+K9*Y(7)*Y(17)+K23*Y(15)*Y(7)+K31*Y(14)*Y(7)-(K01+K02+K03+K04)*Y(1)*Ne-K1*Y(1)*Y(11)-K2*Y(1)*Y(12)-K3*Y(1)*Y(13); f2=K13*Y(8)*Y(6)+K28*Y(14)*Y(12)+K41*Y(4)*Y(13)-K05*Y(2)*Ne-K49*Y(2)*Y(9); f3=(K02+K03)*Y(1)*Ne+K1*Y(1)*Y(11)+K5*Y(7)*Y(13)+K10*Y(8)*Y(11)+K21*Y(15)*Y(11)+K26*Y(14)*Y(11)+K35*Y(11)*Y(11)*M+K39*Y(11)*Y(13)+K44*Y(17)*Y(11)-K40*Y(3)*Y(12)-K42*Y(13)*Y(3); f4=K05*Y(2)*Ne+K12*Y(8)*Y(6)+K15*Y(9)*Y(6)+K17*Y(9)*Y(12)+K24*Y(14)*M+K25*Y(14)*Y(13)+K26*Y(14)*Y(11)+K27*Y(14)*Y(12)+K29*Y(14)*Y(6)+K31*Y(14)*Y(7)+K47*Y(10)*Y(6)+K48*Y(10)*Y(13)+K49*Y(2)*Y(9)-K41*Y(4)*Y(13); f5=K3*Y(1)*Y(13)+K22*Y(15)*Y(13)+K25*Y(14)*Y(13)+K32*Y(13)*Y(11)*M+K37*Y(13)*Y(13)+K42*Y(13)*Y(3)+K46*Y(17)*Y(13); f6=K9*Y(7)*Y(17)+K30*Y(14)*Y(17)+K34*Y(12)*Y(12)*M+K38*Y(12)*Y(13)+K44*Y(17)*Y(11)+K45*Y(17)*Y(12)+K46*Y(17)*Y(13)-K7*Y(7)*Y(6)-K11*Y(8)*Y(6)-K12*Y(8)*Y(6)-K13*Y(8)*Y(6)-K15*Y(9)*Y(6)-K16*Y(9)*Y(6)-K19*Y(16)*Y(6)-K29*Y(14)*Y(6)-K36*Y(11)*Y(6)*M-K47*Y(10)*Y(6); f7=K01*Y(1)*Ne+K1*Y(1)*Y(11)+K2*Y(1)*Y(12)+K3*Y(1)*Y(13)-K4*Y(7)*Y(11)-K5*Y(7)*Y(13)-K6*Y(7)*Y(12)-K7*Y(7)*Y(6)-K8*Y(7)*Y(17)-K9*Y(7)*Y(17)-K23*Y(15)*Y(7)-K31*Y(14)*Y(7); f8=K02*Y(1)*Ne-K10*Y(8)*Y(11)-K11*Y(8)*Y(6)-K12*Y(8)*Y(6)-K13*Y(8)*Y(6)-K14*Y(8)*Y(13); f9=K03*Y(1)*Ne+K10*Y(8)*Y(11)-K15*Y(9)*Y(6)-K16*Y(9)*Y(6)-K17*Y(9)*Y(12)-K49*Y(2)*Y(9); f10=K04*Y(1)*Ne-K47*Y(10)*Y(6)-K48*Y(10)*Y(13); f11=(K01+K03+4*K04)*Y(1)*Ne+K6*Y(7)*Y(12)+K12*Y(8)*Y(6)+2*K13*Y(8)*Y(6)+K14*Y(8)*Y(13)+K17*Y(9)*Y(12)+K18*Y(16)*M+K24*Y(14)*M+K28*Y(14)*Y(12)+K38*Y(12)*Y(13)+K40*Y(3)*Y(12)+K41*Y(4)*Y(13)+K42*Y(13)*Y(3)+K48*Y(10)*Y(13)-K1*Y(1)*Y(11)-K4*Y(7)*Y(11)-K10*Y(8)*Y(11)-K21*Y(15)*Y(11)-K26*Y(14)*Y(11)-K32*Y(13)*Y(11)*M-K33*Y(12)*Y(11)*M-2*K35*Y(11)*Y(11)*M-K36*Y(11)*Y(6)*M-K39*Y(11)*Y(13)-K43*Y(17)*Y(11)-K44*Y(17)*Y(11); f12=K05*Y(2)*Ne+K16*Y(9)*Y(6)+K37*Y(13)*Y(13)+K39*Y(11)*Y(13)+K47*Y(10)*Y(6)-K2*Y(1)*Y(12)-K6*Y(7)*Y(12)-K17*Y(9)*Y(12)-K20*Y(15)*Y(12)-K27*Y(14)*Y(12)-K28*Y(14)*Y(12)-K33*Y(12)*Y(11)*M-2*K34*Y(12)*Y(12)*M-K38*Y(12)*Y(13)-K40*Y(3)*Y(12)-K45*Y(17)*Y(12); f13=K2*Y(1)*Y(12)+K7*Y(7)*Y(6)+K8*Y(7)*Y(17)+K11*Y(8)*Y(6)+K12*Y(8)*Y(6)+K15*Y(9)*Y(6)+K20*Y(15)*Y(12)+K27*Y(14)*Y(12)+K33*Y(12)*Y(11)*M+K40*Y(3)*Y(12)+2*K43*Y(17)*Y(11)+K45*Y(17)*Y(12)-K3*Y(1)*Y(13)-K5*Y(7)*Y(13)-K14*Y(8)*Y(13)-K22*Y(15)*Y(13)-K25*Y(14)*Y(13)-K32*Y(13)*Y(11)*M-2*K37*Y(13)*Y(13)-K38*Y(12)*Y(13)-K39*Y(11)*Y(13)-K41*Y(4)*Y(13)-K42*Y(13)*Y(3)-K46*Y(17)*Y(13)-K48*Y(10)*Y(13); f14=K11*Y(8)*Y(6)+K16*Y(9)*Y(6)+K20*Y(15)*Y(12)+K21*Y(15)*Y(11)+K22*Y(15)*Y(13)+K23*Y(15)*Y(7)+K49*Y(2)*Y(9)-K24*Y(14)*M-K25*Y(14)*Y(13)-K26*Y(14)*Y(11)-K27*Y(14)*Y(12)-K28*Y(14)*Y(12)-K29*Y(14)*Y(6)-K30*Y(14)*Y(17)-K31*Y(14)*Y(7); f15=K5*Y(7)*Y(13)+K6*Y(7)*Y(12)+K7*Y(7)*Y(6)+K14*Y(8)*Y(13)+K18*Y(16)*M+K19*Y(16)*Y(6)+K30*Y(14)*Y(17)-K20*Y(15)*Y(12)-K21*Y(15)*Y(11)-K22*Y(15)*Y(13)-K23*Y(15)*Y(7); f16=K8*Y(7)*Y(17)-K18*Y(16)*M-K19*Y(16)*Y(6); f17=K19*Y(16)*Y(6)+K29*Y(14)*Y(6)+K36*Y(11)*Y(6)*M-K8*Y(7)*Y(17)-K9*Y(7)*Y(17)-K30*Y(14)*Y(17)-K43*Y(17)*Y(11)-K44*Y(17)*Y(11)-K45*Y(17)*Y(12)-K46*Y(17)*Y(13); dYdt=[f1;f2;f3;f4;f5;f6;f7;f8;f9;f10;f11;f12;f13;f14;f15;f16;f17];
% ------------------------------------------------------------------------- function K0 = Krate(A,B,C,T) NA=6.02217e+23; % T: Temperatre, K K0 = A*T^B*exp(-C/T)*NA*1e-06; 
% ------------------------------------------------------------------------- function Kco2 = func(Eev,Eth,Te) sigma=1e-16;Rydberg=1.361e-2; E1=Eev-Eth; % Eev和Eth单位为keV c1=3.160e-1;c2=1.046;c3=1.210e-2;c4=1.320e-1;c5=5.140e-2;c6=9.300e-1; f1=sigma*c1*((E1./Rydberg).^c2); f3=f1./(1+((E1./c3).^(c2+c4))+((E1./c5).^(c2+c6))); Kco2=1e+6*Eev.*f3.*exp(-1000*Eev./Te); % ------------------------------------------------------------------------- function K1= rate1(A,b,E,T,n) % T: Temperatre, K j=4.1840;R=8.31434; K1 = A*T^b*exp(-E*j/(R*T))*(1e-6)^(n-1); % ------------------------------------------------------------------- function K2= rate2(A,b,E,T,n) % T: Temperatre, K R=8.31434; K2 = A*T^b*exp(-E*1000/(R*T))*(1e-6)^(n-1); % ------------------------------------------------------------------- function K3= rate3(A,B,EdR,T,n) % T: Temperatre, K K3 = A*T^B*exp(-EdR/T)*(1e-6)^(n-1);

[此贴子已经被作者于2005-12-24 13:54:52编辑过]

xueyi · 发表于 2005-12-24 14:37

不爱看程序，不过给个建议， 看下帮助是不是2个ode实现过程的要求不一样 也即其收敛条件不一样

aspen · 发表于 2005-12-25 08:42

估计你的方程是刚性的，ode45适合求解非刚性方程，ode23s适合求解刚性方程

账号		自动登录	找回密码
密码			我要加入

[求助]matlab运行了10个小时还不出结果

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

回复：（mathtype）[求助]matlab运行了10个小时还不出...