「R/SでODE」の編集履歴（バックアップ）一覧に戻る

R/SでODE - (2007/03/17 (土) 07:47:21) のソース

*R で ODE

 ex1 <-
 function(t, y, p) {
   dy1 <- -p["ka"] * y[1] # amount
   dy2 <-  p["ka"] * y[1] / p["V"] - p["CL"] * y[2] # concentration
 
   list(c(dy1, dy2)) # list で返す
 }
 
 ka <- 0.7; CL <- 0.1; V <- 1
 parms <- c(ka=ka, CL=CL, V=V)
 
 TIME <- seq(0, 24, by=1)
 Dose <- 1000
 
 require(odesolve)
 
 my.atol <- c(1e-6, 1e-10)
 out <-
 lsoda(c(Dose, 0), TIME, ex1, parms, rtol=1e-5, atol=my.atol)
 
 out

*S-PLUS で ODE

 ex1 <-
 function(x, y, CL, V, ka) {
   dy1 <- -ka * y[1] # amount
   dy2 <-  ka * y[1] / V - CL * y[2] # concentration
   c(dy1, dy2)
 }
 
 TIME <- seq(1, 24, by=1)
 CONC <- double(length(TIME))
 
 CL <- 0.1; V <- 1; ka <- 0.7
 Dose <- 1000
 out <- ivp.ab(TIME[1], c(0, c(Dose, 0)), ex1, aux=list(CL=CL, V=V, ka=ka))
 CONC[1] <- out$values[3]
 
 for (i in seq(2, length(TIME))) {
   out <- ivp.ab(TIME[i], restart=out)
   CONC[i] <- out$values[3]
 }
 
 CONC