See inline remarks.
On 16-07-2013, at 10:07, Rapha?lle Carraud <raphaelle.carraud at
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 at 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.