[R-sig-dyn-mod] NA error in deSolve

[Dabing Chen]

Hi Thomas:
           Thanks a lot for your help. You are right, the integration
stopped at the first step, where the X0 was 0 and the condition was NA.
           To avoid the problem, I included the particle size map in the
code. However, the code itself did not give me any error message, but the
P_diss function did not work as I expected. There must be something wrong,
I could not figure it our.
           Best regards!

Dabing


rm(list=ls())
library (deSolve)


      W_frac <- rep(0.2,5)# the weight fraction
      R_size <- seq(10,50,10) # the corresponding size in um

Dose <- 200 # initial dose in mg



Particle <- rep(0,5)

Size <- R_size

Y <- 0

Tablet = Dose

state <- c(Tablet,Particle,Y,Size)



S = 6 # solubility in mg/ml
den = 1300 # density unit mg/ml
Diff = 7E-6   # diffusion coeffiecient cm2/s
V= 900 # volume of dissolution medium


P_diss <- function(X,X0,r){

  thick <- r
  thick[thick > 30] <- 30
  disso <- 3*60*Diff*(X0)^(1/3)*X^(2/3)/(den*thick/10000*r/10000)*(S-Y/V)
  result <- -pmin(disso,X)

  result[(X<0)|(S<Y/V)]<- 0

  return(result)
}



P_size <- function(X,X0){


    result <- (X/X0)^(1/3)
    result[X0<=0] <- 0
    return(result)
}


Lorenz <- function(t, state, para){
with (as.list(c(state, para)), {


           ddisinte <- 0.02  # rate of tablet disintegration

           dTablet <- -ddisinte*Tablet
           dParticle <- ddisinte*Tablet*W_frac +
P_diss(Particle,(Dose-Tablet)*W_frac,Size)
           dY <- -sum(P_diss(Particle,(Dose-Tablet)*W_frac,Size))
           dSize <- -P_size(dParticle,(Dose-Tablet)*W_frac)*Size


list (c(dTablet, dParticle, dY,dSize))

})
}

times <- seq(0,1440, by = 1)

out <- ode(y = state, t = times, func = Lorenz, parms = 0)

write.csv(out, file = "out.csv")








On Wed, Jan 14, 2015 at 6:15 PM, Thomas Petzoldt <
thomas.petzoldt at tu-dresden.de> wrote:

> Hi,
>
> this has nothing to do with deSolve, it is your function P_diss, that
> needs debugging. Enter some print or cat functions like this:
>
> P_diss <- function(X,X0,r){
>   result <- -pmin(-3*Diff*X0^(2/3)*X^(1/3)/(den*(30E-4)*r)*(S-Y/V)*60,X)
>   result[(X<=0)|(S<Y/V)]<- 0
>   cat("X=", X, "X0=", X0, "\n")
>   cond1 <- r*(X/X0)^(1/3) < 30E-4
>   cat(cond1, "\n")
>   result[cond1] <-
>     -pmin(-3*Diff*X0^(1/3)*X^(2/3)/(den*r^2)*(S-Y/V)*60,X)[cond1]
>   return(result)
> }
>
> ...  and you'll see:
>
> 0/0 is NaN
>
> and NaN < 30E-4 is NA
>
> result[NA] is an error.
>
>
> Hope it helps
>
> thpe
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>

	[[alternative HTML version deleted]]