Rs deSolve与Pythons odeint的区别

Question

我目前正在使用Python和R探索Lorenz系统，并注意到ode包中的细微差别。Python和ode的odeint都表示，他们使用lsoda来计算衍生品。但是，对这两个命令使用lsoda命令似乎会产生截然不同的结果。我对R中的ode函数尝试了ode45，以获得更类似于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

所有求解常微分方程的数值方法都是达到给定精度的近似值。默认情况下，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的详细信息。

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

最后注意:由于Lorenz是一个混沌系统，由于误差放大，局部误差将导致发散行为。这是混沌系统的一个基本特征，从理论上讲，从长远来看，混沌系统是不可预测的。因此，无论你做什么，以及你设置了多高的精度，模拟的轨迹都不是“真实的”，它们只是显示了类似的模式。