See inline remarks.
On 16-07-2013, at 10:07, Raphaëlle Carraud raphaelle.carr...@oc-metalchem.com
wrote:
Hello,
I am creating a program with R to solve a differential equation system.
However, I get the following message I do not understand :
out - ode(y = state, times = z, func = liquide, parms = 0, atol = 0)
DLSODA- EWT(I1) is R1 .le. 0.0
In above message, I1 = 1
In above message, R1 = 0
Error in lsoda(y, times, func, parms, ...) :
illegal input detected before taking any integration steps - see written
message
or this one when I tried modifying the atoll value :
out - ode(y = state, times = z, func = liquide, parms = 0, atol = 10^-14)
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R1 = 0, R2 = 0
DINTDY- T (=R1) illegal
In above message, R1 = 1
T not in interval TCUR - HU (= R1) to TCUR (=R2)
In above message, R1 = 0, R2 = 0
DINTDY- T (=R1) illegal
In above message, R1 = 2
T not in interval TCUR - HU (= R1) to TCUR (=R2)
In above message, R1 = 0, R2 = 0
DLSODA- Trouble in DINTDY. ITASK = I1, TOUT = R1
In above message, I1 = 1
In above message, R1 = 2
Here is my program. I also tried changing the initial values but it does not
work well.
liquide - function(z, state, parameters) {
with(as.list(c(state,parameters)),{
# rate of change
Tr - 273+90
Why are you defining Tr? It is not used anywhere
C - CA + CB + CC + CD + CE + CI + CG + CJ + CK + CH
Same thing. Not used.
K32 - 6.54*10^4
K33 - 1.37*10^4
K34 - 330
K35 - 5.81*10^4
kf2 - 1.37*10^3
kf3 - 1.37*10^3
kf4 - 8.68*10^5
kf5 - 157.2
K2 - 10^1.37
K3 - 10^(-3.35)
r1 - kf4*CD - kf4/K34*CE^2
r2 - kf3*CA*CB - kf3/K33*CD
r3 - kf2*CA^2 - kf2/K32*CC
r4 - kf5*CC - kf5/K35*CE*CI^2
dCA - -r2 # dNO/dt
dCB - -r3 - r2 # dNO2/dt
dCC - r3/2 - r4 # dN2O4/dt
dCD - r2 - r1# dN2O3/dt
dCE - 2*r1 + r4# dHNO2/dt
dCI - r4 #
dHNO3/dt
dCG - -r4 - r1 # dH2O/dz
dCH - (dCE + dCI)/((K2 + K3)*(CE + CI)) # dH/dz
dCJ - (CH*dCI - CI*dCH)/(K3*CH^2) # dNO3-/dz
You are dividing by CH, which is 0 initially. So what value does dCH then get?
dCK - (CH*dCE - CE*dCH)/(K2*CH^2)# dNO2-/dz
Same thing.
list(c(dCA, dCB, dCC, dCD, dCE, dCI, dCG, dCH, dCJ, dCK))
}) # end with(as.list ...
}
Ti - 273+90 # K
You are not using Ti.
Ct - 5100 # mol/m^3
And Ct is also not used.
state -c(CA = 0, # mol/m^3 NO2
CB = 0, # mol/m^3 NO
CC = 0,# mol/m^3 N2O4
CD = 0, # mol/m^3 N2O3
CE = 50, # mol/m^3 HNO2
CI = 50, # mol/m^3 HNO3
CG = 5000, # mol/m^3 H2O
CH = 0, # mol/m^3 H+
0!!
CJ = 0,# mol/m^3 NO3-
CK = 0) # mol/m^3 NO2-
parameters - c(Ct = 5100)
Why not parameters - c(Ct = Ct)?
z - seq(0, 15, by = 1) # en seconde
library(deSolve)
out - ode(y = state, times = z, func = liquide, parms = 0, atol = 10^-14)
head(out)
plot(out)
You will still get messages.
You should really learn to de elementary debugging.
Such as inserting
liquide(0,state,parameters)
after defining state and parameters to check and test.
Berend
Thank you
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
