如何提取R语言deSolve包ode()函数每步的局部截断误差(LTE)值？

当然可以提取lsode求解过程中的LTE估计值！deSolve的lsode底层实现其实会记录每一步积分的误差信息，只是默认不会直接返回给用户，需要通过特定的参数设置来解锁这些数据。下面结合你的Lotka-Volterra模型示例，一步步展示如何操作：

步骤1：修改ode()调用以启用内部信息记录 #

在调用ode()时，需要添加control参数来开启密集输出和求解器内部数据的返回：

rm(list = ls())
install.packages('deSolve')
library('deSolve')

# Example ODE system for the Lotka-V predator prey model
LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion <- rIng * Prey * Predator
    GrowthPrey <- rGrow * Prey * (1 - Prey/K)
    MortPredator <- rMort * Predator
    dPrey <- GrowthPrey - Ingestion
    dPredator <- Ingestion * assEff - MortPredator
    return(list(c(dPrey, dPredator)))
  })
}

# values of the parameters
pars <- c(rIng = 0.2, # /day, rate of ingestion
          rGrow = 1.0, # /day, growth rate of prey
          rMort = 0.2 , # /day, mortality rate of predator
          assEff = 0.5, # -, assimilation efficiency
          K = 10) # mmol/m3, carrying capacity

#initial values for the predator prey state variables
yini <- c(Prey = 1, Predator = 2)

#times at which the ode return state output
times <- seq(0, 200, by = 1)

# 修改调用：启用lsode的内部数据返回
out <- ode(yini, times, LVmod, pars, method = "lsode",
           control = list(dense = TRUE, return.sol = TRUE))

步骤2：提取LTE估计值 #

求解完成后，返回的out对象会带有一个名为lsode的属性，里面包含了求解过程的核心细节，包括每一步的LTE估计：

# 提取lsode的内部求解信息
lsode_diagnostics <- attr(out, "lsode")

# 查看包含的关键字段
str(lsode_diagnostics)

你会看到lsode_diagnostics里包含：

• h：每一步的积分步长
• err：局部截断误差估计值，每一行对应一个积分步，每一列对应一个状态变量（这里就是Prey和Predator）
• ewt：用于误差判断的权重向量（和rtol、atol相关）

可以把这些数据整理成更易读的数据框：

# 整理LTE和步长为数据框
lte_results <- data.frame(
  Integration_Step = seq_along(lsode_diagnostics$h),
  Step_Size = lsode_diagnostics$h,
  Prey_LTE = lsode_diagnostics$err[, 1],
  Predator_LTE = lsode_diagnostics$err[, 2]
)

# 查看前几行结果
head(lte_results)

关键说明 #

• 求解器会根据lsode_diagnostics$err中的值动态调整步长，确保误差尽可能接近（但略高于）你设置的rtol和atol阈值，这是变步长线性多步方法的正常行为。
• 如果只需要全局误差统计，直接调用summary(out)就能看到，但要获取每一步的LTE细节，必须通过提取lsode属性的方式。

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