Lyapunov指数とは

Lyapunov指数は、力学系において初期値が少し異なる二つの軌道がどの程度の速さで離れていくかをあらわす指標です。特に、ある軌道が正のLyapunov指数を持つ場合、その軌道がカオスであることを示唆するため、系を特徴づける重要な量のひとつになっています。Lyapunov指数を解析的に計算することは通常困難で、多くの場合数値計算をすることになります。その際、1次元離散力学系のLyapunov指数は比較的容易に計算できて多くの教科書でも計算例が例示されていますが、連続力学系の場合はいくぶんかの考慮を要します。また、多自由度の力学系においては自由度ごとに異なる速さで軌道が離れて（もしくは近づいて）いくため、系の自由度と同じ数だけLyapunov指数があります。これらはLyapunovスペクトラムと呼ばれています。当然ながら、Lyapunovスペクトラムの計算は最大Lyapunov指数のみの計算に比べて計算量が多く、実用に耐える程度の時間で計算するには工夫が必要となります。

この記事では3次元連続力学系のLorenz方程式を例として、連続力学系のLyapunov指数を数値計算する方法を説明します。

離散力学系の最大Lyapunov指数の数値計算

まずはじめに、離散力学系の最大Lyapunov指数の定義を確認します。

$$x_{i+1} = f(x_i) \tag{1}$$

という一次元写像を考えます。この力学系の軌道上の点$x_n$における関数$f$の微分の絶対値$|f'(x_n)|$は、点$x_n$とその近傍の点との距離が一回の写像によってどの程度離れる（もしくは近づく）かをあらわしています。つまり、点$x_n$とその近傍の点との距離は写像を繰り返すたびにおおよそ$|f'(x_n)|$倍ずつ変化していくことになります。下式のように、軌道 $\lbrace x_1, x_2, \cdots,x_n\rbrace$に沿ってこの拡大縮小率$|f'(x_n)|$を相乗平均したものが最大Lyapunov数$L$です。

$$L = \lim_{n \to \infty}\prod_{i=1}^{n}|f'(x_i)|^{\frac{1}{n}} \tag{2}$$

数値計算の際に総積よりも総和の方が都合がよいため、最大Lyapunov数$L$の対数をとり、最大Lyapunov指数$\lambda$を

$$\lambda = \lim_{n \to \infty}\frac{1}{n}\sum_{i=1}^{n}\ln|f'(x_i)| \tag{3}$$

と定義します。たとえば、$\lambda=\ln2$であれば、一回の写像で軌道間の距離が平均2倍になることを意味しています。

連続力学系の最大Lyapunov指数

離散力学系の最大Lyapunov指数は「写像の反復一回あたりの軌道の平均拡大縮小率」でした。これを連続力学系の場合は「単位時間あたりの軌道の平均拡大縮小率」と定義するのが自然だと考えられます。そこで、初期値が非常に近い二つの軌道の間の距離を $\boldsymbol\delta(t)$ として、連続力学系の最大Lyapunov指数を次のようにあらわします。
$$ \lambda = \lim_{t \to \infty} \frac{1}{t}\ln\frac{|\boldsymbol\delta(t)|}{|\boldsymbol\delta(0)|} \tag{4}$$

この式を書き換えた次の式からも、最大Lyapunov指数が二つの軌道の距離が拡大するスピードとなっていることがわかります。
$$\frac{|\boldsymbol \delta(t)|}{|\boldsymbol \delta(0)|} \sim \exp(\lambda t) \tag{5}$$

連続力学系のLyapunov指数の数値計算

それでは3次元連続力学系のLorenz方程式のLyapunov指数を計算してみます。

1. ナイーブな方法（最大Lyapunov指数）

最も簡単な方法は、定義どおりに軌道の離れていくスピードを計算することです。ただし、この方法で求められるのは最大Lyapunov指数のみです。

以下のコードでは、初期値が$\boldsymbol \delta(0)$だけ離れた二つの軌道 $\boldsymbol x_1(t)$ と $ \boldsymbol x_2(t)$ の間の距離$\boldsymbol \delta(t)$を計算しています。

