R语言deSolve包ODE模型中mf变量始终为0的问题求助

导致mf恒为0的核心原因 #

• 脉冲触发条件无效：deSolve的ode求解器采用连续时间数值积分，几乎不会精确取到t == 156，这导致app1始终为0，无法给mf注入初始农药量。
• ODE方程逻辑混淆：原方程错误地将离散时间步长（dt）的计算逻辑混入连续ODE模型中，手动使用exp(-kf*dt)这类离散衰减项会导致变化率计算异常，即使有输入也无法累积有效数值。
• 方程项冗余矛盾：原方程中-(mf* exp(-kf*dt)) + (yf * (mf* exp(-kf*dt)))化简后为mf*exp(-kf*dt)*(yf-1)，结合连续ODE的逻辑，这部分的物理意义不清晰，且与自然衰减的连续表达冲突。

修复方案与代码调整 #

1. 使用事件机制处理脉冲操作 #

用deSolve的events参数定义精确时间点的农药施加和收获，避免连续时间求解的匹配问题。

2. 重构ODE为连续时间变化率形式 #

聚焦连续时间下的变化率，让求解器自动处理积分：

• 自然衰减：-kf * mf（连续速率）
• 雨水冲刷：-mf * (1 - exp(-fet * P))（简化后的连续冲刷项，因每日步长为1）

修复后的完整代码 #

library(deSolve)

# 定义连续时间ODE模型
mod <- function(t, yini, param, Precip) {   
  jgrow <- 32
  # 计算覆盖度，修正t的边界判断（原t>168改为t>=169避免重叠）
  igrow <- ifelse(t >= 136 & t < 169, t - 136, ifelse(t >= 169, 0.9, 0))
  cover <- 0.9 * igrow / jgrow
  # 匹配对应天数的降水量
  P <- Precip[round(t)]
  
  with(as.list(c(yini, param)), { 
    # 连续变化率：自然衰减 + 雨水冲刷
    d_MF <- -kf * mf - mf * (1 - exp(-fet * P))
    return(list(c(d_MF))) 
  })
}

# 参数配置
params <- list(drift = 0.1, kf = 0.023, sa = 32.4, fet = 0.2, yf = 0.6)
# 时间序列
t <- seq(1, 365, by = 1)

# 降水量数据
Precip <- rep(0, 365)
Precip[155:205] <- c(0, 0.96, 0.71, 0, 0, 0.08, 14.99, 0.2, 0.2, 0.25, 0, 0, 0.1, 0, 3.66, 0, 0.13, 0, 5.28, 0.03, 0, 0, 3.15, 0.84, 1.83, 0.03, 0, 2.29, 0, 0, 0.15, 1.65, 0, 0.91, 0.18, 2.01, 0, 0, 0, 0, 0, 3.17, 0.1, 0, 1.45, 0.25, 0, 0.03, 0.46, 0.03, 0)

# 初始条件
y0 <- c(mf = 0)

# 定义事件：第156天施加农药，第201天收获清零
events <- list(
  data = data.frame(
    time = c(156, 201),
    # 第156天的农药施加量计算：app1*sa*(1-drift)*cover
    value = c(1.12 * params$sa * (1 - params$drift) * (0.9*(156-136)/32), 0),
    method = c("add", "set")  # add为增加mf，set为设置mf为0
  )
)

# 求解模型
sol <- ode(y = y0, times = t, func = mod, parms = params, Precip = Precip, events = events)

# 查看第156天前后的结果验证
print(sol[155:160, ])

关键说明 #

• 事件机制确保在精确时间点执行农药施加和收获，彻底解决连续时间求解器的时间匹配问题。
• 重构后的ODE方程专注于连续时间变化率，让deSolve自动处理积分过程，保证数值求解的正确性。
• 农药施加量直接在事件中计算，确保覆盖度、沉降效率等参数正确代入。

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