R S用事件解决物质平衡问题

判断质量平衡需要close the system.这可以通过添加源和汇的显式状态变量来实现，其中source代表药物的存量，sink代表通过分解或排泄离开血液的累积量.源变量可以设置为零以判断质量平衡，或者设置为任何正值以表示库存中可用的药物数量.

然后在模型中加入相应的过程，即分解/排泄作为状态方程，以及在事件中注入过程中从原料中移除.然后，该模型的内容如下:

library(deSolve)
b <- 0.6
yini <- c(source=0, blood = 0, sink=0)

pharmaco2 <- function(t, yini, p) {
  blood   <-  yini[2]
  dsource <-  0
  dblood  <- -b * blood
  dsink   <- +b * blood
  list(c(dsource, dblood, dsink))
}

injectevents <- data.frame(var = rep(c("blood", "source"), 21),
                           time = rep(0:20, each = 2),
                           value =  rep(c(40, -40), 21),
                           method = "add")

times <- seq(from = 0, to = 10, by = 1/24)

out2 <- ode(func = pharmaco2,
            times = times,
            y = yini,
            parms = NULL,
            method = "lsoda",
            events = list(data = injectevents))

plot(out2, lwd = 2, xlab="days", mfrow=c(1,3))

summary(rowSums(out2[,-1])) # sum of all columns, except time

质量平衡判断可以通过添加所有列的值来完成，时间除外.如果初始值source被设置为零，则余额也应该为零，直到数字四舍五入误差:

      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-1.990e-13 -4.974e-14 -2.309e-14 -2.541e-14  2.132e-14  9.948e-14 

另请注意，我使用默认的lsoda作为求解器.求解器adams或bdf也可以工作，impAdams也可以，但我不建议在这种情况下使用后者.