Использование разных наборов моментов времени (параметр t) для встроенной функции python scipy ODEINT

Question

Я работаю по системе Лоренца. Я использовал встроенную функцию scipy.integrate.odeint для интеграции, как это было предложено в Википедии 1. Система Лоренца имеет три переменные: x, y, z. Когда я сравнил эволюцию x для двух систем Лоренца с одинаковыми начальными условиями, одинаковым дифференциалом времени (dt), но разными заданными точками времени, я получил разные наборы. Обе системы какое-то время развивались одинаково, но позже разошлись 2. Какова причина? Единственная разница в том, что у них разные наборы времени (но с одинаковой разницей во времени). Как установленный момент времени влияет на интеграцию?

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D

rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0

def f(state, t):
    x, y, z = state  # Unpack the state vector
    return sigma * (y - x), x * (rho - z) - y, x * y - beta * z  # Derivatives
##Evolution_1
state0 = [3.0, 1.0, 1.0]
t = np.arange(0.0, 40.0, 0.01)

states = odeint(f, state0, t)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(states[:, 0], states[:, 1], states[:, 2])
plt.show()

#Evolution_2
state0 = [3.0, 1.0, 1.0]
t2 = np.arange(1.0, 41.0, 0.01)

states2 = odeint(f, state0, t2)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(states2[:, 0], states2[:, 1], states2[:, 2])
plt.show()

plt.plot(states2[:, 0],range(len(states2[:, 0])),states[:, 0],range(len(states2[:, 0])),np.absolute(np.subtract(states2[:, 0],states[:, 0])),range(len(states2[:, 0])))`

Я также приложил график, полученный мной для приведенного выше кода: эволюция x в этих двух системах.

Answer 1

То, что вы наблюдаете, совершенно нормально.

В принципе, как вы уже догадались, разницы быть не должно. Хитрость здесь в том, что при сохранении временных точек в виде значений с плавающей запятой (у вас, конечно, нет выбора), интервалы не полностью равны. Чтобы наблюдать это, рассмотрим дополнительный временной массивt2_shifted чье начало установлено обратно в ноль:

t = np.arange(0.0, 40.0, 0.01)
t2 = np.arange(1.0, 41.0, 0.01)
t2_shifted = t2 - t2[0]

Теперь нарисуйте разницу:

plt.plot(t - t2_shifted)

Есть небольшие отличия порядка машинной точности (на моем компьютере разница составляет примерно 3 10^{-14}, архитектура x86_64). Эта разница во временных шагах приведет к крошечным различиям в траекториях.

Поскольку система Лоренца хаотична, крошечная разница в эволюции системы приведет к расхождению траекторий.

Чтобы наблюдать это явление расхождения, вы можете построить координату x (или, конечно, y или z) для двух решений. Сначала они идентичны, затем немного расходятся, а затем полностью разъединяются, как и положено хаотическим системам. Однако в обоих случаях аттрактор будет одинаковым.

plot(states[:,0])
plot(states2[:,0])