def f(t, X, sigma=10, beta=8/3, rho=28):
    x,y,z = X
    return np.array([sigma * (y - x), x * (rho - z) - y, x * y - beta * z])

def runge_kutta(func, t, X, dt):
    k1 = func(t, X)
    k2 = func(t + dt / 2., X + k1 * dt / 2.)
    k3 = func(t + dt / 2., X + k2 * dt / 2.)
    k4 = func(t + dt, X + k3 * dt)
    return X + (k1 + 2. * k2 + 2. * k3 + k4) * dt / 6.

T = 50
dt = 0.001
delta = np.array([0, 0, 1e-9])
x1 = np.array([2, 3, -14])
x2 = x1 + delta
ts = np.arange(0, T+dt, dt)
xs1 = np.empty((len(ts), 3))
xs2 = np.empty((len(ts), 3))
for i, t in enumerate(ts):
    xs1[i] = x1
    xs2[i] = x2
    x1 = runge_kutta(f, t, x1, dt)
    x2 = runge_kutta(f, t, x2, dt)
dist = np.linalg.norm(xs1-xs2, axis=1)

cutoff = int(25 / dt)
slope, intercept = np.polyfit(ts[:cutoff], np.log(dist[:cutoff]), 1)
dist_pred  = slope * ts + intercept
print(slope)

時間経過とともに$\delta(t)$がどのように変化したかを片対数グラフにプロットすると下図のようになります。青線は二つの軌道の間の距離が時間とともに離れていく様子をあらわしています。縦軸は対数になっており、この線の傾きが最大Lyapunov指数の近似値です。ここで、青線の傾きが$t=25$付近からフラットになっていることがわかります。これは、ある程度長い時間が経過すると二つの軌道の距離がLorenzアトラクタ全体のサイズと同程度になり、それ以上は離れられなくなるためです。したがって、直線のフィッティングには$t=25$未満の部分のみを使用しています。

plt.plot(ts, dist, label="distance")
plt.plot(ts, np.exp(dist_pred), label=f"fit (slope={slope:.3f})")
plt.vlines(cutoff*dt, 0, 1e10, ls="dotted", lw=1, color="gray")
plt.xlim(0, T)
plt.ylim(1e-10, 1e10)
plt.legend()
plt.yscale("log")
plt.xlabel("Time")
plt.ylabel("Distance")
plt.show()

結果は$\lambda=0.928$となりました。
この方法は初期条件やフィッティングに使う時間幅によって値が変わるため、安定した結果が得られるものではありません。

2. ヤコビ行列の利用

より安定した結果を得るために、ヤコビ行列を利用します。ヤコビ行列は軌道の各点における伸縮率を表します。それらを軌道に沿って平均することでLyapunov指数が計算できます。また、最大Lyapunov指数だけではなくLyapunovスペクトラムを計算することができます。

まず、変形行列というものを導入します。$\boldsymbol\delta(0)$から$\boldsymbol\delta(t)$への時間発展を式$(6)$のように書いたときに出てくる$\boldsymbol M(t)$が変形行列です。

$$\boldsymbol\delta(t) = \boldsymbol M(t) \boldsymbol\delta(0) \tag{6}$$

変形行列とその転置との積をとった行列 $\boldsymbol M^T \boldsymbol M$ の固有値を${\lbrace \alpha_i \rbrace}$とすると、Lyapunov指数$\lambda_i$は次のように表すことができます。
$$\lambda_i = \lim_{t \to \infty}\frac{1}{t}\ln \alpha_i^\frac{1}{2} \tag{7}$$

時間とともに変化する$M(t)$を計算するためには、$M(t)$の時間発展の方程式を書く必要があります。式$(6)$の両辺を時間$t$で微分すると、次のようになります。
$$\dot {\boldsymbol\delta} (t) = \dot {\boldsymbol{M}}(t) \boldsymbol \delta (0) \tag{8}$$

