Цикл по параметрам для достижения равновесия с 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 srored
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.

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