-2

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

这些是我的初始值 在此处输入图像描述

我将它们初始化为零，因为我不想要任何种类的微生物，只是为了检查没有微生物的化学物质会如何变化。然而，经过 5 或 6 步后，我开始在一些微生物变量中看到不同于零的值： 在此处输入图像描述。

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

4

1 回答 1

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在浮点运算的限制范围内运行良好。
  • 模型应该经过仔细设计和测试。
于 2021-09-06T20:19:48.750 回答