Hi All:
     The error message at ka =0.01 is like this:

1: In lsoda(y, times, func, parms = parms, bandup = bandup, banddown =
banddown,  :
  repeated convergence test failures on a step, but integration was
successful - inaccurate Jacobian matrix?
2: In lsoda(y, times, func, parms = parms, bandup = bandup, banddown =
banddown,  :
  Returning early. Results are accurate, as far as they go

     Best regards!
Dabing

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

ka <- 0.01

N= 100

Dose <- 100 # initial dose in mg

V1 <- 300
V2<-  c(rep(3,30),rep(0.5,(N-30)))


Si = 0.02 # initial solubility in jejunum unit is mg/ml
Sf = 0.02 # final solubility in jejunum is mg/ml
t0 = 20 # time for half phase transformation in jejunum in min
p0 = 0.5

SI <- rep(0.02,N)
GImobility = N/60/8

for (i in 1:N){
  SI[i] <- (Si-Sf)/(1+exp(p0*((i/GImobility)-t0)))+Sf
}

Weibulla <- 2 # weibull function for time factor
Weibullb <-0.9 # weibull fuction for shape factor
Tlag <- 0 # lag time for disintegration
criticT <- (15)^(1/Weibullb)*Weibulla # critical T for disintegration

z <- 50

para <- c(SGi <- 0.8, # inital solubility in Stomach unit is mg/ml
          SGf <- 0.8,# final solubility in Stomach unit is mg/ml
          TG <- 20, #time for half transformation in Stomach
          PG <- 0.2) # rate of half transformation in Stomach

 advModel <- function(t, state, parms) {

     with (as.list(parms), {
         Stomach1 <- state[1]  # undisintegrated
         Stomach2 <- state[2]   # granule
         Stomach3 <- state[3]   # Precipitate
         Stomach4 <- state[4]   # dissolved
         Intestine1 <- state[5:(N+4)]           # undisintegated
         Intestine2 <- state[(N+5):(2*N+4)]     # granule
         Intestine3 <- state[(2*N+5):(3*N+4)]   # precipitate
         Intestine4 <- state[(3*N+5):(4*N+4)]   # dissolved
         T1 <- state[4*N+5]   # total undisintegrated mass
         T2 <- state[4*N+6]   # total granule
         T3 <- state[4*N+7]   # total precipitate
         T4 <- state[4*N+8]   # total in solution

         absorption <- state[4*N+9]
         plasma <- state[4*N+10]

         SG <- SGi

         ddisinte <- 0.03
         # wettability factor

         # define a function for dissolution kinetics
         P_diss <- function(X,Y,V,S){
           # X the solid remaining at time t, X0, the inital dose, r, the
particle size
           if (X <= 0) {
             result <- 0  # avoid numerical error
           } else if (S < Y/V){
             result <- 0
           } else {
             result <- min(z*(X^(1/3))*(S-Y/V),X)
           }
           return(result)
         }

         # define the precipitation rate function

         P_precip <- function(X,Y,V,S){
           # X the solid remaining at time t, assuming the precipitate has
particle size of 1 um
           if (X <= 0) {
             result <- 0  # avoid numerical error
           } else if (Y/V > S){

           result <- V * ( Y/V-S)
           } else {
             result <- -min(z*(X^(1/3))*(S-Y/V),X) # +  (V*( Y/V-S) +
min(z*(X^(1/3))*(S-Y/V),X))/(1+exp(1000*(S-Y/V)))
           }
           return(result)
         }


         emptying <-
0.01+0.05*(tanh(10*(t%%240-210))-tanh(10*(t%%240-220)))
         emptying1 <- 0.05*(tanh(10*(t%%240-210))-tanh(10*(t%%240-220)))


         # dissolution in Stomach
         disso1 <- P_diss(Stomach2,Stomach4,V1,SG)
         precip1 <- P_precip(Stomach3,Stomach4,V1,SG)

         dStomach1 <- -ddisinte*Stomach1 - emptying1*Stomach1
         dStomach2 <- ddisinte*Stomach1 - disso1 - emptying*Stomach2
         dStomach3 <- precip1 -emptying*Stomach3
         dStomach4 <- disso1 - precip1 - emptying*Stomach4

         disso2 <- rep(0,N)
         precip2 <- rep(0,N)

         for (i in 1:N){
           disso2[i] <- P_diss(Intestine2[i],Intestine4[i],V2[i],SI[i])
           precip2[i] <- P_precip(Intestine3[i],Intestine4[i],V2[i],SI[i])
         }

         Tran1 <- advection.1D(C = Intestine1, flux.up =
emptying1*Stomach1, v = GImobility, dx=1,adv.method="super")
         Tran2 <- advection.1D(C = Intestine2, flux.up = emptying*Stomach2,
v = GImobility, dx=1,adv.method="super")
         Tran3 <- advection.1D(C = Intestine3, flux.up = emptying*Stomach3,
v = GImobility, dx=1,adv.method="super")
         Tran4 <- advection.1D(C = Intestine4, flux.up = emptying*Stomach4,
v = GImobility, dx=1,adv.method="super")


        dIntestine1 <- Tran1$dC - ddisinte*Intestine1
        dIntestine2 <- Tran2$dC - disso2 + ddisinte*Intestine1
        dIntestine3 <- Tran3$dC + precip2
        dIntestine4 <- Tran4$dC + disso2 - precip2 - ka*Intestine4

         dT1 <- sum(dIntestine1)
         dT2 <- sum(dIntestine2)
         dT3 <- sum(dIntestine3)
         dT4 <- sum(dIntestine4)


        dabsorption <- sum(ka*Intestine4)
        dplasma <- dabsorption -0.002*plasma


         return (list(c(dStomach1,dStomach2,dStomach3,dStomach4,
                        dIntestine1,dIntestine2,dIntestine3,dIntestine4,
                        dT1,dT2,dT3,dT4,
                        dabsorption,dplasma)))
       })
}

times <- seq(1,1000,by =1)

state <- c(Dose,rep(0,3),rep(0,4*N),rep(0,4),0,0)

systime <- system.time(out <- ode.1D(func = advModel, y=state,times =times,
                                     parms = para,nspec=1))



On Thu, Jul 25, 2013 at 3:36 PM, Dabing Chen <dabing.c@gmail.com> wrote:

> Hi All:
>
>          I am having some problems with the program below. I was able to
> run when ka = 0.04. However, when I change it to 0.01, the program cannot
> finish.
>
>          I also tried using different integration methods in ode.1D, such
> as "vode, daspk", which did not work. radau runs almost 3 min, the program
> runs better when I left it blank, which took about 1 min to run.
>
>          How can I make it faster and more robust?
>          Thanks a lot in advance.
>
> Dabing
>
>
> rm(list=ls())
> library (deSolve)
> library(ReacTran)
>
> ka <- 0.04
>
> N= 100
>
> Dose <- 100 # initial dose in mg
>
> V1 <- 300
> V2<-  c(rep(3,30),rep(0.5,(N-30)))
>
>
> Si = 0.02 # initial solubility in jejunum unit is mg/ml
> Sf = 0.02 # final solubility in jejunum is mg/ml
> t0 = 20 # time for half phase transformation in jejunum in min
> p0 = 0.5
>
> SI <- rep(0.02,N)
> GImobility = N/60/8
>
> for (i in 1:N){
>   SI[i] <- (Si-Sf)/(1+exp(p0*((i/GImobility)-t0)))+Sf
> }
>
> Weibulla <- 2 # weibull function for time factor
> Weibullb <-0.9 # weibull fuction for shape factor
> Tlag <- 0 # lag time for disintegration
> criticT <- (15)^(1/Weibullb)*Weibulla # critical T for disintegration
>
> z <- 50
>
> ka <- 0.04
>
> para <- c(SGi <- 0.8, # inital solubility in Stomach unit is mg/ml
>           SGf <- 0.8,# final solubility in Stomach unit is mg/ml
>           TG <- 20, #time for half transformation in Stomach
>           PG <- 0.2) # rate of half transformation in Stomach
>
>  advModel <- function(t, state, parms) {
>
>      with (as.list(parms), {
>          Stomach1 <- state[1]  # undisintegrated
>          Stomach2 <- state[2]   # granule
>          Stomach3 <- state[3]   # Precipitate
>          Stomach4 <- state[4]   # dissolved
>          Intestine1 <- state[5:(N+4)]           # undisintegated
>          Intestine2 <- state[(N+5):(2*N+4)]     # granule
>          Intestine3 <- state[(2*N+5):(3*N+4)]   # precipitate
>          Intestine4 <- state[(3*N+5):(4*N+4)]   # dissolved
>          T1 <- state[4*N+5]   # total undisintegrated mass
>          T2 <- state[4*N+6]   # total granule
>          T3 <- state[4*N+7]   # total precipitate
>          T4 <- state[4*N+8]   # total in solution
>
>          absorption <- state[4*N+9]
>          plasma <- state[4*N+10]
>
>          SG <- SGi
>
>          ddisinte <- 0.03
>          # wettability factor
>
>          # define a function for dissolution kinetics
>          P_diss <- function(X,Y,V,S){
>            # X the solid remaining at time t, X0, the inital dose, r, the
> particle size
>            if (X <= 0) {
>              result <- 0  # avoid numerical error
>            } else if (S < Y/V){
>              result <- 0
>            } else {
>              result <- min(z*(X^(1/3))*(S-Y/V),X)
>            }
>            return(result)
>          }
>
>          # define the precipitation rate function
>
>          P_precip <- function(X,Y,V,S){
>            # X the solid remaining at time t, assuming the precipitate has
> particle size of 1 um
>            if (X <= 0) {
>              result <- 0  # avoid numerical error
>            } else if (Y/V > S){
>
>            result <- V * ( Y/V-S)
>            } else {
>              result <- -min(z*(X^(1/3))*(S-Y/V),X) # +  (V*( Y/V-S) +
> min(z*(X^(1/3))*(S-Y/V),X))/(1+exp(1000*(S-Y/V)))
>            }
>            return(result)
>          }
>
>
>          emptying <-
> 0.005+0.05*(tanh(10*(t%%240-210))-tanh(10*(t%%240-220)))
>          emptying1 <- 0.05*(tanh(10*(t%%240-210))-tanh(10*(t%%240-220)))
>
>
>          # dissolution in Stomach
>          disso1 <- P_diss(Stomach2,Stomach4,V1,SG)
>          precip1 <- P_precip(Stomach3,Stomach4,V1,SG)
>
>          dStomach1 <- -ddisinte*Stomach1 - emptying1*Stomach1
>          dStomach2 <- ddisinte*Stomach1 - disso1 - emptying*Stomach2
>          dStomach3 <- precip1 -emptying*Stomach3
>          dStomach4 <- disso1 - precip1 - emptying*Stomach4
>
>          disso2 <- rep(0,N)
>          precip2 <- rep(0,N)
>
>          for (i in 1:N){
>            disso2[i] <- P_diss(Intestine2[i],Intestine4[i],V2[i],SI[i])
>            precip2[i] <- P_precip(Intestine3[i],Intestine4[i],V2[i],SI[i])
>          }
>
>          Tran1 <- advection.1D(C = Intestine1, flux.up =
> emptying1*Stomach1, v = GImobility, dx=1,adv.method="super")
>          Tran2 <- advection.1D(C = Intestine2, flux.up =
> emptying*Stomach2, v = GImobility, dx=1,adv.method="super")
>          Tran3 <- advection.1D(C = Intestine3, flux.up =
> emptying*Stomach3, v = GImobility, dx=1,adv.method="super")
>          Tran4 <- advection.1D(C = Intestine4, flux.up =
> emptying*Stomach4, v = GImobility, dx=1,adv.method="super")
>
>
>         dIntestine1 <- Tran1$dC - ddisinte*Intestine1
>         dIntestine2 <- Tran2$dC - disso2 + ddisinte*Intestine1
>         dIntestine3 <- Tran3$dC + precip2
>         dIntestine4 <- Tran4$dC + disso2 - precip2 - ka*Intestine4
>
>          dT1 <- sum(dIntestine1)
>          dT2 <- sum(dIntestine2)
>          dT3 <- sum(dIntestine3)
>          dT4 <- sum(dIntestine4)
>
>
>         dabsorption <- sum(ka*Intestine4)
>         dplasma <- dabsorption -0.002*plasma
>
>
>          return (list(c(dStomach1,dStomach2,dStomach3,dStomach4,
>                         dIntestine1,dIntestine2,dIntestine3,dIntestine4,
>                         dT1,dT2,dT3,dT4,
>                         dabsorption,dplasma)))
>        })
> }
>
> times <- seq(1,1000,by =1)
>
> state <- c(Dose,rep(0,3),rep(0,4*N),rep(0,4),0,0)
>
> systime <- system.time(out <- ode.1D(func = advModel, y=state,times =times,
>                                      parms = para,nspec=1))
>
>

	[[alternative HTML version deleted]]

