[R-sig-dyn-mod] question on ode.1D
Karline Soetaert
Karline.Soetaert at nioz.nl
Tue Jul 30 13:04:37 CEST 2013
Dabeng,
Your complex example does not work, so I cannot really test it.
However, ode.1D only works efficiently if (1) all components are defined in an equal number of compartments (boxes), and (2) if for one component there is transition only between adjacent compartments, e.g. between boxes i-1, and i. (3) Also, interactions between different components can only occur if in the same box.
This is so because ode.1D was meant to efficiently solve 1-D reactive-transport equations.
If these conditions are not valid for your model, then ode.1D will be very inefficient.
If your model is really sparse, you could try to use lsodes.
Hope it helps
Karline
-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Dabing Chen
Sent: donderdag 25 juli 2013 21:36
To: R simulation
Subject: [R-sig-dyn-mod] question on ode.1D
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]]
_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
More information about the R-sig-dynamic-models
mailing list