python - Rs deSolve 和 Pythons odeint 之间的差异

Question

我目前正在使用 Python 探索洛伦兹系统和 R并注意到 ode 中的细微差别包。 odeint 来自 Python和 ode 都说他们使用 lsoda来计算它们的导数。但是，使用 lsoda两者的命令似乎给出了截然不同的结果。我试过 ode45为 ode函数在 R获得更类似于 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 这样的 odepack 求解器。或 vode而不是经典 ode45 .更多详情请见 ode和 lsoda帮助页面和 ode45在 ?rkMethod 帮助页面中。 odeint 也可能存在类似的参数.
最后一点:由于 Lorenz 是一个混沌系统，由于误差放大，局部误差将导致发散行为。这是混沌系统的一个基本特征，从理论上讲，从长远来看，混沌系统是不可预测的。所以无论你做什么，设置多少精度，模拟的轨迹都不是“真实的”，它们只是显示出类似的模式。