Why ode45 can not solve while ode15s can solve?

5 次查看（过去 30 天）
Jidapa Adam
Jidapa Adam 2022-3-18
评论： Steven Lord 2022-3-18
My code can solve with ode15s but I would like to try on ode45, So can I use RelTol with value 1e-8 or can you give me other value of RelTol?
P0 = [0.5] ;%bar (Intial pressure)
Ta0 = [873] ;%K (Intial air heat temperature)
repeat = [1] ;%Repeat condition Ta constant and Ta non constant
sheet = 0 ;%Sheet in excel
condition = 0 ;%Number of condition
for a = 1:length(repeat)
for i = 1:length(P0)
for j = 1:length(Ta0)
T0 = Ta0(j) ;%K (Intial temperature before entering bed)
pj0 = yj0.*P0(i) ;%bar (Partial pressure)
Fj0 = (pj0.*v0_feed)./(R_feed.*T0) ;%mol/s (Molar flowrate)
Lspan = [0:0.01:0.03] ;%m (Bed length range)
dep_var0 = [Fj0' T0 Ta0(j)] ;%[mol/s,K,K] (Intial dependent variable)
condition = condition + 1 ;
var = [P0(i) Ta0(j) condition T0] ;%[bar,K,-,K] (Another variable that nessery)
options = odeset('RelTol',1e-3);
[L ,dep_var] = ode45(@func_NiCo,Lspan,dep_var0,options,var) ;
  2 个评论
Walter Roberson
Walter Roberson 2022-3-18
What happened when you tried with ode45?
Note: in order for us to test we would need your func_NiCo as well
Jidapa Adam
Jidapa Adam 2022-3-18
When I run ode45 it doesn't stop running and when I change RelTol value to 1e-1 the result is complex number.
load Data_import_NiCo
P0 = [0.5] ;%bar (Intial pressure)
Ta0 = [873] ;%K (Intial air heat temperature)
repeat = [1] ;%Repeat condition Ta constant and Ta non constant
sheet = 0 ;%Sheet in excel
condition = 0 ;%Number of condition
for a = 1:length(repeat)
for i = 1:length(P0)
for j = 1:length(Ta0)
T0 = Ta0(j) ;%K (Intial temperature before entering bed)
pj0 = yj0.*P0(i) ;%bar (Partial pressure)
Fj0 = (pj0.*v0_feed)./(R_feed.*T0) ;%mol/s (Molar flowrate)
Lspan = [0:0.01:0.03] ;%m (Bed length range)
dep_var0 = [Fj0' T0 Ta0(j)] ;%[mol/s,K,K] (Intial dependent variable)
condition = condition + 1 ;
var = [P0(i) Ta0(j) condition T0] ;%[bar,K,-,K] (Another variable that nessery)
options = odeset('RelTol',1e-3);
[L ,dep_var] = ode45(@func_NiCo,Lspan,dep_var0,options,var) ;
this is func_NiCo
function dydx = func_NiCo(y,x,op)
load Data_import_NiCo
FCH4 = x(1,1) ;
FCO2 = x(2,1) ;
FCO = x(3,1) ;
FH2 = x(4,1) ;
FH2O = x(5,1) ;
FC = x(6,1) ;
FN2 = x(7,1) ;
T = x(8,1) ;
Ta = x(9,1) ;
P0 = op(1) ;%bar (Intial pressure)
Tac = op(2) ;%K (Intial air heat temperature constant)
condition = op(3) ;
T0 = op(4) ;%K (Intial temperature)
FT = FCH4+FCO2+FCO+FH2+FH2O+FC+FN2 ;%mol/s (Total molar flowrate)
%% Partial pressure (bar)
pCH4 = P0*(FCH4/FT)*(T0/T) ;
pCO2 = P0*(FCO2/FT)*(T0/T) ;
pCO = P0*(FCO/FT)*(T0/T) ;
pH2 = P0*(FH2/FT)*(T0/T) ;
pH2O = P0*(FH2O/FT)*(T0/T) ;
pC = P0*(FC/FT)*(T0/T) ;
pN2 = P0*(FN2/FT)*(T0/T) ;
%% Heat capacity (J/(mol*K))
cpj = heatcp(:,1) + heatcp(:,2)*T + heatcp(:,3)*(T^2) + heatcp(:,4)*(T^-2) ;
%% Del_Heat of reaction
delH1 = hrxi_Tr(1)+((2.*cpj(3)+2*cpj(4)-cpj(1)-cpj(2)).*(T-Tr)) ;
delH2 = hrxi_Tr(2)+((cpj(3)+cpj(5)-cpj(2)-cpj(4)).*(T-Tr)) ;
delH3 = hrxi_Tr(3)+((cpj(6)+2.*cpj(4)-cpj(1)).*(T-Tr)) ;
delH4 = hrxi_Tr(4)+((cpj(3)+cpj(4)-cpj(5)-cpj(6)).*(T-Tr)) ;
delH5 = hrxi_Tr(5)+((2*cpj(3)-cpj(6)-cpj(2)).*(T-Tr)) ;
%% Reaction rate constant (mol/(kg*s))
ki = k0i.*exp(-Eai./T) ;
%% Equilibrium constant
Kpi = Kp0i.*exp(-dHrxi_Tr./T) ;%[bar^2,-,bar,bar,bar]
%% Adsorption equilibrium constant
KCH4_1 = K0j1(1)*exp(hdx1(1)/T) ;%mol/(kg*s*bar)
KCO2_1 = K0j1(2)*exp(hdx1(2)/T) ;%mol/(kg*s*bar)
KCO2_2 = K0j2(2)*exp(hdx2(2)/T) ;%bar^-1
KH2_2 = K0j2(4)*exp(hdx2(4)/T) ;%bar^-1
KCH4_3 = K0j3(1)*exp(hdx3(1)/T) ;%bar^-1
KH2_3 = K0j3(4)*exp(hdx3(4)/T) ;%bar^1.5
KCH4_4 = K0j4(1)*exp(hdx4(1)/T) ;%bar^-1
KH2_4 = K0j4(4)*exp(hdx4(4)/T) ;%bar^1.5
KH2O_4 = K0j4(5)*exp(hdx4(5)/T) ;%-
KCO2_5 = K0j5(2)*exp(hdx5(2)/T) ;%bar
KCO_5 = K0j5(3)*exp(hdx5(3)/T) ;%bar^-1
%% Rate reaction (mol/(kg*s)
r1 = ((ki(1)*pCH4*pCO2)/(KCO2_1*pCO2+KCH4_1*pCH4))*(1-(((pCO*pH2)^2)/(Kpi(1)*pCH4*pCO2))) ;
r2 = ((ki(2)*KCO2_2*KH2_2*pCO2*pH2)/((1+KCO2_2*pCO2+KH2_2*pH2)^2))*(1-((pCO*pH2O)/(Kpi(2)*pCO2*pH2))) ;
r3 = ((ki(3)*KCH4_3)*(pCH4-((pH2^2)/Kpi(3))))/((1+KCH4_3*pCH4+(pH2^1.5/KH2_3))^2) ;
r4 = ((ki(4)/KH2O_4)*((pH2O/pH2)-(pCO/Kpi(4))))/((1+KCH4_4*pCH4+(pH2O/(KH2O_4*pH2))+(pH2^1.5/KH2_4))^2) ;
r5 = ((ki(5)/(KCO_5*KCO2_5))*((pCO2/pCO)-(pCO/Kpi(5))))/((1+KCO_5*pCO+(pCO2/(KCO_5*KCO2_5*pCO)))^2) ;
%% rnet of species (mol/(kg*s))
rn_CH4 = -r1-r3 ;
rn_CO2 = -r1-r2-r5 ;
rn_CO = 2*r1+r2+r4+2*r5 ;
rn_H2 = 2*r1-r2+2*r3+r4 ;
rn_H2O = r2-r4 ;
rn_C = r3-r4-r5 ;
rn_N2 = 0 ;
%%
sumFjcpj = FCH4*cpj(1)+ FCO2*cpj(2)+ FCO*cpj(3)+ FH2*cpj(4)+ FH2O*cpj(5)+ FC*cpj(6)+ FN2*cpj(7) ;%Watt/K
sumrihi = r1*-delH1+ r2*-delH2+ r3*-delH3+ r4*-delH4+ r5*-delH5 ;%Watt/kg
%% ODE equations
dFCH4dL = Rho_bed*Ac*rn_CH4 ;%mol/(m*s)
dFCO2dL = Rho_bed*Ac*rn_CO2 ;%mol/(m*s)
dFCOdL = Rho_bed*Ac*rn_CO ;%mol/(m*s)
dFH2dL = Rho_bed*Ac*rn_H2 ;%mol/(m*s)
dFH2OdL = Rho_bed*Ac*rn_H2O ;%mol/(m*s)
dFCdL = 0 ;%mol/(m*s)
dFN2dL = Rho_bed*Ac*rn_N2 ;%mol/(m*s)
if condition <= 6
dTdL = ((Rho_bed*Ac*sumrihi)-(U*ain*Ac*(T-Tac)))/sumFjcpj ;%K/m
dTadL = 0 ;%K/m
else
dTdL = ((Rho_bed*Ac*sumrihi)-(U*ain*Ac*(T-Ta)))/sumFjcpj ;%K/m
dTadL = (U*aout*Ac_air*(T-Ta))/(mc_air*cpj(8)) ;%K/m
end
%% Calculations
dydx(1,1) = dFCH4dL ;
dydx(2,1) = dFCO2dL ;
dydx(3,1) = dFCOdL ;
dydx(4,1) = dFH2dL ;
dydx(5,1) = dFH2OdL ;
dydx(6,1) = dFCdL ;
dydx(7,1) = dFN2dL ;
dydx(8,1) = dTdL ;
dydx(9,1) = dTadL ;
end

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2022-3-18
You have a "stiff" system of equations. The equations wobble a lot and ode45 has to take very small steps to try to follow the wobbles because they might be important.
  4 个评论
显示 2更早的评论
Walter Roberson
Walter Roberson 2022-3-18
编辑：Walter Roberson 2022-3-18
You can fix the problem by stopping interrupting ode45 and let it run the several years it needs to calculate the system properly.
(It is not a bug.)
Steven Lord
Steven Lord 2022-3-18
If you're not familiar with stiffness in ODEs I recommend reading through this post on Cleve Moler's blog. Run at least the one example Cleve recommends running and read his explanation for its performance.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by