パイこね変換のリアプノフ指数をPythonで計算する

離散力学系でカオス的特徴を持つ代表的な例とされているパイこね変換（Baker’s Transformation）をPythonで実装してみた。

パイこね変換とは

パイこね変換（パイこねへんかん、英語: baker’s transformation）とは、2次元の離散力学系の一種で、カオスを生み出す典型的な仕組みを抜き出した基礎的な系として知られる。名称は料理におけるパイの生地を引き延ばして折り畳む操作に因む。（From Wiki)

今回参考にしたのは平山修著「Excelで試す非線形力学」コロナ社で、元々は各期の数値をセルで計算しながらグラフで描画させるような記載になっている。数値の変化が非常にわかりやすくて助かった。

やってみたこと

２つのごくわずかに異なる初期値をパイこね変換の関数に入力し、それぞれの軌道を追っていく。これを推移させ、２つの軌道がi期において初期と比べて差がどれほど開いていくのかの指標がリアプノフ指数。カオス性を持つ時系列データの場合、リアプノフ指数は正の値になることが知られている。

まずは2つの軌道を50期回し、誤差の推移を観察した。コードは下記の通り。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import *
import math
# 2つの初期値を持つパイこね変換（baker’s transformation）
def bakers_transformation_1(y0):
    y = [y0]
    for i in range(50):
        if y[i] < 0.5:
            y_next = y[i] * 2
            y.append(y_next)
        else:
            y_next = 2*(1-y[i])
            y.append(y_next)
    return y
def bakers_transformation_2(z0):
    z = [z0]
    for i in range(50):
        if z[i] < 0.5:
            z_next = z[i] * 2
            z.append(z_next)
    else:
        z_next = 2*(1-z[i])
        z.append(z_next)
    return z
# 初期値
y0 = 0.324827
z0 = y0 + 0.000001
y = bakers_transformation_1(y0)
z = bakers_transformation_2(z0)
# プロット表示
ax = plt.gca()
ax.set_title(‘Baker\’s Transformation’)
ax.xaxis.set_major_locator(MultipleLocator(5))
plt.plot(y, label=”b_t1")
plt.plot(z, label=”b_t2")
plt.legend()
plt.show()

15期ごろまではほぼ差異がないものの、突如大きくずれていく様子が見て取れる。カオスの特徴として初期値鋭敏性、それに伴う長期予測の困難が挙げられますが、それがよく理解できた。

リアプノフ指数の計算

リアプノフ指数の定義は、以前ライブラリでビットコイン価格のカオス性を調べた記事のこちらに記載している。今回は直接微分はしておらず、二つの差異の絶対値をとり、初期値との比率を対数化するという実装をした。

def lyapunov_exponent():
    lyapunov_exponent = []
    d_0 = abs(y[0]-z[0]) #初期値段階での差異の絶対値
    for i in range(int(len(y))):
        d_i = abs(y[i] - z[i]) # i期での差異の絶対値
        d_rate = (d_i/d_0) # i期と初期値のずれ幅の比率
        ln_d_rate = math.log(d_rate) # その対数
        lmd = (1/(i+1)) * ln_d_rate
        lyapunov_exponent.append(lmd)
    return lyapunov_exponent
lyapunov_exponent = lyapunov_exponent()
ax.set_title('Lyapunov Exponent')
plt.plot(lyapunov_exponent, label="Lyapunov Exponent")
plt.legend()
plt.show()

結果はこちら。リアプノフ指数は正の値をとっており、パイこね変換の力学系にはカオス性が認められることが伺えた。