第20回カオス

フラクタルと並んでとりあげられるのがカオス理論。「カオス」は日本語で「渾沌」といいますが、数学での「カオス」は、ちょっとニュアンスが異なります。その代表格は、ロジスティック写像とストレンジアトラクタ。この回は、ストレンジアトラクタのうち、有名なローレンツアトラクタをCinderella.2で表示します。

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

その微分方程式は次の３つ。（係数の文字は本やWebページによって少しずつ異なります）

計算機で扱うために、これを「差分」という形に直します。

微分の定義を思い出してください。

ここで、Δt が十分に小さければ、最初の式は

と見なせます。すると、初期値のx1,y1からx2が求められます。これを漸化式のようにして、x2,y2からx3、x3,y3からx4・・・と続けて求めていき、得られたx,y,zを座標としてプロットすると図が描けるというわけです。

「Excelで学ぶ理工系シミュレーション入門」（臼田、伊藤、井上共著ＣＱ出版社）ではExcelを使って4000回計算をした表を作り、散布図として図示しています。

ここで一つ問題があります。

ローレンツの式は、空間座標(x,y,z)の式です。しかし、Excelの散布図では３Ｄプロットはできないので、(x,y)と(x,z)に分けて２つ図示しています。Cinderella.2は、もとが平面幾何のソフトですので、そのままでは立体は表示できません。しかし、CindyScriptで行列を使って３Ｄグラフィクスを実現することができます。

まずは上記の本と同様、xyとxzの２つを表示することを考えましょう。３Ｄにするのはそのあとの課題とします。

まず、差分化した３つの方程式を眺めます。

これらは、関数として定義すればよいでしょう。この中の係数は定数として用意しておきます。ローレンツが特に注目したのは、s=10 , b=8/3 , r=28 です。Δtも小さな値を割り当てておきます。

初期値を決めて、繰り返し処理で次々に(x,y,z)の組を求めていきますが、その結果はリストにすればよいでしょう。

リストについては「CindyScript入門」のリストの節を参照してください。

（行の最初は説明のためにつけた行番号です。実際には不要なので、カットアンドペーストした場合は削除すること）

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

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: plist=[[10,12,25]]; //初期値

9: work=plist_1;

10: repeat(4000,

11: x=work_1;

12: y=work_2;

13: z=work_3;

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

15: plist=append(plist,work);

16: )

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

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

8行目の plist が計算結果のリスト。[x,y,z]で空間座標を表すことになります。その[x,y,z]のリストで、最初に初期値として [10,12,15] を設定しました。

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

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

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

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

15行目でその結果をリストに追加します。append()はリストに追加する関数です。

では、ここまでで動作を確かめてみましょう。4000では多すぎるので、3くらいにしましょう。計算した結果を最後に表示します。

空間座標 (x,y,z) からなるリストができています。

つぎに、これをプロット・・・したいのですが、ひとまずはxy,xzの２つに分解して表示しましょう。計算の都度表示してしまえばよいのでrepeatの繰り返しの中に入れてしまいます。また、プロットだけであれば、値を記録していく必要もないのでplistの作成も不要です。repeatの中身を次のようにします。

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

10: repeat(4000,

11: x=work_1;

12: y=work_2;

13: z=work_3;

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

15: listxy=[fx(x,y),fy(x,y,z)];

16: listxz=[fx(x,y),fz(x,y,z)];

17: draw(listxy,size->1,color->[1,0,0],noborder->true);

18: draw(listxz,size->1,color->[0,0,1],noborder->true);

);

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

15,16行目で(x,y) , (x,z) のペアを作り、17,18行目で表示します。これをまとめて

draw([fx(x,y),fy(x,y,z)],size->1,color->[1,0,0],noborder->true);

draw([fx(x,y),fz(x,y,z)],size->1,color->[0,0,1],noborder->true);

としてしますことも可能ですが、どちらが読みやすいかはなんともいえないところです。

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

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

歯車アイコンをクリックして実行すると、次のような図が描かれます。

定数や初期値を変えると図が変わります。

次に、これを３Ｄで表示します。

まずは、３ＤプロットができるMathematicaではどうなるのかやってみました。

MathematicaとCindyScriptでは、括弧の扱いなどが違いますが、考え方は同じです。

まずは、CindyScriptと同様、(x,y)と(x,z)に分けて表示してみます。

ここでは、Mathematicaの ListPlot関数を使ってみました。つまり、計算の都度表示するのではなく、計算結果はリストにため込んでいって、最後にまとめてプロットします。

つぎに、空間座標のリストplistをListPointPlot3D[] という関数で表示します。

これが３次元のローレンツアトラクタです。Mathematicaでは、マウスをドラッグして視点を変えることができます。もう少し下にすると

次に、Cinderella.2で３Ｄ表示をします。そのためには、３Ｄ表示をするためのスクリプトが必要になります。第23回の３Ｄグラフィクスを参照してください。つぎのようになります。Mathematicaと同様マウスドラッグでぐりぐりっと動かすことができます。（この図はスクリーンショットなので動かせません）

第20回 カオス

第20回カオス