循环参数以使用 deSolve 获得平衡

Question

我直觉地与循环作斗争。我有一个简单的消费者资源模型，我想遍历资源增长率g 的值以获取最终状态值，然后将平衡绘制为参数值的函数。这是我目前所拥有的：

param.values = seq(from = 1, to = 10, by = 1)
variable = rep(0,length(param.values))
for (i in 1:length(param.values)){ 
  state <- c(r = 1, n = 1)
  parameters = c(g = variable[i],# resource growth rate
                 d = 0.5, # n mortality rate
                 k = 5, # r carrying capacity
                 c = 1, # consumption rate of n on r
                 e = 1, # conversion efficiency for n on r
                 h = 1 # handling time n on r
  )
  function1 <- function(times, state, parameters) {
    with(as.list(c(state, parameters)),{
      # rate of change
      dr = variable[i]*r*(1 - (r/k)) - (c*n*r/(1+(h*c*r))) 
      dn = (e*c*n*r/(1+(h*c*r)))- n*d
      
      # return the rate of change
      list(c(dr, dn))
    }) 
  }
  times <- seq(0, 100, by = 1)
  
  out <- ode(y = state, times = times, func = function1, parms = parameters)
  
  sol <- out[101, 2:3] # trying to get last equilibrium value to plot against param values...
  print(sol[i])
}

plot(sol[,1] ~ param.values)
plot(sol[,2] ~ param.values)

我想我一直在思考直到 ode 函数 - 在此之后我应该在哪里索引 i？我希望这是有道理的。

Answer 1

你的方法有几个问题，所以我试图重新组织它，让它贯穿始终。但是，由于您的模型显示出一个稳定的循环，它并没有达到平衡。

这里有一些提示

• 循环应该只包含在模拟过程中发生变化的东西。固定代码段应该在循环之前。这更易于维护且速度更快。
• 首先，在没有循环的情况下运行模型，看看它是否有效。
• 然后定义一个数据结构（矩阵或数据框）来存储结果。

这里有一种实现方法：

library("deSolve")

## define as much as possible outside the loop
function1 <- function(times, state, parameters) {
  with(as.list(c(state, parameters)),{
    # rate of change
    dr = g*r*(1 - (r/k)) - (c*n*r/(1+(h*c*r)))
    dn = (e*c*n*r/(1+(h*c*r)))- n*d

    # return the rate of change
    list(c(dr, dn))
  })
}

state <- c(r = 1, n = 1)

parameters = c(g = 1,   # resource growth rate
               d = 0.5, # n mortality rate
               k = 5,   # r carrying capacity
               c = 1,   # consumption rate of n on r
               e = 1,   # conversion efficiency for n on r
               h = 1    # handling time n on r
)
times <- seq(0, 100, by = 1)

## first test single run of model
out <- ode(y = state, times = times, func = function1, parms = parameters)
plot(out)
## It runs and we see a cycling model. I suspect it has no equilibrium!

param.values = seq(from = 1, to = 10, by = 1)

## define a matrix where results can be stored
sol <- matrix(0, nrow=length(param.values), ncol=2)

for (i in 1:length(param.values)){

  ## replace single parameter g with new value
  parameters["g"] <- param.values[i]
  out <- ode(y = state, times = times, func = function1, parms = parameters)

  ## store result of last value in row of matrix.
  ## Note that it may not be an equilibrium
  sol[i, ] <- out[101, 2:3] # trying to get last equilibrium value to plot against param values...
  print(sol[i, ])
}

plot(sol[,1] ~ param.values, type="l")
plot(sol[,2] ~ param.values, type="l")
## We see that the model has no equilibrium.

还有其他几种方法，如前所述，该模型没有平衡。这里是另一个model example，所谓的恒化器具有平衡。