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.