稻田农药ODE模型：dw_t为零时mw无法置零的问题排查

为什么直接赋值会残留极小值？ #

1. 浮点精度限制：数值求解器计算时，dw_t几乎不可能严格等于0，总会有极小的舍入误差（比如1e-16量级），这时候dw_t > 0依然为TRUE，导致mw不会被置零。
2. 求解器不处理非连续点：(mw - sol_diff) * (dw_t > 0)相当于给系统引入了状态突变，但ode45这类自适应步长求解器默认只处理连续光滑的微分方程，不会主动检测这种不连续的突变点，会跳过该时刻继续积分，最终留下残留值。

为什么Events机制没正确计算mout？ #

大概率是两个问题：

• 触发条件错误：用了dw_t == 0直接判断，而没有用**根触发（root-finding）**来捕捉dw_t从正变负的交叉时刻，导致事件没在正确时机触发。
• 事件逻辑错误：在置零mw前，没有将当前mw的数值累加到mout中，直接置零导致这部分农药质量的流出量未被记录。

正确实现方法 #

使用deSolve的根触发事件来精确捕捉dw_t变为0的时刻，并在事件中完成mw置零和mout的累加：

1. 定义微分方程函数 #

去掉原有的条件赋值，只写连续状态下的微分方程：

model <- function(t, state, parms) {
  with(as.list(c(state, parms)), {
    # 替换为你模型中的水深变化率方程
    d_dw_t <- -0.1 * dw_t  # 示例：水深随时间衰减
    
    # 替换为你模型中的水体农药质量变化率方程
    d_mw <- -0.05 * mw - 0.01 * dw_t * mw  # 示例：降解+随水流失
    
    # 正常情况下的农药流出量（不含dw_t=0时的全部流出）
    d_mout <- 0.01 * dw_t * mw
    
    return(list(c(d_dw_t, d_mw, d_mout)))
  })
}

2. 定义根函数（触发事件的条件） #

根函数返回dw_t，当dw_t=0时触发事件，求解器会精确找到这个交叉点：

root_func <- function(t, state, parms) {
  return(state["dw_t"])
}

3. 定义事件处理函数 #

先将当前mw值累加到mout，再把mw置为0：

event_func <- function(t, state, parms) {
  # 把剩余的全部农药质量计入流出量
  state["mout"] <- state["mout"] + state["mw"]
  # 强制置零水体农药质量
  state["mw"] <- 0
  return(state)
}

4. 调用求解器 #

指定rootfun和events参数，启用根触发事件：

library(deSolve)

# 初始状态：水深10cm，水体农药50g，流出量0
init_state <- c(dw_t = 10, mw = 50, mout = 0)
# 模型参数（根据你的需求替换）
parms <- c(k_decay = 0.05, k_runoff = 0.01)
# 时间序列
times <- seq(0, 100, by = 1)

# 求解微分方程
sol <- ode(y = init_state, times = times, func = model, parms = parms,
           method = "ode45",
           rootfun = root_func,
           events = list(func = event_func, root = TRUE))

# 查看结果
head(sol)

关键说明 #

• 根触发事件会让求解器自动调整步长，精确找到dw_t变为0的时刻，避免浮点精度问题。
• 事件处理函数中必须先更新mout再置零mw，确保所有残留农药都被计入流出量。
• 不要在微分方程中保留条件赋值，否则会与事件机制的状态突变逻辑冲突，导致计算错误。

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