如何在B. Raynor的集合种群模型中引入时间依赖mu参数并解决维度错误？

我想在B. Raynor的集合种群模型中引入随时间变化的mu参数，期望mu取向量mu = seq(from = 0, to = 1, length = 100)的值。原代码可正常运行，但替换mu的定义后报错：

Error in eval(substitute(expr), data, enclos = parent.frame()) :
dims [product 3] do not match the length of object [100]

原代码如下：

#加载库
library(dplyr)
library(deSolve)

MODEL <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
  
    #定义初始条件
    S = matrix(state[1:3], ncol=1)
    I = matrix(state[(4:6)], ncol=1)
    R = matrix(state[(7:9)], ncol=1)
    
    #定义参数
    gamma <- matrix(parameters[paste0("gamma", 1:3)], ncol=1)
    alpha <- matrix(parameters[paste0("alpha", 1:3)], ncol=1)
    mu <- matrix(parameters[paste0("mu", 1:3)], ncol=1)

    #mu = seq(from = 0.001, to = 0.0015, length = 100)
    
    dS <- -mu*S + alpha*R
    dI <-  mu*S - gamma*I 
    dR <-  gamma*I - alpha*R
    
    output <- c(dS, dI, dR)
    list(output)
  })
}

init <- c(S1 = 100,  S2 = 100,  S3 = 100,  I1 = 10,  I2 = 10, I3 = 0, R1 = 0,  R2 = 0,  R3 = 0 )

#参数
parms <- c(gamma1 = 0.2000000, gamma2 = 0.2500000, gamma3 = 0.3333333, mu1 = 0.0010000,
           mu2 = 0.0010000, mu3 = 0.0010000,alpha1 = 0.200000, alpha2 = 0.3000000, alpha3 = 0.4000000)

Time = 100
dt = 1 #步长dt
times <- seq(0, Time, by = dt)

#运行模拟
out <- ode(y=init, times=times, func=MODEL, parms=parms)

1. 维度不匹配：原代码中mu是对应3个种群的3行1列矩阵，直接替换为100个值的向量后，和3行1列的S矩阵相乘时维度无法兼容，触发报错。
2. 逻辑错误：deSolve的ode函数会在每个时间步调用MODEL并传入当前时间点，直接生成全时间序列的mu向量，没有和当前时间点对应，不符合微分方程求解的逻辑。

提前生成与时间序列长度匹配的mu序列，在MODEL函数中根据当前时间点获取对应mu值，再分配给3个种群（以下示例假设3个种群的mu随时间同步变化，若需不同规律可调整）。

修改后的代码：

#加载库
library(dplyr)
library(deSolve)

MODEL <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
  
    #定义初始条件
    S = matrix(state[1:3], ncol=1)
    I = matrix(state[(4:6)], ncol=1)
    R = matrix(state[(7:9)], ncol=1)
    
    #定义参数
    gamma <- matrix(parameters$gamma, ncol=1)
    alpha <- matrix(parameters$alpha, ncol=1)
    
    # 根据当前时间匹配序列索引，获取对应mu值
    time_idx <- which(times == time)
    current_mu <- parameters$mu_seq[time_idx]
    # 将当前mu值分配给3个种群
    mu <- matrix(rep(current_mu, 3), ncol=1)

    dS <- -mu*S + alpha*R
    dI <-  mu*S - gamma*I 
    dR <-  gamma*I - alpha*R
    
    output <- c(dS, dI, dR)
    list(output)
  })
}

init <- c(S1 = 100,  S2 = 100,  S3 = 100,  I1 = 10,  I2 = 10, I3 = 0, R1 = 0,  R2 = 0,  R3 = 0 )

Time = 100
dt = 1 #步长dt
times <- seq(0, Time, by = dt)

# 生成与时间序列长度一致的mu序列
mu_seq <- seq(from = 0, to = 1, length = length(times))

# 参数改为列表格式，方便传入序列
parms <- list(
  gamma = c(0.2000000, 0.2500000, 0.3333333),
  alpha = c(0.200000, 0.3000000, 0.4000000),
  mu_seq = mu_seq
)

#运行模拟
out <- ode(y=init, times=times, func=MODEL, parms=parms)

• 若需要3个种群的mu有不同时间变化规律，可分别生成mu1_seq、mu2_seq、mu3_seq，在MODEL中分别获取对应时间点的值即可。
• 将parms改为列表类型，是为了更方便地传入整个时间序列参数，避免原向量格式的局限性。

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