有没有办法强制 ode() [deSolve-R 包] 在 ode 函数中的每个集成步骤提供输出答案

【问题标题】：Is there a way to force ode() [deSolve-R package] to provide output at each integration step in the ode function有没有办法强制 ode() [deSolve-R 包] 在 ode 函数中的每个集成步骤提供输出
【发布时间】：2019-08-21 22:39:36
【问题描述】：

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

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

如果我们以 Lotka Volterra Predator-Prey 模型（有后勤猎物）为例：

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 步（在时间对象中设置）。

对于我正在运行的模型而言，这种差异要大得多，并且要对系统的稳定性进行完整分析，我需要每个集成步骤的输出以及步骤的大小。有没有办法从 ode 中提取这些信息？

【问题讨论】：

您可以在ode() 中使用method = "euler" 作为额外参数，这样您将有200 个步骤（因为您的时间步长为1）。不是一个好的答案，但也许是第一次尝试。

标签： r ode

【解决方案1】：

可以通过在模型函数中放置 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)))
  })
}

这适用于自动和固定步长求解器。但请注意，自动步进器有时可能会丢弃步骤并重试，因此时间并不总是单调的。如果您想保存数据以供以后使用，请使用

【讨论】：

【解决方案2】：

所以，经过一番研究[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)

在这个非常简单的系统中使用显式 Euler 方案可能有意义，但是对于更复杂的系统，建议使用 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 中的初值微分方程。

【讨论】：