ローレンツ方程式を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が低スペックでフリーズしたこともあり、今回はやめておきました。

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



スポンサードリンク

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

アッカーマン関数をpythonで書いたら瞬殺できた話。




スポンサードリンク

アッカーマン関数とは

先日、アッカーマン関数というものを初めて知ったので、それについて書いていきます。なんかのアニメのキャラを連想しますが笑

アッカーマン関数の定義を以下に示します。
\[
\mathrm{Ack}(m,n)=
\left\{
\begin{array}{ll}
n+1 , (m=0)\\
\mathrm{Ack}(m-1,1), (n=0)\\
\mathrm{Ack}(m-1,\mathrm{Ack}(m,n-1)) ,(m\neq n\neq 0)
\end{array}
\right.
\]

関数の中に、自身の関数が含まれますので再帰関数ということになります。計算自体はそんなに難しくありませんが、根気が必要です。というよりこのタイプの関数はプログラム書いた方が、早いですよね。

ついでに手計算するとこんな感じになります。

即興感満載で申し訳ないですが、だらだらと地道にやるだけです。難しい計算もないです。中学生、いや足し算と引き算しか使っていないので小学生でもできます。

ですが、pythonで書いたら一瞬でした。さすが関数型の言語というところでしょうか。

 

特徴

この関数の特徴は「解の発散性」にあります。試しにm=0〜3、n=0〜4の時の計算結果を3次元散布図でプロットして見ました。

少しわかりづらいので、違う角度から。

指数関数とも少し違う増加傾向がありますね。今、mを3までしか計算してませんが、m=4にすると、pythonが

”RuntimeError: maximum recursion depth exceeded ”とか言ってきます。数がでかすぎて計算できねぇよって。ウソーンと思ってちょっと調べると、

wikipedia参照

 

A(3,1)では13なのにA(4,1)では、65533になっています。nの増加率は2nとそんな変わらないですが、mの増加率はえぐいことになってます。この増加傾向がアッカーマン関数の特徴と言えるでしょう(適当笑)



スポンサードリンク

この関数何に使うの?

問題はそこです。こんなでかすぎる関数何に使うんだよって感じですが、グラフ理論でグラフの計量評価に使ったりするらしい(意味不)。他にも、単純な計算手続きというのを利用して、いろんな言語でこのプログラムを書いて、どの言語が一番早いのか試した人もいるそう(ついでにpythonはクソ遅いww)。

関数が定義されているくらいなので、使っている研究者は腐る程いるんでしょう。ついでに考案したのはドイツの数学者の「ヴィルヘルム・アッカーマン」という人。よくこんな関数思いつきますよね。

ついでに進撃のなんちゃらはドイツのどこかの町がモデルになってるみたいですが、何か関係あるかもですね。

最後に

2問だけ問題と回答を載せます。
\begin{eqnarray}
\mathrm{Ack}(1,1)
&=&\mathrm{Ack}(0,\mathrm{Ack}(1,0))\\
&=&\mathrm{Ack}(0,\mathrm{Ack}(0,1))\\
&=&\mathrm{Ack}(0,2)\\
&=&3
\end{eqnarray}
\begin{eqnarray}
\mathrm{Ack}(2,1)&=&\mathrm{Ack}(1,\mathrm{Ack}(2,0))\\
&=&\mathrm{Ack}(1,\mathrm{Ack}(1,1))\\
&=&\mathrm{Ack}(1,\mathrm{Ack}(0,\mathrm{Ack}(1,0)))\\
&=&\mathrm{Ack}(1,\mathrm{Ack}(0,\mathrm{Ack}(0,1)))\\
&=&\mathrm{Ack}(1,\mathrm{Ack}(0,2))\\
&=&\mathrm{Ack}(1,3)\\
&=&\mathrm{Ack}(0,\mathrm{Ack}(1,2))\\
&=&\mathrm{Ack}(0,\mathrm{Ack}(0,\mathrm{Ack}(1,1)))\\
&=&\mathrm{Ack}(0,\mathrm{Ack}(0,\mathrm{Ack}(0,2)))\\
&=&\mathrm{Ack}(0,\mathrm{Ack}(0,3))\\
&=&\mathrm{Ack}(0,4)\\
&=&5
\end{eqnarray}
なんだかゲシュタルト崩壊でも起こしそうですね。



スポンサードリンク



スポンサードリンク

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

Rossler方程式を Runge-Kutta法を用いて計算してみた(言語はpython)。



スポンサードリンク

 

多変量時系列解析の勉強その2
前回の記事で「アトラクター」という単語が出てきたので、今回はその用語を中心に書いていこうと思う。

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

とりあえずググると、アトラクターとは、「ある力学系がそこに向かって時間発展をする集合のことである。カオスな力学系に対してアトラクターを描写することは、現在においてもカオス理論における一つの研究課題である。」とか「実世界における力学系は散逸的であることが多いであろう。すなわち、もし力学系に運動の駆動力が無ければ、運動は停止するものと考えられる(そのような散逸は、様々な原因による内部摩擦熱力学的損失、物質の損失などにより生じうる)。散逸と駆動力が組み合わさることにより、初期の摂動を鎮め、その力学系の振る舞いを典型的なものへと落ち着かせる傾向にある。そのような典型的な振る舞いに対応している力学系からなる位相空間の一部分はattracting section または attractee と呼ばれる。」とかウィキペディアには書いてある。

しかしこれだけではもちろん意味不なので、論文を読んでみた。

>>>読んだ論文はこちら

今回はこちらの内容を元に、勉強したことを書いていこうと思います。

埋め込みとアトラクタ表示

例えばある3変数関数f(x,y,z)があるとする。しかし得られたデータは状態変数xの時系列データ(x(t))であるとする。そしてここから遅れ時間\(/delta\)なる変数を新たに導入し、新たにn次元ベクトル\({\bf V}= (x(t),x(t+\delta),…,x(t+(n-1)\delta))(t = 0, 1, · · · , N − 1 − (n − 1)\delta)\)を生成する。

以下では、例として、連続力学系のRossler方程式を取り上げる。
\[
\left\{
\begin{array}{ll}
\frac{dx}{dt}= -y-z\\
\frac{dy}{dt}= x+ay\\
\frac{dz}{dt}= b+z(x-c)
\end{array}
\right.
\]
なお今回は、a = 0.398, b = 2, c = 4 とし,ルンゲクッタ法を用 いて数値計算を行なった.

まず、全ての状態変数(x、y、z)を用いて軌道を描いたアトラクターが以下である。

図1 全ての状態変数を用いて描いたアトラクター

次に、3つの状態変数のうちxのみが観測可能であると仮定し、遅れ時間\(\delta\)を用いて。アトラクタを再構成したものが以下である。なお遅れ時間\(\delta=1,10,500\)を用いて時系列データを作成した。

図2 遅れ時間[latex]\delta=1[\mathrm{s}][/latex]の再構成アトラクタ
図3 遅れ時間[latex]\delta=10[\mathrm{s}][/latex]の再構成アトラクタ
図4 遅れ時間[latex]\delta=500[\mathrm{s}][/latex]の再構成アトラクタ
 

図からわかる通り、図3のアトラクタの位相は全ての状態変数が観測できた時に描けるアトラクタ(図1)と位相構造が似ている。

そして今回の計算ではTakensの埋め込み定理は満たされていないが、それでも元の力学系の位相的構造が保存されていることがわかる。

本来、アトラクタの再構成を適切におこなうには、時間遅れ\(\delta\)は時系列データの「自己相関関数」や「パワースペクトラム」を用いるらしいが、今回のような状態変数が少なく、データを視覚的に捉えやすい場合は、必ずしもそうする必要がないことがわかる。



スポンサードリンク

Runge-Kutta法によるR\(\ddot{\mathrm{o}}\)ssler方程式の計算

以下に、今回の計算に用いたプログラムを載せる。

なお、上のプログラムは全ての状態変数が観測できた時の、アトラクタを描画する際のプログラムである。

 

本来なら、このプログラムに時間遅れを加味した時系列データをいい感じに元のデータにappendでもしたり、二次元配列のデータの任意の一部のデータを除去したりできればよかったです。ですがそこまでの実力がなかったので(一応数時間奮闘はしました笑)、時間遅れの再構成のアトラクタは別にプログラムを書いた(とは言ってもデータ処理の部分のみ)。

 

xdata.csvは観測できる時系列データです。そして、xdata-2とxdata-3は時間遅れの新規時系列データx(t+\(\delta\))とx(t+2\(\delta\))ということになります。

これらのデータはとても原始的に用意しました笑。xdataの時系列データの時間遅れの分だけ上の行から取り除いたデータを、xdata2にコピペし、さらに、そのxdata2の時系列データから、時間遅れの分だけ上の行から取り除いたデータを、xdata3にコピペする。

最後に時間遅れの分、xdata2,xdata3は、要素数(ファイルの行数)が少ないので、一番データ行数の少ない、xdata3と行数が合うようにデータを消して、3つのファイルの行数を揃えます。

こんな原始的にやらなくても、pythonならこれくらいのこと瞬殺できそうですが、ここはまた今度にします(pythonの無駄遣い笑)。

 

最後に

今回は「多変量時系列解析」に関してしっかりと理解するため、論文のような計算結果が得られるか確認しました。ルンゲクッタとかはググったら瞬殺できたのでよかったのですが、pythonに慣れていなく、データ処理の所に時間がかかってしまいました。

今後はpythonのデータ処理と、ローレンツ方程式ではどんな計算結果が得られるのか勉強したいと思います。

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