如何解决在 ODE 系统中添加 if then 语句来条件初始值;解开
我正在尝试添加一个 if then 语句来调节我的状态变量之一的初始值,并且正在使用 deSolve。本质上,我想在模拟开始后引入第 3 个 ODE(在本例中为第 3 个物种)。
没有条件的代码如下所示:
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 语句之后,说在 time = 50 之前,第三个状态变量的初始值将为 0,并且在 time = 50 或以上,初始值第三个状态变量的值为 1。
Antia_3sp_Model <- function(t,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,verbose = TRUE)
有什么建议吗?
如果我应该添加任何其他信息,请告诉我,非常感谢您的阅读! :)
解决方法
对我来说,“第三个状态变量的初始值”在 t >= 50 时应为 1 的含义并不完全清楚。初始值定义了状态变量的开始,然后由微分方程演化。下面,我将展示以下方法:
- 状态变量 Pj 在 t = 50 时初始化为给定值。这可以通过事件处理。
- 状态变量 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,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,dI))
})
}
signal <- approxfun(x=c(0,50,max(TT)),y=c(0,1,1),method="constant",rule=2)
results <- lsoda(N0,import=signal,3))
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。