一方、ヤコビ行列$\boldsymbol{J}(\boldsymbol x(t))$を使うと、
$$\dot{\boldsymbol\delta}(t) = \boldsymbol J(\boldsymbol x(t)) \boldsymbol{\delta}(t) = \boldsymbol J(\boldsymbol x(t)) \boldsymbol{M}(t) \boldsymbol{\delta}(0) \tag{9}$$

なので、式$(7)$と式$(8)$を比較すると、

$$\dot {\boldsymbol{M}}(t) = \boldsymbol J(\boldsymbol x(t)) \boldsymbol{M}(t) \tag{10}$$

という微分方程式が得られます。

元の力学系の時間発展方程式と併せて、

\left\{
    \begin{align*}
        &\dot{\boldsymbol{x}}(t) = \boldsymbol{f}(\boldsymbol{x}(t)) \\
        &\dot{\boldsymbol{M}}(t) = \boldsymbol J(\boldsymbol{x(t)}) \boldsymbol{M}(t)
    \end{align*}
\right.

という方程式を一緒に解くことで、変形行列$\boldsymbol{M}(t)$を計算することができます。十分長い時間をとれば初期条件の違いによる差異や時間幅にほとんど依存しないLyapunovスペクトラムが得られます。

しかし、単純にこれを解こうとするとすぐに$\boldsymbol{M}$の要素が大きくなり過ぎて数値計算に困難が生じます。

3. QR分解による方法

数値計算におけるオーバーフローを回避しつつ$\boldsymbol{M(t)}$の固有値を計算するために、QR分解を使います。

まず、微分方程式$(10)$に現れる時間を離散化して、

$$\frac{ \boldsymbol{M}(t_{n}) - \boldsymbol{M}(t_{n-1})}{\delta t} = \boldsymbol J(\boldsymbol x(t_{n-1})) \boldsymbol M(t_{n-1}) \tag{11}$$

と書きます。すると、

\begin{align}
\boldsymbol{M}(t_{n}) 
&= \lbrack \boldsymbol I + \boldsymbol J(\boldsymbol x(t_{n-1}))\delta t \rbrack \boldsymbol{M}(t_{n-1})\\
&= \lbrack \boldsymbol I + \boldsymbol J(\boldsymbol x(t_{n-1}))\delta t \rbrack \lbrack \boldsymbol I + \boldsymbol J(\boldsymbol x(t_{n-2}))\delta t \rbrack \boldsymbol{M}(t_{n-2})\\
&= \lbrack \boldsymbol I + \boldsymbol J(\boldsymbol x(t_{n-1}))\delta t \rbrack \cdots \lbrack \boldsymbol I + \boldsymbol J(\boldsymbol x(t_{0}))\delta t \rbrack \boldsymbol{M}(t_{0}) \tag{12}
\end{align}

となり、$\boldsymbol M(t_n)$ を各時刻におけるヤコビ行列で表される行列$\boldsymbol M^{i} \equiv \boldsymbol I + \boldsymbol J(\boldsymbol x(t_i))\delta t$の積で

$$ \boldsymbol M(t_n) = \boldsymbol M^{n-1} \boldsymbol M^{n-2} \cdots \boldsymbol M^{0} \tag{13}$$

と表すことができます。

ここでQR分解を使って、

$$\boldsymbol M^{0} = \boldsymbol Q^{0} \boldsymbol R^{0}　\tag{14}$$

のように、$\boldsymbol M$を直交行列$\boldsymbol Q$と上三角行列$\boldsymbol R$の積に分解します。

さらに、$\boldsymbol M^{1} \boldsymbol Q^{0}$をQR分解して、

\begin{align}
\boldsymbol M^{1} \boldsymbol M^{0} &= 
\boldsymbol M^{1} \boldsymbol Q^{0} \boldsymbol R^{0}\\
&= \boldsymbol Q^{1} \boldsymbol R^{1} \boldsymbol R^{0} \tag{15}
\end{align}

とします。これを繰り返すと、

\begin{align}
\boldsymbol M(t_n) &=
\boldsymbol M^{n-1} \boldsymbol M^{n-2} \cdots \boldsymbol M^{0}\\
&= \boldsymbol Q^{n-1} \boldsymbol R^{n-1} \cdots \boldsymbol R^{1} \boldsymbol R^{0}\tag{16}
\end{align}

と書くことができます。

$\boldsymbol M^T \boldsymbol M$の固有値$\alpha_i$は、

\begin{align}
\alpha_i &= \hat{\boldsymbol \delta_i}^T \boldsymbol M(t_N)^T \boldsymbol M(t_N) \hat{\boldsymbol \delta_i}\\
&= \lvert \boldsymbol R \hat{\boldsymbol \delta_i} \rvert ^2\\
&= \lvert R_{ii} \rvert ^2 \tag{17}
\end{align}

となります。ここで$\boldsymbol R = \boldsymbol R^{n-1} \cdots \boldsymbol R^{1} \boldsymbol R^0$です。また、$\boldsymbol Q$は直交行列なので$\boldsymbol Q^T \boldsymbol Q= \boldsymbol I$となることを使っています。$\hat{\boldsymbol \delta_i}$は$i$番目の固有ベクトルと同じ方向の単位ベクトルです。

上述のLyapunovスペクトラムの式$(7)$より、

\begin{align}
\lambda_i &= \lim_{t \to \infty}\frac{1}{t}\ln \alpha_i^\frac{1}{2}\\
&= \lim_{N \to \infty}\frac{1}{N\delta t}\ln \lvert R_{ii}\rvert \\
&= \lim_{N \to \infty}\frac{1}{N\delta t} \sum_{n=0}^{N-1} \ln \lvert R_{ii}^n\rvert
\end{align} \tag{18}

となります。

このように各時刻の$\ln R_{ii}$を足し合わせていくことで、オーバーフローを回避しながらLyapunovスペクトラムを計算できます。

次のコードは実際にLorenz方程式のLyapunovスペクトラムをQR分解を利用して計算しています。

def jacobian(X):
    x,y,z = X
    return np.array([[-sigma, sigma, 0],[rho-z, -1, -x],[y, x, -beta]])

T = 100
dt = 0.001
x = [1, 5, -14]
ts = np.arange(0, T+dt, dt)
xs = np.empty((len(ts), 3))
lambda_ = np.empty((len(ts), 3))
Q = np.identity(3)
for i, t in enumerate(ts):
    J = jacobian(x)
    M = np.identity(3) + J * dt
    Q,R = np.linalg.qr(np.dot(M, Q))
    lambda_[i] = np.squeeze(np.asarray(np.log(np.abs(R.diagonal()))))
    x = runge_kutta(f, t, x, dt)
lambda_ = lambda_.cumsum(axis=0) / (np.expand_dims(ts, axis=1) + dt)

print(lambda_[-1])

この例では$T=100$で$\lambda_1=0.9315, \lambda_2=0.011, \lambda_3=-14.645$となりました。下図は時間を横軸にとってLyapunovスペクトラムをプロットしたものです。きちんと収束していることがわかります。

plt.plot(ts, lambda_[:, 0], label="$\lambda_1$")
plt.plot(ts, lambda_[:, 1], label="$\lambda_2$")
plt.plot(ts, lambda_[:, 2], label="$\lambda_3$")
plt.grid()
plt.xlabel("Time")
plt.ylabel("Lyapunov Exponents")
plt.legend()
plt.ylim(-20, 5)
plt.show()

参考文献

Chaos: An Introduction to Dynamical Systems (Kathleen T.Alligood, Tim D.Sauer, James A.Yorke)
Nonlinear Dynamics And Chaos (Steven Strogatz)
Govindan Rangarajan, Salman Habib, and Robert D. Ryne (1998). “Lyapunov Exponents without Rescaling and Reorthogonalization”. Physical Review Letters 80: 3747–3750.
Chalmers University of Technology Lecture Materials