Трудность запуска модели ODE в R, включая параметр, который изменяется во времени (функция форсирования)

Question

Я пытаюсь подобрать достаточно простую модель ODE, используя deSolve, включая в модели параметр, который изменяется во времени (сила заражения; FOI). Хотя модель без этого параметра работает нормально, включение параметра, зависящего от времени, приводит к ошибке (см. Ниже).

Я относительно новичок в R и математическом моделировании и уже некоторое время пытаюсь решить эту проблему.

Я создал параметр FOI в качестве матрицы значений, а затем использовал функцию приближения для интерполяции (как я уже видел, это работает с форсирующими функциями, например, https://rdrr.io/rforge/deSolve/man/forcings.html).

Модель без этого зависящего от времени параметра работает без каких-либо ошибок, но при попытке включить его выдает ошибку:

Error in checkFunc(Func2, times, y, rho) : 
  The number of derivatives returned by func() (200) must equal the 
length of the initial conditions vector (2)

Я не могу понять, как исправить эту ошибку, поскольку у меня есть только 2 начальных условия, и кажется, что включение этого зависящего от времени параметра FOI генерирует гораздо больше производных.

Мне известно, что другие задавали аналогичный вопрос, но я не нашел этот вопрос в отношении принудительных функций.

Большое спасибо заранее за любые советы.

# Forcing function data

foi <- matrix(ncol=2,byrow=TRUE,data=c(
  0, 0.003, 2, 0.03, 3, 0.08, 4,0.1,  5, 0.12, 6, 0.15, 
  8, 0.16, 10, 0.14,12, 0.12,14,0.08,15, 0.06,16, 0.03,
  17, 0.01,18,0.003,19,0.003,20,0.003,30,0.003,40,0.003,
  50,0.003,60,0.003,65,0.01, 70,0.08,  72,0.095,74,0.10, 
  76,0.1, 78,0.08,  80,0.06))

age <- seq(0, 80, by = 1)

input <- approxfun(x = foi[,1], y = foi[,2], method = "constant", rule = 2)

# Function
ab <- function(time, state, pars) {
  with(as.list(c(state, pars)), {

import<-c(input(t))

diggP<- (import *iggN) - iggR*iggP   
diggN<- (-import*iggN) + iggR*iggP

return(list(c(diggP, diggN))) 
  })
}


# Initial values
yini  <- c(iggP=0, iggN=1) 

# Parameters
pars  <- c(iggR = 0, import)

# ODE solver
results<- ode(y=yini, times=age, func=foi_model, pars)

Я надеюсь создать модель, в которой в каждый момент времени (или в данном случае возраст) FOI варьируется в зависимости от значений, которые я ввел в матрицу FOI. Поэтому я хотел бы увидеть, как изменение FOI с возрастом влияет на результаты дифференциальных уравнений.

Answer 1

Ваша основная проблема заключалась в том, что вы передавали аргумент t в input, но эта переменная не существует в вашем коде.Время передается вашей модели в качестве аргумента time.(Кроме того, ваша модель называется ab, а не foi_model, как указано в вызове к ode, плюс pars не требуется import и ее следует передавать на ode.)

# Load library
library(deSolve)

# Create FOI matrix
foi <- matrix(ncol=2,byrow=TRUE,data=c(
  0, 0.003, 2, 0.03, 3, 0.08, 4,0.1,  5, 0.12, 6, 0.15, 
  8, 0.16, 10, 0.14,12, 0.12,14,0.08,15, 0.06,16, 0.03,
  17, 0.01,18,0.003,19,0.003,20,0.003,30,0.003,40,0.003,
  50,0.003,60,0.003,65,0.01, 70,0.08,  72,0.095,74,0.10, 
  76,0.1, 78,0.08,  80,0.06))

# Times for model solution
age <- seq(0, 80, by = 1)

# Linear interpolation function from FOI data
input <- approxfun(x = foi[,1], y = foi[,2], method = "constant", rule = 2)

# Model to be integrated
ab <- function(time, state, parms) {
  with(as.list(c(state, parms)), {

    ##### IMPORTANT #####
    import<-input(time) #<- 'time' was previously 't'
    #####################

    # Derivatives
    diggP<- (import *iggN) - iggR*iggP   
    diggN<- (-import*iggN) + iggR*iggP

    # Return results
    return(list(c(diggP, diggN))) 
  })
}

# Initial values
yini  <- c(iggP=0, iggN=1) 

# Parameters
pars  <- c(iggR = 0)

# Solve model
results<- ode(y=yini, times=age, func=ab, parms = pars)

# Plot results
plot(results)

^{Создано в 2019-03-27 с помощью представительного пакета (v0.2.1)}