用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.

还有其他几种方法，就像前面说的，这个模型没有均衡。在这里，另一个模型实例，一种所谓的具有平衡的恒化器.