R语言使用deSolve遍历参数求解消费者-资源模型平衡值

核心问题梳理 #

• 参数赋值错误：你要遍历的资源增长率g应该取param.values[i]，而非你初始化为全0的variable[i]
• 结果存储逻辑错误：每次循环中sol会被新的迭代结果覆盖，循环结束后仅保留最后一次运行的结果，无法对应所有参数值
• 微分方程内不应直接引用外部循环变量i，你已经将g传入参数列表，直接调用g即可
• 打印语句print(sol[i])索引逻辑错误，每次sol只有2个值（r和n的终值），i超过2时会返回NA

修正后的完整代码 #

# 需提前安装并加载deSolve包
library(deSolve)

# 定义参数序列：资源增长率g的取值
param.values = seq(from = 1, to = 10, by = 1)
# 初始化结果存储矩阵，每行对应一个g值，两列分别存储r和n的平衡态
sol_matrix = matrix(NA, nrow = length(param.values), ncol = 2,
                    dimnames = list(NULL, c("r_eq", "n_eq")))

for (i in 1:length(param.values)){ 
  state <- c(r = 1, n = 1)
  parameters = c(g = param.values[i],# 替换为你定义的参数序列值
                 d = 0.5, 
                 k = 5, 
                 c = 1, 
                 e = 1, 
                 h = 1 
  )
  function1 <- function(times, state, parameters) {
    with(as.list(c(state, parameters)),{
      dr = g*r*(1 - (r/k)) - (c*n*r/(1+(h*c*r))) # 直接用传入的参数g即可
      dn = (e*c*n*r/(1+(h*c*r)))- n*d
      list(c(dr, dn))
    }) 
  }
  times <- seq(0, 100, by = 1)
  
  out <- ode(y = state, times = times, func = function1, parms = parameters)
  
  # 把当前迭代的结果存到矩阵的第i行，用nrow(out)取最后一行比硬写101更稳妥
  sol_matrix[i, ] <- out[nrow(out), 2:3]
}

# 平衡态关系图
plot(sol_matrix[,1] ~ param.values, xlab = "资源增长率g", ylab = "资源平衡态r", type = "b")
plot(sol_matrix[,2] ~ param.values, xlab = "资源增长率g", ylab = "消费者平衡态n", type = "b")

核心索引逻辑说明 #

提前创建和参数序列长度一致的结果存储矩阵sol_matrix，每次循环迭代时，将当前运行得到的终值赋值给sol_matrix[i, ]，用i对应参数序列的位置，循环结束后sol_matrix的每一行就和param.values的每个g值一一对应，直接画图即可。

内容的提问来源于stack exchange，提问作者Clara