dynamical systems

In [26]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib nbagg

II. Dynamical systems

  1. Dynamical systems
  2. Simple simulation example 1
  3. Simple simulation example 2 (Euler's method, lsode)
  4. Chaotic example: Lorenz system (multivariate)
  5. Takens embedding theorem demonstration on Lorenz-system

1. Dynamical systems

1.1. Simple simulation example 1: Logistic map

The basic generating formula:

$ x_{t+1} = r x_t (1-x_t)$

In [27]:
def x_next(x):
    return 4. * x * (1. - x)
In [28]:
n = 200
x = [0.3]
for i in range(n):
    x.append(x_next(x[-1]))
x = np.array(x)
print x.shape
(201,)
In [29]:
plt.figure()
plt.plot(x)

plt.title('Logistic map timeseries')
plt.xlabel(r'$t$')
Out[29]:
<matplotlib.text.Text at 0x7fcb8b4b0c50>
In [30]:
fs = 18
ts = 20

plt.figure()
plt.plot(x[:-1], x[1:], '.')

plt.title('Return map for logistic map', fontsize=ts)
plt.xlabel(r'$x_t$', fontsize=fs)
plt.ylabel(r'$x_{t+1}$', fontsize=fs)
Out[30]:
<matplotlib.text.Text at 0x7fcb8b0ecc90>

2. Simple simulation example 2: Exponentional decay

The differential equation for exponentional decay:

$ \dot{x} = -a x$


$ \frac{1}{x} \mathrm{d}x = -a \mathrm{d}t \\
\int \frac{1}{x} \mathrm{d}x = \int-a \mathrm{d}t \\
\log(x) = - C a t \\ \\
x = x_0 \mathrm{e}^{-at}
$

In [31]:
def dx(x, a=0.1):
    return -a * x
In [32]:
a = 0.1
n = 300
dt = 0.1
t = np.arange(0, (n+1)*dt, dt)
x = [10.]
for i in range(n):
    x_new = x[-1] + dx(x[-1], a=a) * dt
    x.append(x_new)
x = np.array(x)
print x.shape
(301,)
In [33]:
plt.figure()
plt.plot(t, x, lw=2., alpha=0.8, label='Euler method')
plt.plot(t, x[0] * np.exp(-a*t), 'r--', lw=2., label='analytical')

plt.title('Exponentional decay', fontsize=ts)
plt.xlabel(r'$t$', fontsize=fs)
plt.ylabel(r'$x$', fontsize=fs)
plt.legend(loc='upper right')
Out[33]:
<matplotlib.legend.Legend at 0x7fcb8afe0290>

3. Lorenz sysytem

In [34]:
from scipy.integrate import odeint

def Lorenz(state,t):
  # unpack the state vector
  x = state[0]
  y = state[1]
  z = state[2]

  # these are our constants
  sigma = 10.0
  rho = 28.0
  beta = 8.0/3.0

  # compute state derivatives
  xd = sigma * (y-x)
  yd = (rho-z)*x - y
  zd = x*y - beta*z

  # return the state derivatives
  return [xd, yd, zd]

state0 = [2.0, 3.0, 4.0]
t = np.arange(0.0, 30.0, 0.01)

state = odeint(Lorenz, state0, t)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(state[:,0],state[:,1],state[:,2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

Time delay embedding:

In [35]:
d = 3 # Embedding dimesnion
tau = 10  # embedding-delay
fs = 25 # font size for axis labels

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(state[2*tau:, 0],state[tau:-tau, 0],state[:-2*tau, 0])
ax.set_xlabel(r'$x_t$', fontsize=fs)
ax.set_ylabel(r'$x_{t-\tau}$', fontsize=fs)
ax.set_zlabel(r'$x_{t-2\tau}$', fontsize=fs)
plt.show()
In [ ]: