[数値計算] 連続時間力学系のリアプノフ指数の計算まとめ

背景

$X(t+1)=F(X(t))X(t+1) = F(X(t))$のような離散時間力学系に対しては、QR分解の方法とかがいろんな所で書かれているけど、 $X˙=F(X,t)\dot{X} = F(X, t)$のような連続時間力学系に対しては、あまりリアプノフ指数を計算する方法が書かれていないので、メモしておく。

リアプノフ指数

$X˙=F(X,t)\dot{X} = F(X, t)$に対して$M˙=DF(X(t),t)M\dot{M} = DF(X(t), t) M$で定義される、$MM$の特異値 $αi\alpha_i$ を用いて、リアプノフ指数 $λi\lambda_i$ は以下で定義される。

$$λi=lim⁡T→∞1Tlog⁡(αi) \lambda_i = \lim_{T\rightarrow \infty} \frac{1}{T} \log(\alpha_i)$$

最大リアプノフ指数だけ計算する方法

適当な$vv$を用意して、 $v(t)v(t)$の拡大率を計算すればいい。

$$v˙=DF(X(t),t)v \dot{v} = DF(X(t), t) v$$

で計算できる。少しでもリアプノフベクトル方向の成分を持っていれば、十分時間が経った時に、その方向だけにずば抜けて指数的に拡大するので、そこからわかる。

上位k個ののリアプノフ指数だけ計算する方法

島田,長島(1979)の方法

k本のランダムベクトル $d0(1)(0),…,d0(2)(0)d^{(1)}_0(0),\dots,d^{(2)}_0(0)$ を用意する。ただし、$d0(i)(0)Td0(j)(0)=δijd^{(i)}_0(0)^T d^{(j)}_0(0) = \delta_{ij}$である。
$τ\tau$秒時間発展させて、 $d0(1)(τ),…,d0(2)(τ)d^{(1)}_0(\tau),\dots,d^{(2)}_0(\tau)$ を得る。ここで、ノルム降順に並び替える。
そこで、Gram-Schmitの方法で、直行化してを $d0⊥(1)(τ),…,d0⊥(k)(τ)d^{(1)}_{0\perp}(\tau), \dots, d^{(k)}_{0\perp}(\tau)$を得る。
これをつかって、$d1(i)(0)=∣d0(i)(0)∣∣d0⊥i(τ)∣d0⊥i(τ)d_1^{(i)}(0) = \frac{|d_0^{(i)}(0)|}{|d_{0\perp}^{i}(\tau)|}d_{0\perp}^{i}(\tau)$を計算して、ワンステップが終わる。

$$∑mλm=lim⁡n→∞1n∑kn−1ln⁡dk(1)(τ),…,dk(m)(τ)が作る超平行m面体の体積dk(1)(0),…,dk(m)(0)が作る超平行m面体の体積 \sum_m \lambda_m = \lim_{n\rightarrow\infty} \frac{1}{n} \sum_k^{n-1} \ln \frac{d_k^{(1)}(\tau),\dots,d_k^{(m)}(\tau)が作る超平行m面体の体積}{d_k^{(1)}(0),\dots,d_k^{(m)}(0)が作る超平行m面体の体積}$$

リアプノフ指数を全部計算する方法

QR分解を使った方法を挙げておく。計算の中で指数関数的に拡大する項がないので、計算が安定する。
$M=QRM=QR$とQR分解をする。この時、$QQ$と$RR$をそれぞれ別々に時間発展することを考える。

$$QTQ˙−QTDF(X(t),t)Q=−R˙R−1 Q^T\dot{Q} - Q^T DF(X(t), t) Q = - \dot{R} R^{-1}$$

この時間発展方程式では、右辺は上三角行列であり、$QTQ˙Q^T\dot{Q}$は反対称行列になっている($QTQ=IQ^TQ = I$を微分すればわかる)。そこから$S:=QTQ˙S:=Q^T\dot{Q}$としたときに、$SS$の下三角成分は$QTDF(X(t),t)QQ^T DF(X(t), t) Q$ と一致していて、反対称なので上三角成分もこれで決まる。よって、$S=triu(QTDF(X(t),t)Q)−triu(QTDF(X(t),t)Q)TS = triu(Q^T DF(X(t), t) Q) - triu(Q^T DF(X(t), t) Q)^T$なので、

$$Q˙=QS \dot{Q} = QS$$

から$QQ$を時間発展できる。

また、右辺の対角成分は $Rii˙/Rii=log⁡(Rii)˙\dot{R_{ii}}/R_{ii} = \dot{\log(R_{ii})}$なので、

$$d(log⁡(Rii))/dt=diag(QTDF(X(t),t)Q) d(\log(R_{ii}))/dt = diag(Q^T DF(X(t), t) Q)$$

で $log⁡(Rii)\log(R_{ii})$を時間発展させることができる。$RiiR_{ii}$は$MM$の特異値に一致するため、これから、$λi=lim⁡T→∞1Tlog⁡(Rii)\lambda_i = \lim_{T\rightarrow \infty}\frac{1}{T}\log(R_{ii})$でリアプノフ指数を計算できる。

まとめると、

1. 目的の力学系を数値積分して、$X(t)X(t)$を解く
2. $QQ$と$log⁡(Rii)\log(R_{ii})$を上の方程式に従って数値積分する
3. $λi=1Tlog⁡(Rii(T))\lambda_i = \frac{1}{T} \log(R_{ii}(T))$で十分大きな＄T＄を取ればリアプノフ指数に収束させることができる。

参考

[1] Lyapunov Exponents for Continuous-Time Dynamical Systems; Janaki T M, et al; India; 2003