R语言问题：用purrr::map结合deSolve求解多参数ODE报错

问题描述 #

我需要基于嵌套数据框test_df_nested中的不同参数批量求解ODE，数据框包含两列：

• PP：核心变量
• data：存储由PP计算得到的14个参数，展开后每行对应一组独立参数

尝试用deSolve的ode函数结合purrr的map2遍历嵌套数据框每行求解ODE，获取变量随时间的变化，但仅当数据框只有1行时可行，多行时报错：ode函数输出元素数多于初始向量y的元素数。推测问题源于参数未随行正确传递，同时还需要处理时间依赖参数（需结合嵌套数据框内的时间参数动态计算），想实现用map迭代带时间依赖逻辑的自定义函数。

附嵌套数据框的dput输出：

structure(list(PP = c(0, 0.1, 0.2), data = list(structure(list(
    u_croiss = 30000, kUpeak = 560000, kUstable = 480000, w_0 = 300, 
    kWpeak = 9700, kWstable = 4500, k_m = 0.84, k_c = 4.74, chi_M = 9.31204630137287e-09, 
    chi_C = 2.91010985906386e-08, t_low = 105, t_kpeak = 155, 
    t_kstable = 255), row.names = c(NA, -1L), class = c("tbl_df", 
"tbl", "data.frame")), structure(list(u_croiss = 42000, kUpeak = 588000, 
    kUstable = 504000, w_0 = 420, kWpeak = 10185, kWstable = 4725, 
    k_m = 0.956, k_c = 5.409, chi_M = 9.24967155645443e-09, chi_C = 2.90483478901307e-08, 
    t_low = 105, t_kpeak = 152.5, t_kstable = 252.5), row.names = c(NA, 
-1L), class = c("tbl_df", "tbl", "data.frame")), structure(list(
    u_croiss = 54000, kUpeak = 616000, kUstable = 528000, w_0 = 540, 
    kWpeak = 10670, kWstable = 4950, k_m = 1.072, k_c = 6.078, 
    chi_M = 9.19322958436019e-09, chi_C = 2.90004488568525e-08, 
    t_low = 105, t_kpeak = 150, t_kstable = 250), row.names = c(NA, 
-1L), class = c("tbl_df", "tbl", "data.frame"))), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -3L), groups = structure(list(
    PP = c(0, 0.1, 0.2), .rows = structure(list(1L, 2L, 3L), ptype = integer(0), class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -3L), .drop = TRUE))

原代码中的make_ODE和微分方程函数：

make_ODE <- function(PP,data){

  Unnested_df <- test_df_nested %>% 
    unnest(data)
  
  yini  <- c(V = 1)
  times <- seq(0, 200, by = 1)
  
  parms  <- c(
    k_V = 100)
  
  actual_ode <- function(yini, times, parms){
    out   <- ode(y = yini,
                 times,
                 equa_diff_sp_test_nest,
                 parms,
                 method = "rk4") %>% 
      as_tibble()
    
    return(out)
    
  }
  
  res <-actual_ode(yini, times, parms)
  return(res)
  
}

equa_diff_sp_test_nest <- function(t,y,parms){
  
  V  <- y[1]

  
  with(as.list(c(y, parms)), {
    

    dVdt <- v_croiss * V * (1 - V/k_V)
    

        return(list(c(dVdt)))
  })
}

错误原因分析 #

1. 全局参数污染：make_ODE中错误地全局读取整个test_df_nested并展开，导致每次调用都拿到所有行的参数，而非当前行的独立参数
2. 参数缺失：微分方程中用到的v_croiss未正确传入参数列表，依赖未定义的全局变量
3. 分组干扰：原数据框是分组状态，可能影响map迭代的上下文

修正方案与代码 #

1. 修正make_ODE函数 #

直接使用传入的当前行参数，避免全局读取，同时将所有必要参数传入ode：

make_ODE <- function(PP, data) {
  # 将当前行的参数子数据框转为命名向量
  current_pars <- data %>% unlist()
  yini <- c(V = 1)
  times <- seq(0, 200, by = 1)
  
  # 合并当前行参数与PP，统一传入ode
  parms <- c(current_pars, PP = PP)
  
  # 直接调用ode求解，简化嵌套逻辑
  out <- ode(y = yini,
             times = times,
             func = equa_diff_sp_test_nest,
             parms = parms,
             method = "rk4") %>% 
    as_tibble()
  
  return(out)
}

2. 修正微分方程函数（支持时间依赖参数） #

根据当前时间t结合参数中的时间阈值动态计算参数：

equa_diff_sp_test_nest <- function(t, y, parms) {
  V <- y[1]
  
  with(as.list(parms), {
    # 示例：根据时间区间动态调整增长速率v_croiss
    if (t < t_low) {
      v_croiss <- u_croiss * 0.5
    } else if (t >= t_low && t < t_kpeak) {
      v_croiss <- u_croiss
    } else {
      v_croiss <- u_croiss * 0.8
    }
    
    # ODE核心方程，这里用kUpeak作为承载能力示例
    dVdt <- v_croiss * V * (1 - V / kUpeak)
    
    return(list(c(dVdt)))
  })
}

3. 执行批量求解 #

先取消数据框分组，再用map2迭代：

# 取消分组避免迭代干扰
test_df_nested <- test_df_nested %>% ungroup()

# 批量求解ODE
test <- test_df_nested %>% 
  mutate(outputs = map2(PP, data, make_ODE))

关键修正点说明 #

• 参数正确传递：make_ODE仅使用当前行的data参数，确保每组参数独立求解
• 时间依赖处理：在微分方程内根据当前时间t和参数中的时间阈值动态计算参数，满足时间依赖需求
• 简化逻辑：移除不必要的嵌套函数，直接调用ode，减少代码层级
• 分组处理：提前取消数据框分组，避免map迭代时的分组上下文干扰

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