ローレンツ方程式をRunge-Kutta法(python)で計算



スポンサードリンク

前回は、アトラクターの勉強ということで、非線形力学系のレスラー方程式のアトラクター再構成の数値実験を行いました。

>>>前回の記事「Rossler方程式を Runge-Kutta法を用いて計算してみた(言語はpython)。」

>>>前々回の記事「多変量時系列解析に関する論文」

今回はローレンツ方程式について実施した数値実験の結果と、時間遅れの大きさについて論じて行きます。

時間遅れ座標系を用いた、アトラクタの再構成

まず、初めにローレンツ方程式を示します。
\[
\left\{
\begin{array}{ll}
\frac{dx}{dt}= -px+py\\
\frac{dy}{dt}= -xz+rx-y\\
\frac{dz}{dt}= xy-bz
\end{array}
\right.
\]
x,y,z 全ての変数の時系列データが観測できた場合のアトラクタを以下に示します。

図1 (x,y,z)の時系列データから得られるアトラクター

アトラクタ再構成のやり方には

  • 微分座標系(V(t)=(\(x(t),\dot x(t),\ddot y(t)\))を用いる
  • 時間遅れ座標系(V(t)=(\(x(t),x(t+\delta),x(t+2\delta)\)))を用いる

方法があるが、ノイズがある場合には、微分座標系はノイズを拡大させる原因にもなるので、時間遅れ座標系を用いる場合が多いです。

以下に観測できるデータがx(t)、時間遅れの大きさ\(\delta=1,10,100,500\)の場合での再構成したアトラクタを示します。

 

図2 時間遅れの大きさδ=1
図3 時間遅れの大きさδ=10
図4 時間遅れの大きさδ=100

 

 

図5 時間遅れの大きさδ=500

 

 

 

 

 

 

 

 

 

 

そして力学系のアトラクタは平衡点、リミットサイクル、kトーラス、ストレンジアトラクタ、ランダムに分類されます。

状態空間(位相構造)を観察する限りでは、ローレンツ方程式のアトラクタはストレンジアトラクタというものに分類されます。

このアトラクターが「ストレンジアトラクターになっている」というのが大切で、これが観察できれば、この非線形力学系は「カオス的な振る舞いをする」と言えます。

ターケンスの埋め込み定理

しかし、上記のような時間遅れ座標系を用いてアトラクタを再構成するには、どうやら以下のような定理が成り立たないといけないみたいです。

「d 次元のコンパクトな多様体 M と C 1 級微分同相写像 f : M → M,C 1 級関数 g : M → R 1 が与えられたとき, m = 2 d + 1 であ れば,次式の写像 V : M → R m は,埋め込みである」

だめだ、わけワカンねぇ。本当に同じ言語なのかと疑うほどであるが、m=2d+1が成り立てば良いらしい。

ここでmとはアトラクタの次元のことで、ここではm=3になります。

問題はdです。この「ターケンスの埋め込み定理」が成り立つには、d=1であれば良いのですが、僕はこのdが何を意味するのかわかりません。可観測な状態変数の数であるならば、今回は可観測なデータがx(t) のみを仮定しているので、成り立ちますがどうなのでしょう。

ここら辺、追求したいですが、体力が切れてきたので、ここはまた今度にします。

しかし前回読んだ論文では、ターケンスの埋め込み定理はあくまで十分条件らしいので、必ず成り立っていないといけないみたいではないようです。


スポンサードリンク

時間遅れの大きさの最適化

時間遅れの大きさの最適化には以下の二つの方法があります。

  • 線形自己相関関数が最初に0になる点を求める
  • 平均相互情報量が最初に0なる点を求める

以下では平均相互情報量の方法について述べて行きます。まずは定義から

\begin{eqnarray}
\hat{I}_{AB}(T)=\sum_{a_{i},b_{k}} P_{AB}(a_{i},b_{k})I_{AB}(a_{i},b_{k})
\end{eqnarray}
A={x(n)},B={x(n+T)}とすれば、
\begin{eqnarray}
\hat{I(T)}=\sum_{n=1}^N P(x(n),x(n+T))\log2(\frac{P(x(n),x(n+T))}{P(x(n))P(x(n+T))})
\end{eqnarray}
となります。そしてこの\(\hat{I(T)}\) の最初(シグマなので)の最小値をとる時のTを「最適な遅れ時間」として選択します。

そして、この遅れ時間の最適化の計算はプログラムを書いて計算し、また次回以降載せます。

計算するとδ=18くらいになるみたいです。上の再構成アトラクタ(図2〜図5)もδ=10が一番元の位相構造に近いので、適切な手法と言えそうです。

計算に用いたプログラム

以下に計算に用いたプログラムを載せます。

今回は、遅れ時間の時系列データもプログラムできました。しかし、遅れ時間のデータをfor文で別々に書いているので、ここを一つにまとめられたらよかったのです。

それと、遅れ時間も毎回、δの部分を書き換えて計算しているので、δ の配列を作り、配列の要素ごとに計算するというのも、やってみたいです。

アトラクタの軌跡のアニメーションも作ろうと思って計算したのですが、PCが低スペックでフリーズしたこともあり、今回はやめておきました。

そこらへんも、次回以降書いて行きたいと思います。



スポンサードリンク

記事が役に立ったらシェア!