deSolve：在 (t-tau) 处获取值

Question

我正在尝试将已发布的模型运行转化为 deSolve，我想知道如何在 t-tau 处获得第三个方程 dH/dt <- phi*B(t-tau)-(amin+amax*(b^2/(b^2+f^2))-sigma*(F/(F+H)))*H 处 B 的值] 我不能得到它运行所以现在我已经删除了 (t-tau) 部分并且可以运行 R 中的模型没有它，但这是我想要的机制的关键部分学习。

library(deSolve)
# Full model
# rigidode <- function(t, y, parms) {
#   df <- c*N*F-ya*(F+H)-yb*B
#   dB <- L*(f^2/(f^2+b^2))*(H/(H+v))-phi*B
#   dH <- phi*B(t-tau)-(amin+amax*(b^2/(b^2+f^2))-sigma*(F/(F+H)))*H
#   dF <- T(a)*(amin+amax*(b^2/(b^2+f^2))-sigma*(F/(F+H)))*H-mr*M(a)*F
#   list(c(dy1, dy2, dy3))
# }

#Now here is the numerical simulation
rigidode <- function(t, y, parms) {
  dy1 <- c*N*y[4]-ya*(y[4]+y[3])-yb*y[2]
  dy2 <- L*(y[1]^2/(y[1]^2+b^2))*(y[3]/(y[3]+v))-phi*y[2]
  dy3 <- phi*y[2]-(amin+amax*(b^2/(b^2+y[1]^2))-sigma*(y[4]/(y[4]+y[3])))*y[3] #y[2] should be "y[2](t-tau) as B(t-tau)"
  dy4 <- Ta*(amin+amax*(b^2/(b^2+y[1]^2))-sigma*(y[4]/(y[4]+y[3])))*y[3]-m*y[4]
  list(c(dy1, dy2, dy3, dy4))
}

yini <- c(y1=1000, y2=0, y3=16000, y4=8000)
parms <- c(c<-0.033,N<-9,ya<-0.007,yb<-0.0018,L<-2000,b<-500,v<-5000,phi<-1/21,amin<-0.25,amax<-0.25,sigma<-1.3,m<-0.05,Ta<-0.75,tau<-12)
times <- seq(from = 0, to = 600, by = 1)
out <- ode (times = times, y = yini, func = rigidode, parms = parms)
head (out, n = 3)
plot(out)

Answer 1

老实说，我不确定我的回答是否正确，但至少给你有用的信息。
deSolve 包有一个函数 lagvalue() 来访问状态变量的过去（滞后）值。
?lagvalue 表示它将与函数 dede 一起使用（但可能与 ode() 一起使用）。

rigidode <- function(t, y, parms) {

  if (t < tau)
    lagged_y2 <- 0
  else
    lagged_y2 <- lagvalue(t - tau)[2]

  if(t == 100) print(lagged_y2)  # for check

  dy1 <- c*N*y[4]-ya*(y[4]+y[3])-yb*y[2]
  dy2 <- L*(y[1]^2/(y[1]^2+b^2))*(y[3]/(y[3]+v))-phi*y[2]
  dy3 <- phi*lagged_y2-(amin+amax*(b^2/(b^2+y[1]^2))-sigma*(y[4]/(y[4]+y[3])))*y[3]  # modified line
  dy4 <- Ta*(amin+amax*(b^2/(b^2+y[1]^2))-sigma*(y[4]/(y[4]+y[3])))*y[3]-m*y[4]
  list(c(dy1, dy2, dy3, dy4))
}


# parameters are the same as yours.
yini <- c(y1=1000, y2=0, y3=16000, y4=8000)
parms <- c(c<-0.033,N<-9,ya<-0.007,yb<-0.0018,L<-2000,b<-500,v<-5000,phi<-1/21,
           amin<-0.25,amax<-0.25,sigma<-1.3,m<-0.05,Ta<-0.75,tau<-12)
times <- seq(from = 0, to = 600, by = 1)

out <- dede(times = times, y = yini, func = rigidode, parms = parms) # use dede

# printed_value
[1] 37022.74
[1] 37022.74

out[100 - 12 + 1,]
#  time        y1        y2        y3        y4 
# 88.00 157425.09  37022.74  59146.39  12862.78

deSolve：在 (t-tau) 处获取值

deSolve: get value at (t-tau)

modeling

r