[R-sig-dyn-mod] NA error in deSolve
Dabing Chen
dabing.c at gmail.com
Fri Jan 16 22:43:00 CET 2015
Hi All:
I figured that the numbers were not updated during the
integration, so I added a few lines inside the integration. The P_diss
function seems to work, but I got another issue:
Warning messages:
1: In lsoda(y, times, func, parms, ...) :
repeated convergence test failures on a step, but integration was
successful - inaccurate Jacobian matrix?
2: In lsoda(y, times, func, parms, ...) :
Returning early. Results are accurate, as far as they go
There are NAs in the output, How do I remove the NAs? Thanks,
rm(list=ls())
library (deSolve)
W_frac <- rep(0.2,5)# the weight fraction
R_size <- seq(10,50,10) # the weight fraction corresponding size in um
Dose <- 50 # initial dose in mg
Particle <- rep(0,5)
Size <- R_size
Y <- 0
Tablet <- Dose
state <- c(Tablet,Particle,Y,Size)
S = 0.1 # 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)), {
Tablet <- state[1]
Particle <- state[2:6]
Y <- state[7]
Size<- state[8:12]
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 Fri, Jan 16, 2015 at 4:36 PM, Dabing Chen <dabing.c at gmail.com> wrote:
> 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