R- ode 函数（deSolve 包）：将参数值更改为时间函数

我正在尝试使用该函数求解一阶微分方程ode来自deSolve包裹。问题如下：药物在某些时间（输注时间）以恒定的输注速率给药，并以一阶速率消除。因此，该过程可以描述为：

if(t %in% Infusion_times){Infusion <- Infusion_rate} else{Infusion <- 0}
    dC <- -Ke*C + Infusion

where t是时候了，Infusion_times是包含给药时间的载体，C是药物的量，Ke是它的消除常数，Infusion是一个变量，有输液时取输液速率的值，否则取值0。 因此，假设我们想要注射 3 剂，从时间 0、24 和 40 开始，每次输注持续两个小时，假设我们想要deSolve每 0.02 个时间单位计算一次答案。 我们想要deSolve例如，以 0.02 次单位的步长求解 0 到 48 次之间的微分方程。所以我们的向量指定了时间ode函数将是：

times <- seq(from = 0, to = 48, by = 0.02)

输注时间由下式给出：

Infusion_times <- c(seq(from = 0, to = 2, by = 0.02), seq(from = 24, to = 26, by = 0.02), 
                    seq(from = 40, to = 42, by = 0.02))

起初，我担心问题可能出在浮点数的比较上。为了防止这种情况，我将两个向量四舍五入到小数点后两位：

times <- round(times, 2)
Infusion_times <- round(times, 2)

所以现在，希望所有Infusion_times都包含在times vector:

(sum(Infusion_times %in% times)/length(Infusion_times))*100
[1] 100

如您所见，中的所有值Infusion_times（100%）包含在向量中times，因此变量Infusion应该取值Infusion_rate在指定的时间。然而，当我们解方程时，它不起作用。 让我们证明一下，但首先，让我们指定该函数所需的其他值ode功能：

parameters <- c(Ke = 0.5)
amounts <- c(C = 0) #Initial value for drug is 0
Inf_rate <- 5

现在，让我们根据需要编写一个函数来说明变化率：

OneComp <- function(t, amounts, parameters){
  with(as.list(c(amounts, parameters)),{
      if(t %in% Infusion_times){Infuse =Inf_rate} else{Infuse = 0}
      dC <- -Ke*C + Infuse
  list(c(dC))})
}

对于那些不熟悉的人deSolve包，参数t函数的值将说明方程应该积分的时间，amounts将指定 C 的初始值和parameters将给出参数的值（在本例中，只是Ke）。 现在，让我们解方程：

out <- ode(func = OneComp, y = amounts, parms = parameters, times = times)

让我们绘制结果：

library(ggplot2)

ggplot(as.data.frame(out)) + geom_line(aes(x = times, y = C))

这与我们得到的结果完全相同Infusion始终等于 0。 然而，请注意，如果我们只模仿单剂量，并且我们要尝试类似的方法，它会起作用：

OneComp <- function(t, amounts, parameters){
      with(as.list(c(amounts, parameters)),{
          if(t < 2){Infuse =Inf_rate} else{Infuse = 0}
          dC <- -Ke*C + Infuse
      list(c(dC))})
    }
out <- ode(func = OneComp, y = amounts, parms = parameters, times = times)
ggplot(as.data.frame(out)) + geom_line(aes(x = times, y = C))

在这里，我们设置了变量Infuse取值Inf_rate仅当时间少于2小时时，它有效！因此，我对这些行为完全感到困惑。这不是更改变量值的问题，也不是浮点数之间比较的问题......知道这些可能是什么吗？谢谢

大多数求解器deSolve使用自动内部时间步长，根据系统的粗糙度或平滑度进行自我调整。指某东西的用途if陈述或if()模型函数中的函数不是一个好主意，原因有两个：(i) 时间步长可能无法准确命中，(2) 模型函数（即导数）理想情况下应为连续且可微并避免逐步行为，即使求解器在这种情况下非常强大。

The deSolve包提供了解决您的问题的两种方法：“强制函数”和“事件”。两者都有其优点和缺点，但如果“事件”（例如注射）的时间与积分时间步相比非常短，“事件”特别有用。

有关此内容的更多信息，请参阅deSolve帮助页面?forcings and ?events, in deSolve：强制函数和事件来自 userR!2017 会议，以及slides来自用户R!2014。

请检查以下内容是否适合您：

library("deSolve")

OneComp <- function(t, y, parms){
  with(as.list(c(y, parms)),{
    dC <- -Ke * C
    list(c(dC))
  })
}

eventfunc <- function(t, y, parms) {
  with(as.list(c(y, parms)),{
    C + Inf_rate
  })
}

parms <- c(Ke = 0.5, Inf_rate = 5)

y0 <- c(C = 0)            # Initial value for drug is 0

Infusion_times <- c(seq(from =  0, to =  2, by = 0.02), 
                    seq(from = 24, to = 26, by = 0.02), 
                    seq(from = 40, to = 42, by = 0.02))

# time step can be made bigger to speedup simulation
times <- round(seq(from = 0, to = 48, by = .1) , 2)

# check that all events are in 'times', but no duplicates
# this check is also done by the solver and may print a warning
# to ensure that the user is very careful with this
unique_times <- cleanEventTimes(times, Infusion_times)
times        <- sort(c(unique_times, Infusion_times))

out <- ode(func = OneComp, y = y0, parms = parms, times = times, 
           events = list(func = eventfunc, time = Infusion_times))

plot(out)
rug(Infusion_times)

The two cleanEventTimes线是确保模拟命中所有事件时间的一种可能方法。它通常由求解器自动完成，并且可能会发出警告以提醒用户对此非常小心。

我使用了“基本”情节rug来指示注射时间。

我对这些条款有些疑问Infusion_times and Inf_rate。在基于事件的方法中，“量”在离散时间点被添加到状态变量 C，而“速率”则表示在一个时间点内连续添加。时间间隔。这通常称为强制函数。

强制函数会更简单并且在数值上更好。