有没有办法强制 ode() [deSolve-R package] 在 ode 函数的每个积分步骤提供输出

Question

我想在数值方案的每一步提取状态变量的值。

deSolve R 包中的 ode() 函数使用其中一个已实现的 ODE 求解器对常微分方程组进行数值求解。为此，它根据每次积分结束时的局部截断误差值使用动态调整的积分步骤。用户基本上指定了定义时间步长网格的所需输出时间，其步长可以等于或小于输出时间。

如果我们以Lotka Volterra Predator-Prey模型（带logistic猎物）为例：

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)))
  })
}

参数和状态变量定义为：


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



yini  <- c(Prey = 1, Predator = 2)

并在每日时间步请求输出：

times <- seq(0, 200, by = 1)

out   <- ode(yini, times, LVmod, pars)
diagnostics(out)

查看诊断，我们可以看到求解器总共使用了 282 个步骤，而生成的输出为 200 个步骤（在 times 对象中设置）。

这个差异对于我运行的模型来说要大得多，为了对我的系统的稳定性进行完整的分析，我需要每个集成步骤的输出以及步骤的大小。有没有办法从 ode 中提取此信息？

Answer 1

所以，经过一些研究 [1,2]：

不适应时间步长的显式方法有两种：欧拉法和rk4法。它们以两种方式实现：

作为通用rk求解器的rkMethod。在这种情况下，可以通过设置参数 hini 独立于时间参数指定使用的时间步长。功能ode使用这个通用代码
作为特殊求解器代码 euler 和 rk4。这些实现被简化并且具有更少的选项以避免开销。使用的时间步长由 times 参数中的时间增量决定。

应用于 LV 示例，接下来的两个语句都触发 Euler 方法，第一个使用时间步长 = 1 的“特殊”代码，由 times 参数强加，第二个使用通用方法hini 设置的时间步长。


    out.euler  <- euler(y = state, times = times, func = LVmod, parms = parameters)
    out.rk <- ode(y = state, times = times, func = LVmod, parms = parameters,
                 method = "euler", hini = 0.01)

在这个非常简单的系统中，使用显式欧拉方案可能是有意义的，但是对于更复杂的系统，建议使用 rk4，即使它仍然是一个显式方案。

总结：

似乎没有办法提取每一步的状态值在 ode() 函数中使用动态时间步时。
Euler 和 rk4 这两个显式方案可以设置常量步长。

[1] Soetaert, K. E. R., Petzoldt, T., & Setzer, R. W. (2010)。在 R 中求解微分方程：包 deSolve。统计软件杂志，33。

[2] Soetaert, K., Petzoldt, T., & Setzer, R. W. (2010)。 deSolve 包：在 R. J Stat Softw 中求解初始值微分方程，33(9), 1-25.

Answer 2

可以通过在模型函数中放置 print 或 cat 来观察积分器的工作方式：

library("deSolve")

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
    cat("Time=", Time, "dPrey"=dPrey, "dPredator=", dPredator, "\n")
    return(list(c(dPrey, dPredator)))
  })
}

这适用于自动和固定步长求解器。但请注意，自动步进器有时可能会放弃步骤并再次尝试，因此时间并不总是单调的。如果要保存数据供以后使用，请将数据保存到带有 <<- 的列表中，最好是围绕模拟构建闭包。

有没有办法强制 ode() [deSolve-R package] 在 ode 函数的每个积分步骤提供输出

Is there a way to force ode() [deSolve-R package] to provide output at each integration step in the ode function

r

ode