2011/10/25

神話の翼を描く【カオスでCG】

先日サイエンス社の「CによるカオスCG」(川上博/上田哲史=共著)をAmazonで取り寄せた。

こちら初版が1994年と17年前のc言語を扱っており、ソース自体は古いが納められている内容はカオスCG向けのものばかりなので、大いに参考になった。

今回は息抜き程度にこの本の内容の一つ「神話の翼を描く」を実行した。

設定やソースをを説明する前に、まずは結果として描いた図形の一つを下に紹介する。


3枚の翼

3枚の翼 (a=-0.96,b=0.97)



3枚の翼が合わさったようななんとも不思議な図形である。

しかしこれを描くのに複雑な数式は必要ない。以下に示すたった2行の数式でこの図形は形成されるのである。

翼状図形を描く写像の例
 
 



上に示した図形はこの数式で(a=-0.96, b=0.97)として反復計算を行った結果をプロットしたものである。

非常に簡単だが、c言語のソースとgnuplotによって翼状図形を描く手順を以下に記す。

c言語ソース
/*神話の翼を描く*/
#include<stdio.h>
#include<math.h>

int main(void){

double x0=0.1, y0=0, x1, y1 ;//変数
double a, b ;//パラメータ
int i ;
FILE *fp;

if ((fp = fopen("tsubasa.dat", "w")) == NULL) {
printf("file open error!!\n");//tsubasa.datというファイルに出力する

}else{

printf("a > ");
scanf("%lf",&a);
printf("b > ");
scanf("%lf",&b);//端末からのパラメータa,bの入力

for(i=0;i<50000;i++){

x1=y0+a*x0-5.0/(x0*x0+1)+6.0+0.2*exp(-y0*y0);
y1=-b*x0; //数式部分
fprintf(fp,"%lf %lf\n",x1+y1,x1-y1); //現在のx,yを出力(描画のため引き伸ばし)
x0=x1;
y0=y1; //現在のx,yを次回の入力値にしてループ
}
}
fclose(fp);
return 0;
}


gnuplotでの描画手順
端末上で上記のプログラムを実行したあと、その場所でgnuplotを起動する。あとは以下のように入力していけば図形が描かれたファイル「tubasa.jpeg」が生成される。

$ gnuplot
gnuplot> set size square
gnuplot> set term jpeg
gnuplot> set output "tsubasa.jpeg"
gnuplot> plot "tsubasa.dat" with points pt 5 ps 0 lt 3


パラメータであるa,bの値を様々に変えることによって描かれる図形も異なってくる。では数例パラメータを変えて出来上がった図形を紹介しよう。

4枚の翼5枚の翼
a=0.02,b=0.97             a=-1.60,b=0.99

7枚の翼9枚の翼
a=-1.79,b=0.99             a=-1.844,b=0.97


順に4枚、5枚、7枚、9枚の翼が見て取れる。
最初にも述べたが、参考文献は「サイエンス社 CによるカオスCG」である。

神話の翼」と銘打っているわけだが、この程度の構造ならばごく普通の鳥の翼にも見られる。あの美しい羽の並びが、実は単純な法則から組み上がっているのかもしれないと考えるととても興味深い。

カオス理論では、「創発」という用語がキーワードの一つになっている。「創発」とは、「構成要素間の局所的な相互作用が系全体の大域的構造を生成することで、全体は部分の総和ではないという非線形現象の特性」のことである。

この翼のプログラムの例で言うと、xとyという2つの要素の前後の関係さえ与えてやれば、反復とともに要素が翼の形を生み出すということだ。出来上がった図のプロットのうち一つの点を取り出しても、それは座標上の1点に過ぎない。重要なのは点と点の間の相互関係である。

複雑系ーカオスにおいては、上の文で度々出てきた「相互作用」もキーワードのひとつとなっている。
系の要素を細分化して見ていくのではなく、系の要素間の相互作用を考えるのである。相互作用の法則さえ見つけてしまえば、全体の挙動を知ることが可能になるのだ。

今回の記事を書きながら私も、部分どうしの単純な相互作用が思いもかけず複雑な構造を生み出すのだということを実感した。今後もカオスCGを描きつつカオスを学んでいきたい。ではまた次の記事で。ありがとうございました。


スポンサーサイト
2011/10/11

Maximaでカオスを描こうという試み2

前回の記事では有名な「ローレンツ・アトラクタ」をMaximaという数式計算ソフトで描こうと試み、よくわからないので結局c言語で描いてお茶をにごしました。

まだMaximaの使い方を十分把握しておらず、ネット上の情報をつぎはぎしてもなかなかうまく行きませんでした。

