如何在 R 中求解具有时间相关参数的 ODE 系统？

Question

我正在尝试通过 deSolve 求解这个 ODE 系统，dX/dt = -X*a + (Y-X)b + c 和 dY/dt = -Ya + (X-Y)*b 对于时间 [0,200]，a=0.30，b=0.2 但对于时间 [50,70] c 为 1，否则为 0。我一直在使用的代码是，

time <- seq(0, 200, by=1)
parameters <- c(a=0.33, b=0.2, c=1)
state <- c(X = 0, Y = 0)

    two_comp <- function(time, state, parameters){
      with(as.list(c(state, parameters)), {
        dX = -X*a + (Y-X)*b + c
        dY = -Y*a + (X-Y)*b
        return(list(c(dX, dY)))
      })
    }

out <- ode(y = state, times = time, func = two_comp, parms = parameters)
out.df = as.data.frame(out)

我遗漏了 c 参数的时变部分，因为我想不出一种方法来顺利地包含它并运行它。我尝试将它包含在函数定义中，但无济于事。

Answer 1

c 等于 1 的区间限制可以作为 parameters 传递。然后，在微分函数内部，使用它们创建一个逻辑值

time >= lower & time <= upper

由于 FALSE/TRUE 被编码为整数 0/1，每次此条件为假时，c 乘以零，技巧就完成了。

library(deSolve)

two_comp <- function(time, state, parameters){
  with(as.list(c(state, parameters)), {
    dX = -X*a + (Y-X)*b + c*(time >= lower & time <= upper)
    dY = -Y*a + (X-Y)*b
    return(list(c(dX, dY)))
  })
}

time <- seq(0, 200, by=1)
parameters <- c(a=0.33, b=0.2, c=1, lower = 50, upper = 70)
state <- c(X = 0, Y = 0)

out <- ode(
  y = state, 
  times = time, 
  func = two_comp, 
  parms = parameters
)

out.df <- as.data.frame(out)
head(out.df)

matplot(out.df$time, out.df[-1], type = "l", lty = "solid", ylim = c(0, 3))
legend("topright", legend = names(out.df)[-1], col = 1:2, lty = "solid")

Answer 2

标准方法是使用approxfun，即创建一个时间相关信号，我们也称之为forcing variable:

library("deSolve")
time <- seq(0, 200, by=1)
parameters <- c(a=0.33, b=0.2, c=1)
state <- c(X = 0, Y = 0)

two_comp <- function(time, state, parameters, signal){
  cc <- signal(time)
  with(as.list(c(state, parameters)), {
    dX <- -X * a + (Y - X) * b + cc
    dY <- -Y * a + (X - Y) * b
    return(list(c(dX, dY), c = cc))
  })
}

signal <- approxfun(x = c(0, 50, 70, 200), 
                    y = c(0, 1,  0,  0), 
                    method = "constant", rule = 2)

out <- ode(y = state, times = time, func = two_comp, 
           parms = parameters, signal = signal)
plot(out)

另请注意 deSolve 特定 plot 函数，并且时间相关变量 cc 用作附加输出变量。

可以找到更多相关信息：

在 ?forcings 帮助页面和
在 Github 上的 short tutorial 中。

如何在 R 中求解具有时间相关参数的 ODE 系统？

How to solve a system of ODE with time dependent parameters in R?

r

ode

differential-equations

equation-solving

desolve