Различия между Rs deSolve и Pythons odeint

Question

В настоящее время я изучаю систему Лоренца с помощью и и заметил тонкие различия в пакетах. odeint из и оба говорят, что используют их для расчета своих производных. Однако с помощью lsodaкоманда для обоих, кажется, дает совершенно разные результаты. я пытался ode45 для ode функция, чтобы получить что-то более похожее, но мне интересно, почему я не могу получить точно такие же результаты:

      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относительная терпимость. Кроме того, он имеет некоторые дополнительные параметры для точной настройки алгоритма автоматического размера шага и может использовать интерполяцию.

Чтобы увеличить допуск, установите, например:

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

Наконец, я бы порекомендовал использовать решатель odepack, например или vodeвместо классического. Более подробную информацию можно найти в ode а также lsoda страницы помощи и для ode45 на странице справки?rkMethod.

Подобные параметры могут также существовать для odeint.

Последнее замечание: поскольку Лоренц представляет собой хаотическую систему, локальные ошибки приведут к расхождению в поведении из-за увеличения ошибки. Это существенная особенность хаотических систем, которые теоретически непредсказуемы в долгосрочной перспективе. Итак, что бы вы ни делали и какую точность вы устанавливали, смоделированные траектории не являются «настоящими», они просто показывают похожую картину.