lsoda 创造虚构价值？答案

【问题标题】：lsoda creating fictional values?lsoda 创造虚构价值？
【发布时间】：2021-10-30 23:52:04
【问题描述】：

我目前正在开发一个集成微生物学和地球化学的模型，该模型使用 lsoda 来求解大量微分方程。该模型太大了，无法在此处发布，因为它由多个模块组成，但我发生了一些非常奇怪的事情。

这些是我的初始值 enter image description here

我已将它们初始化为零，因为我不想要任何种类的微生物，只是为了检查没有微生物时化学会如何变化。然而，经过 5 或 6 步后，我开始在一些微生物变量中看到不为零的值： enter image description here.

我想知道 lsoda 是否正在做某种回合，这就是我得到这些值的原因，否则我无法解释这些值是从哪里冒出来的。如果是这种情况，有谁知道如何阻止这种围捕？

【问题讨论】：

欢迎来到 Stack Overflow。请不要使用代码或数据的图像，因为没有大量不必要的努力就无法使用它们。如果您的问题是可重现的，这真的很有帮助。 How to ask a good question
我很欣赏你的模型很复杂，但如果不看更多细节，基本上是不可能回答这个问题的。浮点数学本质上是不精确的；虽然（例如）乘以 0 应该总是给出零，并且添加 0 永远不应该改变一个值，但大多数其他计算将导致舍入错误（例如 sqrt(2)^2-2 !=0 和其他更微妙的例子......）可能值得尝试欧拉集成（ode(..., method="euler") 以尽量减少其他地方发生任何复杂事情的可能性。
如果您收紧容错，请检查该行为是否会发生变化。绝对公差与组件的比例是否兼容？沿这些点生成 ODE 函数的结果，预期为零的非零分量可能暗示一些编程错误。
请提供足够的代码，以便其他人更好地理解或重现问题。

标签： r differential-equations desolve

【解决方案1】：

简答

不，lsoda 不会创建“虚构”值。

详细解答

正如@Ben Bolker 所说，“浮点运算本质上是不精确的”。我想补充一点，包括lsoda 在内的所有求解器都会计算近似值，并且非线性模型会放大误差。

您可以通过单独调用模型函数（不使用求解器）检查您的模型，如果所有导数都完全为零，请参见以下 Lorenz 系统示例，其中所有状态都设置为零：

library(deSolve)

Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <-  a * X + Y * Z
    dY <-  b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}

parameters <- c(a = -8/3, b = -10, c = 28)
state      <- c(X = 0, Y = 0, Z = 0)
times      <- seq(0, 500, by = 1)

range(Lorenz(0, state, parameters))

给了

[1] 0 0

反例

如果我们修改模型，使得一个导数由于舍入误差而与零略有不同，例如通过使用 Ben 的例子：

Fancy  <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <-  a * X + Y * Z + 27 * (sqrt(X + Y + Z + 2)^2 - 2)
    dY <-  b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}
range(Fancy(0, state, parameters))

然后它给出：

[1] 0.000000e+00 1.199041e-14

根据公差，这样的模型可能会波动甚至崩溃。令人惊讶的是，较小的公差有时可能会更糟：

out1 <- ode(y = state, times = times, func = Lorenz, parms = parameters, atol=1e-6)
out2 <- ode(y = state, times = times, func = Fancy, parms = parameters, atol=1e-6)
out3 <- ode(y = state, times = times, func = Fancy, parms = parameters, atol=1e-8)

par(mfrow=c(3,1))
plot(out1[,1:2], type="l", main="Lorenz, derivatives = 0")
plot(out2[,1:2], type="l", main="Fancy, small error in derivatives, atol=1e-6")
plot(out3[,1:2], type="l", main="Fancy, small error in derivatives, atol=1e-8")

结论

lsoda在浮点运算的限制下工作得很好。
应仔细设计和测试模型。

【讨论】：