背景
X(t+1) = F(X(t))のような離散時間力学系に対しては、QR分解の方法とかがいろんな所で書かれているけど、 \dot{X} = F(X, t)のような連続時間力学系に対しては、あまりリアプノフ指数を計算する方法が書かれていないので、メモしておく。
リアプノフ指数
\dot{X} = F(X, t)に対して\dot{M} = DF(X(t), t) Mで定義される、Mの特異値 \alpha_i を用いて、リアプノフ指数 \lambda_i は以下で定義される。
\lambda_i = \lim_{T\rightarrow \infty} \frac{1}{T} \log(\alpha_i)
最大リアプノフ指数だけ計算する方法
適当なvを用意して、 v(t)の拡大率を計算すればいい。
で計算できる。少しでもリアプノフベクトル方向の成分を持っていれば、十分時間が経った時に、その方向だけにずば抜けて指数的に拡大するので、そこからわかる。
上位k個ののリアプノフ指数だけ計算する方法
島田,長島(1979)の方法
k本のランダムベクトル d^{(1)}_0(0),\dots,d^{(2)}_0(0) を用意する。ただし、d^{(i)}_0(0)^T d^{(j)}_0(0) = \delta_{ij}である。
\tau秒時間発展させて、 d^{(1)}_0(\tau),\dots,d^{(2)}_0(\tau) を得る。ここで、ノルム降順に並び替える。
そこで、Gram-Schmitの方法で、直行化してを d^{(1)}_{0\perp}(\tau), \dots, d^{(k)}_{0\perp}(\tau)を得る。
これをつかって、d_1^{(i)}(0) = \frac{|d_0^{(i)}(0)|}{|d_{0\perp}^{i}(\tau)|}d_{0\perp}^{i}(\tau)を計算して、ワンステップが終わる。
\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=QRとQR分解をする。この時、QとRをそれぞれ別々に時間発展することを考える。
Q^T\dot{Q} - Q^T DF(X(t), t) Q = - \dot{R} R^{-1}
この時間発展方程式では、右辺は上三角行列であり、Q^T\dot{Q}は反対称行列になっている(Q^TQ = Iを微分すればわかる)。そこからS:=Q^T\dot{Q}としたときに、Sの下三角成分はQ^T DF(X(t), t) Q と一致していて、反対称なので上三角成分もこれで決まる。よって、S = triu(Q^T DF(X(t), t) Q) - triu(Q^T DF(X(t), t) Q)^Tなので、
からQを時間発展できる。
また、右辺の対角成分は \dot{R_{ii}}/R_{ii} = \dot{\log(R_{ii})}なので、
d(\log(R_{ii}))/dt = diag(Q^T DF(X(t), t) Q)
で \log(R_{ii})を時間発展させることができる。R_{ii}はMの特異値に一致するため、これから、\lambda_i = \lim_{T\rightarrow \infty}\frac{1}{T}\log(R_{ii})でリアプノフ指数を計算できる。
まとめると、
- 目的の力学系を数値積分して、X(t)を解く
-
Qと\log(R_{ii})を上の方程式に従って数値積分する
-
\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
Discussion