しかし、試行錯誤を繰り返すうち、とりあえず前回の記事にあったローレンツアトラクタを描くことができました!

いきなりですが、その画像を下に示します。

maxima_atracter.jpeg



...まあ前回と同じようなグラフが描けたわけですが、今回はちゃんとMaximaを用いて描きました。
前回の記事を読んでいない方は以下のリンクで概要を確認できます。
Lorenz_atracter概要

では、Maximaでどのように操作すればよいのか順を追って説明します。

1.wxmaximaの導入
Maximaという大変便利な数式計算ソフトは、Ubuntu11.04には標準搭載されていません

Maximaは端末上で操作するソフトですが、数式が見づらいので「wxmaxima」をインストールすることをおすすめします。
「Ubuntuソフトウェアセンター」から「wxmaxima」を検索すると出てきます。

wxmaximaを起動すると、毎回、tipsを表示する小窓も現れますが、無視して閉じましょう。
表示されるのは真っ白なウィンドウであり、カーソルが点滅しないので不安ですが、すぐに文字を打つことが可能です。

ここで注意すべきことが1点あります。それは、実行のために「Enter」キーを押してもただ改行になってしまうということです。
maximaの端末上での実行は文の後にセミコロン「;」または「$」をつけて「Enter」ですが、wxmaxima上では「Shift + Enter」が実行の役割を果たします。最も良く使う操作なので注意しましょう。

2.wxmaximaでの実行

実行は以下の入力を1行ずつ行います。

(%i1) load("dynamics");
(%i2) f(x,y,z):=-10*(x-y);
(%i3) g(x,y,z):=-y-x*z+28*x;
(%i4) h(x,y,z):=x*y-8/3*z;
(%i5) rkdat:rk([f(x,y,z),g(x,y,z),h(x,y,z)],[x,y,z],[1,1,1],[t,0,50,0.01]);
(%i6) load(numericalio);
(%i7) write_data(rkdat,"rk.dat");


(%i3)」などは「input」を表しており、実行するたびに(%o3)の「output」が実際には表示されます。実際のwxmaxima上での見た目は以下の画像。

wxmaximaの画面

ではmaximaソースの説明をしたいと思います。

(%i1)では、"dynamics"というパッケージをロードしています。これは(%i5)で出てくる「rk」関数を使うために必要です。

(%i2)から(%i4)までは、前の記事で出てきたローレンツのアトラクタの連立微分方程式の右辺をそれぞれ関数f(x,y,z)、g(x,y,z)、h(x,y,z)とおいて定義しています。

(%5)で「rkdat」というものを「rk」関数のアウトプットと定義しています。
rk関数というのは、微分方程式を解くための「ルンゲ・クッタ法」を行ってくれる関数です。書式は、
rk([関数],[変数],[変数の初期値],[t,初期時刻,終了時刻,刻み幅])
というふうになります。

(%i6)で"numericalio"というパッケージをロードしています。次の(%i7)の「write_data」を用いるために必要です。

(%i7)この関数で、(%i5)で定義した「rkdat」を、"rk.dat"というファイルに出力しています。

3.gnuplotでの描画
できあがったdatファイルをgnuplotで描画します。手順は前回の記事と同様にします。
端末を開き、以下の流れを打ち込んでいきます。

$ gnuplot
gnuplot>set xlabel "x"
gnuplot>set ylabel "y"
gnuplot>set zlabel "z"
gnuplot>set view 75,30
gnuplot>set term jpeg
gnuplot>set output "maxima_atracter.jpeg"
gnuplot>splot "rk.dat" using 2:3:4 w l


これは前回とほぼ同じ操作ですが、最後の行だけ少し補足します。
using 2:3:4」というのはdatファイルの2,3,4列目をそれぞれx,y,zの値として読み込むということです。確認してみるとわかりますが、maximaで出力させた「rk.dat」というdatファイルは、1列目にtの値が入っています。使用したいのは2,3,4列のx,y,zの値なので、2:3:4のように書きます。

そして、その後についている「w l」ですが、これは「with lines」の略で同じ、プロットを線で結ぶという意味を表します。

以上で、ローレンツのアトラクタが描かれた画像ファイル「maxima_atracter」が出来上がります。


Maximaでカオスを描くという目的をなんとか遂げるできました。なかなか試行錯誤しましたが、次のカオス図形をMaximaで描くときはもっとスムーズにいけるだろうと思います。

今後の記事も「Mathematicaで学ぶシリーズ4 カオス入門」に載せられているMathematica向けのソースをMaximaに置き換えていく内容が多くなると思います。

ここまで呼んでくださりありがとうございます。ではまた次の記事で。


2011/10/10

Maximaでカオスを描こうという試み

