r-在ODE系统中，向条件初始值添加if-then语句;deSolve

Question

我正在尝试添加一个if-then语句来调节我的一个状态变量的初始值，并且正在使用deSolve。从本质上讲，我想在模拟开始后引入第三个ODE(在本例中是种群中的第三个物种(。

以下是没有条件的代码：

Antia_3sp_Model <- function(t,y,p1){
# Parms
ri <- p1[1]; rj <- p1[2]; k <- p1[3]; p <- p1[4]; o <- p1[5] 
# State vars
Pi <- y[1]; Pj <- y[2]; I <- y[3]
# ODEs
dPi = ri*Pi - k*Pi*I
dPj = rj*Pj - k*Pj*I
dI = p*I*(Pi/(Pi + o) + Pj/(Pj + o))
list(c(dPi,dPj,dI))
}
# Parm vals
ri <- 0.3; rj <- 0.2; k <- 0.001; p <- 1; o <- 1000 # Note that r can range btw 0.1 and 10 in this model
parms <- c(ri,rj,k,p,o)
# Inits
Pi0 <- 1; Pj0 <- 1; I0 <- 1
N0 <- c(Pi0,Pj0,I0)
# Time pt. sol'ns
TT <- seq(0.1,200,0.1)
# Sim
results <- lsoda(N0,TT,Antia_3sp_Model,parms,verbose = TRUE)

到目前为止，在尝试添加if-then语句后，我得到了这样的结果，即在时间=50之前，第三个状态变量的初始值将为0，而在时间=5或以上，第三状态变量的起始值将为1。

Antia_3sp_Model <- function(t,y,p1){
# Parms
ri <- p1[1]; rj <- p1[2]; k <- p1[3]; p <- p1[4]; o <- p1[5] 
# State vars
Pi <- y[1]; Pj <- y[2]; I <- y[3]

if (t[i] < t[50]){
Pj0 = 0
}
else if (t[i] >= t[50]){
Pj0 = 1
}

# ODEs
dPi = ri*Pi - k*Pi*I
dPj = rj*Pj - k*Pj*I 
dI = p*I*(Pi/(Pi + o) + Pj/(Pj + o))
list(c(dPi,dPj,dI))
}
# Parm vals
ri <- 0.3; rj <- 0.2; k <- 0.001; p <- 1; o <- 1000 # Note that r can range btw 0.1 and 10 in this model
parms <- c(ri,rj,k,p,o)
# Inits
Pi0 <- 1; Pj0 <- 1; I0 <- 1
N0 <- c(Pi0,Pj0,I0)
# Time pt. sol'ns
TT <- seq(0.1,200,0.1)
# Sim
results <- lsoda(N0,TT,Antia_3sp_Model,parms,verbose = TRUE)

有什么建议吗？

请让我知道我是否应该添加任何其他信息，非常感谢您的阅读！：(

Answer 1

对我来说，"；第三状态变量"0"的初始值；对于t>=应为150.初始值定义了一个状态变量的开始，然后通过微分方程进行演化。在下文中，我展示了以下方法：

1. 状态变量Pj被初始化为t=50时的给定值。这可以通过事件来处理
2. 状态变量Pj在t>=50。这可以通过外部信号处理，也称为强制函数

第一个示例显示了作为数据帧eventdat实现的事件机制。它也可以通过事件函数以更灵活的形式来实现。

在这里我增加了"；首字母"；t＝50到100时的状态值，以使效果更加明显。对时间矢量TT进行舍入是为了避免出现警告(如果您想知道原因，请询问(。

library("deSolve")
Antia_3sp_Model <- function(t, y, p1){
# Parms
ri <- p1[1]; rj <- p1[2]; k <- p1[3]; p <- p1[4]; o <- p1[5]
# State vars
Pi <- y[1]; Pj <- y[2]; I <- y[3]
# ODEs
dPi <- ri*Pi - k*Pi*I
dPj <- rj*Pj - k*Pj*I
dI <- p*I*(Pi/(Pi + o) + Pj/(Pj + o))
list(c(dPi, dPj, dI))
}
parms <- c(ri = 0.3, rj = 0.2, k = 0.001, p = 1, o = 1000)
N0 <- c(Pi = 1, Pj = 1, I = 1)
TT <- round(seq(0.1, 200, 0.1), 1)
## An "initial value" is the value at the beginning. We call the value during
## simulation the "state". If it is meant that the state should be changed at
## a certain point of time, it can be done with an event
# tp: initial value at t=50 set to 100 to improve visibility of effect (was 1)
eventdat <- data.frame(var = "Pj", time = 50, value = 100, method = "rep")
results <- lsoda(N0, TT, Antia_3sp_Model, parms, events=list(data=eventdat), verbose = TRUE)
plot(results, mfcol=c(1, 3))

强制函数可用于实现与时间相关的参数或向状态连续添加常数值。还要注意ODE模型的紧凑样式。是否使用with函数是一个品味问题。两者各有利弊。

但是，使用事件还是强制函数会产生很大的不同。

Antia_3sp_Model <- function(t, y, p, import){
with(as.list(c(y, p)), {
dPi <- ri*Pi - k*Pi*I
dPj <- rj*Pj - k*Pj*I + import(t)
dI <- p*I*(Pi/(Pi + o) + Pj/(Pj + o))
list(c(dPi, dPj, dI))
})
}
signal <- approxfun(x=c(0, 50, max(TT)), y=c(0, 1, 1), method="constant", rule=2)
results <- lsoda(N0, TT, Antia_3sp_Model, parms, import=signal, verbose = TRUE)
plot(results, mfcol=c(1, 3))