deSolve:获得(t-τ)的价值

Question

我试图在deSolve中运行一个已发布的模型，我想知道如何才能在t-tau上得到第三个方程dH/dt <- phi*B(t-tau)-(amin+amax*(b^2/(b^2+f^2))-sigma*(F/(F+H)))*H的值--我无法让它运行--所以现在我已经删除了(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