Using different time point sets (parameter t) for python scipy builtin function ODEINT

Question

I am working on the Lorenz system. I used the 'scipy.integrate.odeint' builtin function to integrate as suggested on Wikipedia1. Lorenz system has three variables: x, y, z. When I compared the evolution of x for two Lorenz systems with same initial conditions, same time differential(dt), but different time set points, I obtained different sets. Both the systems evolved similarly for some time, but later on, diverged 2. What is the reason? The only difference is they have different time sets (but with same time differential). How is the time point set influencing the integration?

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])))`

I have also attached a graph I obtained for the above code: the evolution of x in these two systems.

score 0 · Accepted Answer · answered Dec 12 '19 at 12:40

What you observe is perfectly normal.

In principle, as you guessed, there should be no difference. The trick here is that by storing the time points as floating point values (you have no choice of course), the intervals are not perfectly equal. To observe this, consider an extra time array t2_shifted whose origin is set back to zero:

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

Now, plot the difference:

plt.plot(t - t2_shifted)

There are small differences, of the order of the machine precision (the difference is up to about 3 10^{-14} on my computer, x86_64 architecture). This difference in the timesteps will generate tiny differences in the trajectories.

As the Lorenz system is chaotic, the tiny difference in the evolution of the system will lead to a divergence of the trajectories.

To observe this divergence phenomenon, you can plot the x coordinate (or y or z of course) for the two solutions. They are identical at first, then diverge a bit, then are completely disconnected as is expected for chaotic systems. The attractor will be the same in both cases, however.

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

Using different time point sets (parameter t) for python scipy builtin function ODEINT

1 Answers1