先日、「Mahtematicaで学ぶシリーズ カオス入門」という本を購入した。

Mathemaicaで学ぶシリーズ4 カオス入門


これは、Mathematicaという、数式を計算するソフトを用いて様々なカオスにまつわるグラフを書いていこう、というコンセプトの書籍だ。

大学の先生に紹介してもらってその存在を知ることができた。

「Mathematica」は、団体向けの数式計算ソフトだ。個人向けのものもあるらしいが、十万円単位の価格とのこと。

とてもじゃないが手を出せないので、ネット上で無料で手に入れることのできる、「Maxima」というソフトをインストールして、Mathematicaの代用にした。

そして、いざ冒頭のMathematicaソースをMaxima用に置き換えようとして行き詰まった。ネット上をあちこちさまよい、情報を集めたがなかなかうまく行かない。Maximaに慣れればいいのだろうが、「これ、C言語で組んだ方がはやくないか?」というありさま。

以下にその問題の概要を示す。

ローレンツのアトラクタ

  

  

  

 

上記の連立微分方程式をMathematicaを用いて、グラフにプロットせよ。



結局のところ現時点でMaximaを用いて図を描くことはできていない。少し悔しかったのでc言語を用いて図を描いた。それを下に示す。

c言語で描いたLorenz_atracter



現時点ではMaximaの知識が少なく、明らかにc言語で組むほうが早かった。しかし、いずれ慣れてくれば、また問題が複雑になればMaximaで描いた方が早くなるかもしれない。

一応、グラフを描いたc言語のソースと、gnuplotでグラフを表示させる手順を述べる。

ちなみに私が使っているOSはUbuntu11.04である。Windowsユーザーの方も、c言語とgnuplotが使える環境であれば同様の操作で描けるのではないだろうか。もちろん端末からのコマンドは違ってくるだろうが…

c言語のソース

/*Lorenzのアトラクタ*/

#include

int main(void){
int i;

double x=1;
double y=1;
double z=1;//x,y,zの初期値を設定

double dx;
double dy;
double dz;
double dt=0.01;//時間tの刻み幅

double s=10;
double r=28;
double b=8/3;//定数部分

for(i=0;i<5000;i++){
dx=(-s*(x-y))*dt;
dy=(-y-x*z+r*x)*dt;
dz=(x*y-b*z)*dt;//微分方程式

x+=dx;
y+=dy;
z+=dz;//変化量を足してループ

printf("%f %f %f\n",x,y,z);//毎回x,y,zを出力
}
return 0;
}




このソースをコンパイルしたのち以下の手順でgnuplotにより描画した。

まず端末を開き以下のように入力する。
$ ./a.out -> Lorenz_atracter.dat
$ gnuplot

1行目はそのままだと端末上に表示される実行結果を"Lorenz_atracter.dat"というファイルに出力する操作。
2行目でgnuplotを起動。


以下gnuplotの操作。

gnuplot>set xlabel "x"
gnuplot>set ylabel "y"
gnuplot>set zlabel "z"

これで、出力したときにx軸、y軸、z軸にそれぞれx,y,zの名前がつく。


gnuplot>set view 75,30

これは、3次元グラフを出力したときの見る角度を操作するコマンド。デフォルトでは[set view 60,30]となっている。ローレンツアトラクタの広がりをよく見るために少し調整した。


gnuplot>set term jpeg
gnuplot>set output "Lorenz_atracter.jpeg"
gnuplot>splot "Lorenz_atracter.dat" with lines lt 3

上の2行で、jpeg形式で出力することを宣言し、出力先の"Lorenz_atracter.jpeg"というファイルを作成し、開いている。

そして3行目で出力をし、これでプロットされたグラフができる。
[splot "ファイル名"]で3次元プロット、[with lines]でプロット点を線で結ぶ、[lt 3]で、線の色の指定(1は赤色、2は緑色、3は青色・・・というように番号によって色が変わる)という操作になっている。

ちなみにこの上の2行を飛ばして入力すれば、gnuplotの別ウィンドウが開き、グラフの確認をすることができる。

以上で操作の説明は終わりである。


冒頭で述べた、「Mahtematicaで学ぶシリーズ カオス入門」という本は、様々なカオスにまつわるグラフを描画するのがコンセプトで、カオス理論について詳しく説明しているわけではない。しかし、本のタイトルにもある通り、カオス理論の入門者にとってみれば、カオスの全体をコンパクトにまとめてくれているので非常にありがたい。

この本をもとにカオス理論を学ぶ骨組みをつくっていくことができそうだ。



Maximaで描画する方法についてはまた後日記述したいと思う。

ではまた次の記事で。
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。