背景
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)の拡大率を計算すればいい。
で計算できる。少しでもリアプノフベクトル方向の成分を持っていれば、十分時間が経った時に、その方向だけにずば抜けて指数的に拡大するので、そこからわかる。
リアプノフ指数を全部計算する方法
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