r - 使用 deSolve R 包中的 ode() 函数时，在每个集成步骤中提取局部截断误差 (LTE) 的值

Question

deSolve R 包中的 ode() 函数使用的默认数值方案是 lsode 方法，它实现了 BDF 和隐式 Adams 线性多步方案来求解 ODE 系统。积分以可变步长完成，该步长通过估计每一步的局部截断误差 (LTE) 来控制。在实践中，步长的调整方式使得 LTE 的估计值小于某个预设值 rtol 和 atol。ode 方法中的默认值为：

ode(rtol = 1e-6, atol = 1e-6)

并且可以适应。

然而，实际评估的 LTE 总是高于这些值。有没有办法从求解器中提取这个值？

下面给出了一个示例系统：

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)


# the solver using lsode
out   <- ode(yini, times, LVmod, pars)

summary(out)

str(out)
 deSolve [1:201, 1:3] 0 1 2 3 4 5 6 7 8 9 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:3] "time" "Prey" "Predator"
 - attr(*, "istate")= int [1:21] 2 282 517 NA 1 1 0 52 22 NA ...
 - attr(*, "rstate")= num [1:5] 1 1 201 0 143
 - attr(*, "lengthvar")= int 2
 - attr(*, "type")= chr "lsoda"

标签： rodenumerical-integrationlte