How can I check which parameters is invalid in my code?

1 次查看（过去 30 天）
Dursman Mchabe
Dursman Mchabe 2018-9-24
Hi everyone, I am trying to run an ODE15i code. I get an error message:
"
Error using mupadengine/feval (line 187)
Invalid parameter '1770887431076117/1180591620717411303424'.
Error in sym/daeFunction (line 117)
A = feval(symengine, 'daetools::daeFunction', eqs, vars, params{:});
Error in SlurryAbsorptionModel1 (line 60)
f = daeFunction(eqns,vars, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, kga, DCa2, DSO2, DHSO3, DSO32, kLa_SO2,
kLa_CO2, HCO2, DCO2, DHCO3, DCO32, KSPCaSO3, BETCaSO3, ktot, BETCaCO3, MWCaCO3, Kad, KCO2, KHCO3, KSO2, KHSO3, Kw);
"
Debugging points me here:
"
f = daeFunction(eqns,vars, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, kga, DCa2, DSO2, DHSO3, DSO32, kLa_SO2, kLa_CO2, HCO2, DCO2, DHCO3, DCO32, KSPCaSO3, BETCaSO3, ktot, BETCaCO3, MWCaCO3, Kad, KCO2, KHCO3, KSO2, KHSO3, Kw);
" The mupadengine.m points me to: "
[stmt,fcn2,args2] = feval2eval(fcn,varargin); %#ok extra args are for holding sym references
[S,err] = evalin(engine,stmt,'message');
if nargout < 2 && err ~= 0
if isa(S,'message') && contains(S.Identifier,':Singularity')
S = sym(NaN);
elseif isa(S,'message')
error(S);
else
error(message('symbolic:mupadengine:feval:FevalError',S));
end
end
if isa(S,'message')
S = getString(S);
end
end
"
To me, all the parameters are valid. How can I check which one is not valid?
The complete code is given below:
syms CSO2_gas(t) CCO2_gas(t) S_total(t) C_total(t) Ca2_total(t) CCaCO3(t) CCaSO3(t) CH(t) V_Headspace F_rate CSO2_in CCO2_in R T HSO2 kga DCa2 DSO2 DHSO3 DSO32 kLa_SO2 kLa_CO2 HCO2 DCO2 DHCO3 DCO32 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad KCO2 KHCO3 KSO2 KHSO3 Kw CCa2
global V_Headspace F_rate CSO2_in CCO2_in CCa2 R T HSO2 kga DCa2 DSO2 DHSO3 DSO32 kLa_SO2 kLa_CO2 HCO2 DCO2 DHCO3 DCO32 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad KCO2 KHCO3 KSO2 KHSO3 Kw
V_Headspace = 1.5e-6;
F_rate = 1.66667e-5;
CSO2_in = 6.51332e-2;
CCO2_in = 0;
CCa2 = 10;
R = 8.314;
T = 323.15;
HSO2 = 149;
kga = 4.14e-6;
DCa2 = 1.39e-9;
DSO2 = 2.89e-9;
kLa_SO2 = 8.4e-4;
kLa_CO2 = 9.598e-4;
HCO2 = 5.15e3;
DCO2 = 3.53e-9;
KSPCaSO3 = 1.07e-7;
BETCaSO3 = 10;
ktot = 8.825e-3;
BETCaCO3 = 12.54;
MWCaCO3 = 100.0869;
Kad = 0.84;
KCO2 = 1.7e-3;
KHCO3 = 6.55e-8;
KSO2 = 6.24;
KHSO3 = 5.68e-5;
Kw = 5.3e-8;
eqn1 = diff(CSO2_gas(t),t)== F_rate/ V_Headspace * ( CSO2_in - CSO2_gas(t)) - (CSO2_gas(t) * R * T - HSO2 * ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / ((1/kga) + 1/(kLa_SO2 * (((DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) + DHSO3 * (((KSO2 * HSO2 * CSO2_gas(t) * R * T/CH(t)) - ((S_total(t) * KSO2 * CH(t))/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 *HSO2* CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3))))))));
eqn2 = diff(CCO2_gas(t) ,t)== F_rate/ V_Headspace * ( CCO2_in - CCO2_gas(t)) - (kLa_CO2 * ((((DCO2*(HCO2 * CCO2_gas(t) * R * T - ((C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3))) + DHCO3*(((KCO2 * HCO2 * CCO2_gas(t) * R * T/CH(t)) - ((C_total (t) * KCO2 * CH(t))/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * CCO2_gas(t) * R * T/CH(t)^2) - ((C_total(t) * KCO2 * KHCO3)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * CCO2_gas(t) * R * T - ((C_total (t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3)))))))*(CCO2_gas(t) * R * T/HCO2 - ((C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3))));
eqn3 = diff(S_total(t),t)== (CSO2_gas(t) * R * T - HSO2*((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3))) / ((1/kga) + 1/(kLa_SO2 * (((DSO2*(HSO2* CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3))) + DHSO3*(((KSO2 * HSO2 * CSO2_gas(t) * R * T/CH(t)) - ((S_total(t) * KSO2 * CH(t))/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 *HSO2* CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))) - 0.162 * exp(-5153/T) * (((CCa2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3))) / KSPCaSO3) - 1)^3 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn4 = diff(C_total(t),t)== (kLa_CO2 * ((((DCO2*(HCO2 * CCO2_gas(t) * R * T - ((C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3))) + DHCO3*(((KCO2 * HCO2 * CCO2_gas(t) * R * T/CH(t)) - ((C_total (t) * KCO2 * CH(t))/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * CCO2_gas(t) * R * T/CH(t)^2) - ((C_total(t) * KCO2 * KHCO3)/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3)))))) / (DCO2 * (HCO2 * CCO2_gas(t) * R * T - ((C_total (t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3)))))))*(CCO2_gas(t) * R * T /HCO2 - ((C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3)))) + (ktot * BETCaCO3 * MWCaCO3 * CCaCO3(t) * CH(t) *(1 - (Kad * CH(t))/(1 + Kad * CH(t))));
eqn5 = diff(Ca2_total(t),t)== (-1) * (ktot * BETCaCO3 * MWCaCO3 * CCaCO3(t) * CH(t) *(1 - (Kad * CH(t))/(1 + Kad * CH(t)))) - (0.162 * exp(-5153/T) * (((CCa2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3))) / KSPCaSO3) - 1)^3 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn6 = diff(CCaCO3(t),t)== (-1) * (ktot * BETCaCO3 * MWCaCO3 * CCaCO3(t) * CH(t) * (1 - (Kad * CH(t))/(1 + Kad * CH(t))));
eqn7 = diff(CCaSO3(t),t)== 0.162 * exp(-5153/T) * ((((CCa2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^3 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn8 = diff(CH(t),t)== CH(t) + 2 * CCa2 - ((S_total(t) *KSO2 * CH(t))/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))- 2*((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)) - ((C_total(t) * KCO2 * CH(t))/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)) - 2 * ((C_total(t) * KCO2 * KHCO3)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))- Kw / CH(t);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8];
vars = [CSO2_gas(t); CCO2_gas(t); S_total(t); C_total(t); Ca2_total(t); CCaCO3(t); CCaSO3(t); CH(t)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, kga, DCa2, DSO2, DHSO3, DSO32, kLa_SO2, kLa_CO2, HCO2, DCO2, DHCO3, DCO32, KSPCaSO3, BETCaSO3, ktot, BETCaCO3, MWCaCO3, Kad, KCO2, KHCO3, KSO2, KHSO3, Kw);
F = @(t, Y, YP) f(t, Y, YP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, kga, DCa2, DSO2, DHSO3, DSO32, kLa_SO2, kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, ktot, BETCaCO3, MWCaCO3, Kad, KCO2, KHCO3, KSO2, KHSO3, Kw);
vars;
y0est = [0; 0; 0; 0; 0; 0; 0; 0];
yp0est = zeros(8,1);
opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [0, 10], y0, yp0, opt);
for k = 1:origVars
S{k} = char(vars(k));
end

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

回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by