Rs deSolve 和 Pythons odeint 的区别

Posted 2023-02-18

【中文标题】Rs deSolve 和 Pythons odeint 的区别

【英文标题】：Differences between Rs deSolve and Pythons odeint

【发布时间】：2021-08-09 22:12:14

【问题描述】：

我目前正在使用 Python 和 R 探索 Lorenz 系统，并注意到 ode 包中的细微差别。 odeint 来自Python 和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 之间的舍入问题？

提前致谢。

【参考方案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 帮助页面以及 ?rkMethod 帮助页面中的 ode45。

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

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

【讨论】：

谢谢，我能够接受您的建议，并为这两个软件包使用 lsoda，并通过将公差数字调整为 odeint 的默认值来重现 X 变量的相同结果.但是，为什么Y 和Z 仍然给出不同的结果？

是odeint 使用与ode() 相同的容差？

我根据odeint 公差更改了公差，即1.49012e-8

ode45 的实现在 Python 和 R 之间有所不同，它们只是使用相同的 Runge-Kutta 对（Dormand-Prince 4、5）。相比之下，lsoda 是 odepack 算法，但它也针对 R 进行了一些修改以允许附加功能。不确定这是否会改变“精确”公差，但这通常是不感兴趣的，因为所有求解器都是近似值。好消息是，deSolve 中的所有求解器都经过了严格的基准测试，请参阅包 deTestSet 中的示例。

是的，由于四舍五入、算法的细微差异以及（最重要的）混沌系统的误差放大。