ローレンツアトラクタ

ローレンツアトラクタはカオスの中で最もポピュラーなものでしょう。蝶の羽のような図はWeb上にもたくさんあります。Excelを用いたもの、Flashを用いたもの、POV-Rayを用いたもの、JAVAを用いたものなど、アトラクタを描くためのソフトウェアもさまざま。Mathematicaにはローレンツアトラクタを描く例題がありますし、Mac OSXに付属のGrapherにも例があります。

Mathematicaの例

Grapherの例

また、「chaoscope」というフリーソフトでは条件(係数や初期値を)入力するだけで、ローレンツアトラクタに限らずいろいろなアトラクタが描けます。しかし、図を鑑賞するだけならそれでよいでしょうが、ここはやはりプログラムを作ってシミュレーションをしてみたいところです。そうすれば、ローレンツアトラクタと似たレスラーアトラクタもすぐに描けます。

この節では、常微分方程式であるローレンツの式から出発し、どのようにローレンツアトラクタを描くのか、その意味は何かを考えていきます。

ローレンツの式

ローレンツアトラクタとは、気象学者のローレンツが大気変動モデルを研究していて導き出した微分方程式の解が描く図形です。微分方程式を解くときは、初期値が必要になりますが、その初期値がわずかに変わっただけで結果が大きく変動し、まさに予測がつかないというわけです。地球の反対側のチリの海岸で誰かがたき火をすると小さな上昇気流が生まれ、日本の台風の進路に影響する、そんなたとえ話ができます。

ローレンツの式は次の連立常微分方程式です。(常微分方程式という言葉は高校数学の範囲では聞きなれない言葉でしょう。高校数学でいう微分方程式と同じものと思ってください)

ここで、xは対流運動の速度、y,zは水平方向と垂直方向の温度変位、係数sはプラントル数、rはレイリー数、bは系の物理的サイズに関わる数です。したがって、空間座標のx,y,z ではないのですが、これを空間座標と見なして図を描くのです。このようなものを「相空間」といいます。

計算機で扱うために、これを差分方程式に直します。(差分方程式については「微分方程式を解く」を参照してください)

    

「Excelで学ぶ理工系シミュレーション入門」(臼田、伊藤、井上共著 CQ出版社)ではExcelを使って4000回計算をした表を作り、散布図として図示しています。しかし、Excelでは3次元の図は描けないので、xyとxzの 2通りを表示しています。CindyScriptなら、3次元の描画も可能なので、3次元で表示しましょう。そのために、まず3D基本パッケージをダウンロードしてください。(ナビゲータの基本パッケージの中にあります)

3つの差分方程式をそれぞれ、関数として定義します。この中の係数は定数として用意しておきます。ローレンツが特に注目したのは、s=10 , b=8/3 , r=28 です。Δtも小さな値を割り当てておきます。

初期値を決めて、繰り返し処理で次々に(x,y,z)の組を求め、draw()で点を打っていきます。

次のスクリプトを、3D基本パッケージのDrawスロットに書きます。(3D基本パッケージではDrawスロットに既にコードが書かれていますので、そのあとに書きます)行の最初は説明のためにつけた行番号です。実際には不要なので、カットアンドペーストした場合は削除してください。

-------------------------------------------------------------------------------------------------------------------------

1: s=10;

2: b=8/3;

3: r=28;

4: dt=0.004;

5: fx(x,y):=x-dt*s*(x-y);

6: fy(x,y,z):=y+dt*(-x*z+r*x-y);

7: fz(x,y,z):=z+dt*(x*y-b*z);

8: work=[10,12,25];

9: repeat(4000,

10:   x=work_1;

11:   y=work_2;

12:   z=work_3;

13:   work=[fx(x,y),fy(x,y,z),fz(x,y,z)];

14:   draw(map(work),size->1,color->[1,0,0],noborder->true);

15:   );

--------------------------------------------------------------------------------------------------------------------------

1行目から4行目は定数の設定。5行目から7行目が差分化した微分方程式を関数としたものです。

8行目、最初に初期値として [10,12,15] を設定しました。

9行目 work は差分化した微分方程式から計算するときの作業用の変数。アンダーバーは、リストの要素を取り出すための演算子です。work=plist_1 は、リストplistの1番目の要素すなわち [12,12,25]を表します。それをworkとします。

9行目から15行目までが繰り返し処理。4000回繰り返します。

10,11,12行目はx,y,zにworkのそれぞれの要素を代入しています。最初はx=10,y=12,z=25 となります。

13行目で、fx() , fy() , fz() で計算した結果を作業用変数workに代入しています。(x1,y1,z1)から(x2,y2,z2)が求められたことになります。

14行目で3次元空間に点を打ちます。map()という関数は3D基本パッケージに入っている関数です。

color->[1,0,0]は表示色の指定。[1,0,0]は赤。RGBの順に指定しますので、[0,0,1]は青

noborder->true は点の境界線(周囲)を書かないという指定。このオプションをつけなければ、点の境界線が黒で表示されます。

 

この節はまだ書きかけです