时间戳

首页 > 解决方案 > R- ode 函数（deSolve 包）：将参数的值更改为时间的函数

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

问题描述

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

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

其中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向量中：

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

如您所见，Infusion_times(100%) 中的所有值都包含在 vectortimes中，因此变量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))

如果总是等于 0，这与我们得到的结果完全相同。Infusion但是，请注意，如果我们只模拟一个剂量，并且我们要尝试类似的方法，它会起作用：

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))

在这里，我们让变量只在时间小于 2 小时的时候Infuse取值，就可以了！ 因此，我对这些行为完全感到困惑。这不是改变变量值的问题，也不是浮点数之间的比较问题......知道这些可能是什么吗？谢谢Inf_rate

标签： rodeequation-solving

解决方案

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

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

有关这方面的更多信息，请参见deSolve帮助页面?forcings和deSolve： userR!2017 会议中的强制函数和事件，以及 userR!2014 中的幻灯片 。?events

请检查以下是否适合您：

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)

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

我使用“基础”图并rug指示注射时间。

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

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

推荐阅读