[R-sig-dyn-mod] NA error in deSolve
Dabing Chen
dabing.c at gmail.com
Fri Jan 16 22:36:26 CET 2015
Hi All:
On Thu, Jan 15, 2015 at 11:20 AM, Dabing Chen <dabing.c at gmail.com> wrote:
> 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]]
More information about the R-sig-dynamic-models
mailing list