レスラーダイナミカルシステム (Rössler Dynamical Sysmtem)のシミュレーション

はじめに

「統計科学のフロンティア 4　階層ベイズモデルとその周辺―時系列・画像・認知への応用―」の「第II部　非線形ダイナミカルシステムの再構成と予測」を読んでいた際に，分析対象の時系列データとして，レスラー(Rössler)系のダイナミカルシステムが出てきたので，そのシミュレーション方法を調べた。

www.iwanami.co.jp

問題設定

本書のテーマは，「非線形の時系列モデルを，階層ベイズモデルの枠組みを使って予測する」というものである。
この非線形時系列モデルの例として，レスラー(Rössler)系のダイナミカルシステムが出てくる。

レスラー(Rössler)系の微分方程式

レスラー(Rössler)系の微分方程式の微分方程式は，以下のように表される。

$\begin{align} \dot{x} &= -y - z \\ \dot{y} &= x + ay \\ \dot{z} &= bx - cz + xz \\ \end{align}$

ただし本書中では， $a = 0.36, b=0.4, c = 4.5$ としている。
上式のように， $xz$ の項が存在するので，非線形な微分方程式であることが分かる。

レスラー(Rössler)系の微分方程式の差分化

数値計算などをしやすくするために，次のように差分化する。

$\begin{align} x_{(t+1)\delta} &= -y_{t\delta} - z_{t\delta} \\ y_{(t+1)\delta} &= x_{t\delta} + ay_{t\delta} \\ z_{(t+1)\delta} &= bx_{t\delta} - cz_{t\delta} + x_{t\delta} z_{t\delta} \\ \end{align}$

なお，次のステップにおける値を求めるためには，ルンゲ・クッタ(Runge-Kutta)法を用いて数値積分している。

また本書では，具体的なデータの生成方法として以下のようにしている。

ステップ幅 $\delta$ は $\delta = 0.01$ としている。
サンプリング周期 $\eta$ は $\eta=70$ としている。
サンプリング後の学習データのサンプル数は500点である。

シミュレーションのサンプルコード

上記のレスラー(Rössler)ダイナミカルシステムのサンプルコードは以下のようになる。

コード

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

def runge_kutta_4_vector(f, y0, t0, t_end, h):
    """
    4次のルンゲ・クッタ法でベクトルの微分方程式を解く関数
    f: 微分方程式 dy/dt = f(t, y), yはベクトル
    y0: 初期値ベクトル y(t0)
    t0: 開始時刻
    t_end: 終了時刻
    h: 時間ステップ幅
    """
    n = int((t_end - t0) / h)  # ステップ数
    t_values = np.linspace(t0, t_end, n)
    y_values = np.zeros((n, len(y0)))  # 各時刻のベクトル値を格納する配列

    y = np.array(y0, dtype=np.float64)
    for i, t in enumerate(t_values):
        y_values[i] = y
        k1 = np.array(f(t, y))
        k2 = np.array(f(t + h / 2, y + h * k1 / 2))
        k3 = np.array(f(t + h / 2, y + h * k2 / 2))
        k4 = np.array(f(t + h, y + h * k3))
        y += h * (k1 + 2 * k2 + 2 * k3 + k4) / 6

    return t_values, y_values

# レスラー方程式
def rossler(t, y_, a=0.36, b=0.4, c=4.5):
    x, y, z = y_
    dxdt = -y - z
    dydt = x + a * y
    dzdt = b + z * (x - c)
    return [dxdt, dydt, dzdt]

# 初期条件とパラメータ
y0 = [1.0, 0.0, 0.0]  # 初期値 [x, y, z]

N_train = 500  # 学習用データのサンプル数
eta = 70 # サンプリング周期
h = 0.01  # 時間ステップ幅
t0 = 0.0  # 開始時刻
t_end = N_train * eta * h  # 終了時刻

# 解を計算
t_values, y_values = runge_kutta_4_vector(lambda t, y_: rossler(t, y_), y0, t0, t_end, h)

# サンプリング周期で間引く
t_obs = t_values[::eta]
y_obs = y_values[::eta]

実装上のポイント

実装上のポイントとして挙げられるのは，runge_kutta_4_vector(f, y0, t0, t_end, h) 関数は，関数fを引数に取る(汎関数になっている)。こうすることで，数値積分したい微分方程式が別のものになってもルンゲ・クッタ法を適用することができるようになっている。

レスラー系の微分方程式 rossler(t, y_, a, b, c) は，時間の変数tと状態量の変数y_を含んでいる。そのため runge_kutta_4_vector 関数の引数にそのまま入れてしまうと「t, y_が定義されていない」というエラーが出るので，ラムダ式を用いて変数t, y_を定義してから代入している。

可視化結果

間引きをしない場合の結果を可視化すると，以下のようになった。

コード

# 間引き無しデータの3次元プロット
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot(y_values[:, 0], y_values[:, 1], y_values[:, 2], lw=0.5, color='b')
ax.set_title("Rössler Attractor", fontsize=14)
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_zlabel("z", fontsize=12)
plt.show()