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