对每个网格单元应用一组常微分方程 [英] Applying a set of ordinary differential equations to each grid cell

问题描述

我正在开发一种基于主体的模型，以模拟传染病在由栖息地多边形(或相连细胞团)组成的异质景观中的传播.为了简化模型，我考虑了包含每个像元的多边形ID的栖息地网格(或栅格).此外，我还具有与每个多边形ID相关的流行病学参数.在每个时间步中，参数值在面中都会变化.因此，数据帧 landscape (请参见下文)在每个时间步都更新.这是t = 0时的示例:

I am developing an agent-based model to simulate the spread of infectious diseases in heterogeneous landscapes composed of habitat polygons (or clumps of connected cells). To simplify the model, I consider a habitat grid (or raster) containing the polygon ID of each cell. In addition, I have epidemiological parameters associated with each polygon ID. At each time step, the parameter values change in the polygon. Thus, the data frame landscape (see below) is updated at each time step. Here is an example at t = 0:

landscape <- data.frame(polygon_ID = seq(1, 10, by = 1), 
                        beta = sample(c(100, 200, 400, 600), 10, replace = TRUE), 
                        gamma = sample(c(25, 26, 27, 28), 10, replace = TRUE))

为了研究疾病动态，我还基于常微分方程(ODE)系统开发了一个隔间模型.这是一个代表ODE系统的示例:

To study the disease dynamics, I also am developing a compartmental model based on a system of ordinary differential equations (ODEs). Here is an example to represent the system of ODEs:

solve_sir_model <- function (times, parameters) { 

  sir_model <- function (times, states, parameters) { 

    with(as.list(c(states, parameters)), { 

      dSdt <- -beta*S*I 
      dIdt <- beta*S*I-gamma*I 
      dRdt <- gamma*I 
      dNdt <- dSdt + dIdt + dRdt 

      return(list(c(dSdt, dIdt, dRdt, dNdt))) 

    }) 
  } 

  states <- c(S = 99, I = 1, R = 0, N = 100) 
  return(ode(y = states, times = times, func = sir_model, parms = parameters)) 
} 

require(deSolve) 
output <- as.data.frame(solve_sir_model(times = seq(0, 5, by = 1), parameters = c(beta = 400, gamma = 28)))

是否可以在每个时间步将ODE系统应用于数据框 landscape 中的每个栖息地多边形(因此每一行)?我正在使用lsoda作为ODE求解器.我是否需要在每个时间步骤上使用另一个求解器来应用ODE?

At each time step, is it possible to apply the system of ODEs to each habitat polygon (thus each row) in the data frame landscape? I am using lsoda as an ODE solver. Do I need to use another solver to apply the ODEs at each time step?

编辑

在我的情况下，函数 ode 中的方法 iteration 似乎很有用:

It seems that the method iteration in the function ode can be useful in my case:

方法迭代"的特殊之处在于，此处的函数func应该返回状态变量的新值，而不是比率改变.可以用于基于个人的模型，以区别等式，或在其中进行积分的情况下功能).

Method "iteration" is special in that here the function func should return the new value of the state variables rather than the rate of change. This can be used for individual based models, for difference equations, or in those cases where the integration is performed within func).

我已经测试了该方法，但是我不明白为什么它不能在单个时间步中起作用:

I have tested the method but I don't understand why it doesn't work with a single time step:

solve_sir_model <- function (times, parameters) { 

  sir_model <- function (times, states, parameters) { 

    with(as.list(c(states, parameters)), { 

      dSdt <- -beta*S*I 
      dIdt <- beta*S*I-gamma*I 
      dRdt <- gamma*I 
      dNdt <- dSdt + dIdt + dRdt 

      return(list(c(dSdt, dIdt, dRdt, dNdt))) 

    }) 
  } 

  states <- c(S = 99, I = 1, R = 0, N = 100) 
  return(ode(y = states, times = times, func = sir_model, parms = parameters, method = "iteration")) 
} 

require(deSolve) 
output <- as.data.frame(solve_sir_model(times = 1, parameters = c(beta = 400, gamma = 28)))

 Error in iteration(y, times, func, parms, ...) : 
   times should be equally spaced In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

推荐答案

类似的东西.

# Library *not* require
library(deSolve)

# Parameters: number of polygons, beta, & gamma
n.polys <- 10 
beta <- sample(c(100, 200, 400, 600), n.polys, replace = TRUE)
gamma <- sample(c(25, 26, 27, 28), n.polys, replace = TRUE)

# Model defintion
sir_model <- function (t, Y, pars) { 
  # Break up state variable into parts
  S <- Y[1:n.polys]
  I <- Y[(n.polys+1):(2 * n.polys)]
  R <- Y[(2 * n.polys+1):(3 * n.polys)]

  # Calculate rate of change
  dSdt <- -beta * S * I 
  dIdt <- beta * S * I - gamma * I 
  dRdt <- gamma * I 

  # Return rates of change after concatenating them 
  return(list(c(dSdt, dIdt, dRdt))) 
} 

# Initial conditions
Y.ini <- c(S = rep(99, n.polys), I = rep(1, n.polys), R = rep(0, n.polys)) 

# Solve the model
ode(y = Y.ini, times = 0:5, func = sir_model, parms = NULL) 

这篇关于对每个网格单元应用一组常微分方程的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！