deSolve: получить значение в (т-тау)

Question

Я пытаюсь запустить опубликованную модель в deSolve, и мне интересно, как я могу получить значение B при 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

Честно говоря, я не уверен, что мой ответ правильный, но, по крайней мере, даю вам полезную информацию. * В пакете 1001 * 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