Добавление оператора if then к начальному значению условия в системе ODE; deSolve

Question

Я пытаюсь добавить оператор if then для определения начального значения одной из моих переменных состояния и использую deSolve. По сути, я хочу ввести 3-е ОДУ (в данном случае 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 начальное значение 3-й переменной состояния будет 0, и что в момент времени = 50 или выше начальное значение 3-я переменная состояния будет 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

Для меня не совсем ясно, что подразумевается под утверждением, что «начальное значение 3-ей переменной состояния» должно быть 1 для t>= 50. Начальное значение определяет начало переменной состояния, которая затем изменяется дифференциальные уравнения. Ниже я покажу следующие подходы:

1. Переменная состояния Pj инициализируется заданным значением при t = 50 . Это может быть обработано событием .
2. Переменная состояния Pj получает дополнительный внешний вход при t> = 50 . Это можно сделать с помощью внешнего сигнала, также называемого функцией принуждения .

В первом примере показан механизм событий, реализованный в виде фрейма данных. . Он также может быть реализован в более гибкой форме с функцией события.

Здесь я увеличил значение «начального» состояния при t = 50 до 100, чтобы эффект был более выраженным. Округление вектора времени сделано, чтобы избежать предупреждения (спросите, хотите ли вы знать почему).

      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. Следует ли использовать функция или нет - дело вкуса. У обоих есть свои плюсы и минусы.

Но, использовать ли функцию события или принудительную функцию, имеет большое значение.

      
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))