Rs deSolve 和 Pythons odeint 的区别

Question

我目前正在使用 Python 和 R 探索 Lorenz 系统，并注意到 ode 包中的细微差别。 odeint from Python and ode 都说他们使用 lsoda 来计算他们的导数。但是，对两者都使用 lsoda 命令似乎会产生截然不同的结果。我已经尝试 ode45 为 R 中的 ode 函数得到更类似于 Python 的东西，但我想知道为什么我不能得到完全相同的结果：

from scipy.integrate import odeint
def lorenz(x, t):
    return [
        10 * (x[1] - x[0]),
        x[0] * (28 - x[2]) - x[1],
        x[0] * x[1] - 8 / 3 * x[2],
    ]


dt = 0.001
t_train = np.arange(0, 0.1, dt)
x0_train = [-8, 7, 27]
x_train = odeint(lorenz, x0_train, t_train)


x_train[0:5, :]
array([[-8.        ,  7.        , 27.        ],
       [-7.85082366,  6.98457874, 26.87275343],
       [-7.70328919,  6.96834721, 26.74700467],
       [-7.55738803,  6.95135316, 26.62273959],
       [-7.41311133,  6.93364263, 26.49994363]])

library(deSolve)
n <- round(100, 0)
# Lorenz Parameters: sigma, rho, beta
parameters <- c(s = 10, r = 28, b = 8 / 3)
state <- c(X = -8, Y = 7, Z = 27) # Initial State
# Lorenz Function used to generate Lorenz Derivatives
lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dx <- parameters[1] * (state[2] - state[1])
    dy <- state[1] * (parameters[2] - state[3]) - state[2]
    dz <- state[1] * state[2] - parameters[3] * state[3]
    list(c(dx, dy, dz))
  })
}
times <- seq(0, ((n) - 1) * 0.001, by = 0.001)
# ODE45 used to determine Lorenz Matrix
out <- ode(y = state, times = times,
           func = lorenz, parms = parameters, method = "ode45")[, -1]
out[1:nrow(out), , drop = FALSE]
             X        Y        Z
 [1,] -8.00000000 7.000000 27.00000
 [2,] -7.85082366 6.984579 26.87275
 [3,] -7.70328918 6.968347 26.74700
 [4,] -7.55738803 6.951353 26.62274
 [5,] -7.41311133 6.933643 26.49994

我必须调用 out[1:nrow(out), , drop = FALSE] 才能获得完整提供的小数位，head 似乎四舍五入到最接近的五位。我知道它非常微妙，但希望得到完全相同的结果。有谁知道这是否不仅仅是 R 和 Python 之间的舍入问题？

提前致谢。

Answer 1

所有求解 ODE 的数值方法都是在给定精度范围内工作的近似值。默认情况下，deSolve 求解器的精度设置为 atol=1e-6, rtol=1e-6，其中 atol 是绝对精度，rtol 是相对公差。此外，ode45还有一些额外的参数来微调自动步长算法，它可以利用插值。

增加公差，设置例如：

out <- ode(y = state, times = times, func = lorenz, 
  parms = parameters, method = "ode45", atol = 1e-10, rtol = 1e-10)

最后，我建议使用像 lsoda 或 vode 这样的 odepack 求解器，而不是经典的 ode45。更多详细信息可以在 ode 和 lsoda 帮助页面以及 ode45 的 ?rkMethod 帮助页面中找到。

odeint 也可能存在类似的参数。

最后一点：由于洛伦兹是一个混沌系统，局部误差会因误差放大而导致发散行为。这是混沌系统的一个基本特征，从理论上讲，混沌系统在长期运行上是不可预测的。所以无论你做什么，设置多少精度，模拟的轨迹都不是“真实的”，它们只是表现出相似的模式。

Rs deSolve 和 Pythons odeint 的区别

Differences between Rs deSolve and Pythons odeint

python

r

odeint

desolve