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

Question

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

deSolve R包中的ode（）函数使用其中一个实现的ODE求解器来数值求解常微分方程组。为此，它使用动态调整的集成步骤，该步骤基于每次集成结束时的Local Truncation Error值。用户基本上指定了定义时间步长网格的所需输出时间，其可以具有等于或小于输出时间的步长。

如果我们以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个步骤（在times对象中设置）。

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

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

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

Answer 2

经过一番研究[1,2]：

有两种显式方法不适应时间步骤：euler方法和rk4方法。它们以两种方式实现：

1. 作为一般rk求解器的rkMethod。在这种情况下，通过设置参数hini，可以独立于times参数指定使用的时间步长。函数ode使用此通用代码
2. 作为特殊求解器代码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：package deSolve中求解微分方程。统计软件杂志，33。

[2] Soetaert，K.，Petzoldt，T。，＆Setzer，R.W。（2010）。 Package deSolve：在R. J Stat Softw，33（9），1-25中求解初值微分方程。