この文章は、力学系研究のためのプログラムパッケージ GAIO を使う上で私の気づいた事をつれづれなるままに書き連ねたものです。 あくまで個人的メモですので、内容については完全に無保証です。悪しからず。
GAIO とは Global Analyis of Invarian Objects の頭文字であり、 Michael Dellnitz、 Oliver Junge を中心とする力学系研究者のグループにより開発された、力学系研究のためのプログラム群の総称です。
GAIO は subdivision argorithm と彼らが呼ぶ方法によって、アトラクターや不変多様体の covering を見つけだしてくれます。 また計算結果を Matlab のグラフィックス機能を用いて図示することができるので、lorenz アトラクターをマウスでぐりぐり回転させたり、アニメーションをさせるといったことが簡単に出来ます。
現在のバージョンの GAIO は Matlab への追加パッケージ (Toolbox) という形で実装されており、GAIO を使うためには Matlab が必要になります。 GAIO の開発は Linux 上で行なわれているようですが、Matlab が動作するたいがいの Unix クローン上で使えると思います。 また、一部の機能を用いるためには有料のライブラリ LEDA が必要になりますが、不変測度関係の計算をしない限りは LEDA 抜きでも十分使えます。
GAIO を入手するためには開発者にメールで連絡してソースファイルをもらう必要があります。 ウェブページ からたどってみてください。
論文を読め。
ディレクトリ "GAIODIR" 以下に GAIO のソースアーカイブを展開してあるとします。 GAIO の実体は GAIODIR 以下の各ディレクトリにある *.m や *.c という名前のファイルたちです。
.m という拡張子はそのファイルが Matlab のスクリプトファイルであることを表し、 拡張子を除いたファイル名を Matlab に打ち込むとこのファイルが呼ばれます。 要するに Matlab のコマンドが単に並んでいるテキストファイルですが、先頭に簡単な解説がついているので一読をおすすめします。
.c という拡張子を持つファイルは C言語のソースファイルです。 主にスクリプトでは速度的に問題がある関数が C言語で書かれています。 これらのファイルを Matlab から呼び出せるようにするには、次節で解説するように、コンパイルして shared object に変換しなくてはなりません。
README を読んでがんばれ。 gcc にコンパイルオプション "-lc" を渡さないと shared object の作成でコケるかもしれない。
とりあえず GAIODIR/doc/matlab/dansebook.ps の11ページにある lorenz アトラクターの例を実行してみましょう。環境変数 GAIODIR が正しくセットされていることを確認して matlab を起動します。 ユーザのホーム以下に GAIO を置いている場合は addpath /home/ユーザ名/gaio/gaio と入力して GAIODIR を matlab の検索パスに加えてください。 ここで次のコマンドたちを順次実行します。 先頭の数字は以下の説明のためにつけた行番号です。
| 01: | lorenz = Model('lorenz') |
| 02: | rk4 = Integrator('RungeKutta4') |
| 03: | rk4.model = lorenz |
| 04: | rk4.tFinal = -0.1 |
| 05: | rk4.h = -0.01 |
| 06: | tree = Tree([0,0,0], [120,120,160]) |
| 07: | tree.integrator = rk4 |
| 08: | tree.domain_points = Points('Edges', 3, 100) |
| 09: | tree.image_points = Points('Center',3) |
| 10: | x = [0;0;0] |
| 11: | depth = 21 |
| 12: | tree.insert(x,depth) |
| 13: | steps = 10 |
| 14: | gum(tree,depth,steps) |
エラーがなければ step 1: 40 boxes ... step 10: 19351 boxes とかなんとか出力されるでしょう。これで原点の安定多様体を含む19531個のボックスが得られました。 plotb(tree) と入力すれば matlab の plot 機能を呼び出してこれらのボックスを図示することができます。
上の例において各行が何をしているのか簡単に見てみます。 実際に計算を実行しているのは14行目だけで、それまでは全て状況の設定です。 各コマンドについてより詳しくは後の節を参照してください。
まず1行目では、変数 lorenz によって lorenz 方程式のモデルが呼び出せるように設定します。 右辺で GAIODIR/gaio/@Model/Model.m が呼び出され、引数が 'lorenz' なので Model.m は GAIODIR/models/アーキテクチャ名/lorenz.so というファイルを変数 lorenz に渡します。 ここで左辺の変数 lorenz と右辺のモデル名 lorenz を混同しないようにしてください。 左辺は変数なので本当は何でも良いのです。 例えば1行目を hogehoge = Model('lorenz') にした場合は、3行目も rk4.model = hogehoge にすれば同じ結果が得られます。
2行目から5行目では常微分方程式の数値積分をどのように行なうかを設定します。 2行目では変数 rk4 によって4次のルンゲ・クッタ法が呼び出されるように設定します。 GAIODIR/integrators に他の積分メソッドもありますし、自分で追加することもできます。 3行目では rk4 に先に定義した lorenz を積分するように指示します。 4行目は積分を止める時間、5行目は差分するときの1ステップの長さの定義です。 これらが負の数値になっているのは、ここで求めようとしているのが安定多様体だからであり、不安定多様体を見たい場合は正負を反転させます。
6行目から9行目において、"tree" の設定を行ないます。 まずは6行目において一番大元になる box を、中心が原点で各座標軸方向の大きさがそれぞれ 120、120、160 となるような箱として定めます。 7行目では tree を力学系によって分割していくときに、rk4 を呼び出すように設定します。 rk4 は先に設定した数値積分法であり、rk4.model = lorenz なので、変数 tree に対して、それに作用する力学系として lorenz が関連づけられたことになります。 8行目と9行目は、力学系によって変形された box がもとの box たちとどのように交わるかを判定するための方法を決定しています。
10行目と11行目は単に12行目で用いる数値を変数に代入しているだけです。 数値を代入すると 12行目は tree.insert([0;0;0],21) となりますが、これは6行目で設定した box を21回分割して、その中から原点を含む box を tree に追加しろという命令です。 box の分割は各座標方向に順に行なわれ、今の場合だと相空間は3次元なのでそれぞれの座標方向に7回分割が行なわれます。 すなわち、原点を中心とする大きな box が各座標方向に 2^7 = 128 分割されて出来た 2097152 個の小さな box たちの中から、原点を含む box のみに注目して、tree に取り込んだわけです。 便宜上この box を b_0 と名付けましょう。
14行目で計算が行なわれます。 現在 tree の 21回目の分割において選ばれているのは b_0 のみなので、b_0 を rk4 で送ります。 すなわち、ルンゲ・クッタ法で -0.1秒の数値積分を行ないます。 この積分により b_0 は変形され、もとの b_0 からはみ出して周囲の box たちと交わりを持ちます。 どの box と交わるかを8行目と9行目の設定に基づいて判定し、交わると判定された box を tree に追加します。 これで1ステップが終わりです。 次のステップでは、一回目のステップで選ばれた box たちに対して同じように数値積分を行ない、交わる box の判定をします。 これを steps = 10回くり返して最終的に出来上がった tree の box の数が 19351 ということになります。
作りたいモデルに似ているモデルを GAIODIR/models から見つけて、少しずつ変形していくとよいでしょう。 多少 C言語の知識が必要になります。 モデル名.c という名前でファイルを保存して、GAIODIR/models で make を実行すると、GIODIR/models/アーキテクチャ名/ というディレクトリに モデル名.so というファイルができます。 これで作ったモデルを dansebook.ps にある既製のモデルと同様に Model('モデル名') で呼び出せます。