问题描述
我目前正在使用 Python 和 R 探索洛伦兹系统，并注意到 ode 包中的细微差别。 Python 和 odeint 中的 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 之间的舍入问题吗？
提前致谢。
解决方法
所有求解 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 也可能存在类似的参数。
最后一点：由于洛伦兹是一个混沌系统，局部误差会因误差放大而导致发散行为。这是混沌系统的一个基本特征，从理论上讲，从长远来看，混沌系统是不可预测的。因此，无论您做什么，设置多少精度，模拟轨迹都不是“真实轨迹”，它们只是显示出类似的模式。