使用deSolve迭代法求解差分方程始终返回初始值的问题

使用deSolve的method="iteration"时，模型函数需要返回下一个时间步的状态值，而非微分方程的导数。你的代码存在几个关键问题导致输出始终为初始值：

核心问题点 #

1. 模型函数逻辑错误：你用了类似微分方程的dXXX命名，但iteration方法要求返回离散时间的下一个状态，而非导数。
2. 初始状态不合理：Village = 0会导致所有依赖Village的计算结果为0，状态无法更新。
3. 状态变量更新逻辑混乱：部分变量的计算没有基于当前状态正确推导下一步值，且返回的状态顺序与初始状态的映射关系不清晰。

修正后的代码 #

library(deSolve)

travelmodel <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    # 提取当前时刻的状态变量，避免计算时覆盖
    current_ImVillage <- ImVillage
    current_ImField <- ImField
    current_SmField <- SmField
    current_Village <- Village
    
    # 计算下一个时刻的各状态变量
    next_ImVillage <- Im * current_Village * (1 - a) + current_ImField
    next_ImField <- Im * current_Village * a + Ih * (1 - Im) * current_Village * a
    next_ImFlying <- next_ImField + Im * current_Village * a + Ih * (1 - Im) * current_Village * a
    next_ImStaying <- 1 - next_ImFlying
    
    next_SmVillage <- (1 - Im) * current_Village * (1 - a) + current_SmField
    next_SmField <- (1 - Ih) * (1 - Im) * current_Village * a
    next_SmFlying <- next_SmField + (1 - Ih) * (1 - Im) * current_Village * a
    next_SmStaying <- 1 - next_ImFlying  # 若此处逻辑需调整，请根据模型定义修改
    
    next_Village <- next_ImField + Im * current_Village * a + Ih * (1 - Im) * current_Village * a + 
                    next_SmField + (1 - Ih) * (1 - Im) * current_Village * a
    
    # 返回下一时刻的状态向量，顺序必须与初始状态完全一致
    list(c(next_ImVillage, next_ImField, next_ImFlying, next_ImStaying,
           next_SmVillage, next_SmField, next_SmFlying, next_SmStaying,
           next_Village))
  })
}  

parmstrav <- c(Ih = 0.5, 
               Im = 0.5, 
               a = 0.15)

times <- seq(0, 365, by = 1)   
# 修正Village初始值为非0（示例设为1，需根据模型实际含义调整）
statetrav <- c(ImVillage = 0.9,
               ImField = 0.1,
               ImFlying = 0,
               ImStaying = 0,
               SmVillage = 0.9,
               SmField = 0.1,
               SmFlying = 0,
               SmStaying = 0,
               Village = 1)     

difference_output <- ode(y = statetrav, times = times, func = travelmodel, 
                         parms = parmstrav, method = "iteration")

# 查看前几行结果验证
head(difference_output)

关键修正说明 #

• 明确状态更新逻辑：用current_XXX提取当前状态，next_XXX计算下一时刻状态，避免计算过程中变量被覆盖导致逻辑错误。
• 修正初始值：将Village初始值设为非0，确保所有依赖该变量的计算能正常触发状态变化。
• 保持状态顺序一致：返回的状态向量顺序必须和statetrav的顺序完全匹配，否则会出现变量映射错位。
• 移除导数式命名：放弃dXXX的命名方式，改用next_XXX明确表示是下一时刻的状态值，符合离散模型的逻辑。

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