[R] Ode error message

[Berend Hasselman]

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.