2012/03/14

チョーサ氏とゴルビツキー氏のカオス

前回、前々回とほぼ同じ体裁の記事となりますがご了承ください。

さて、ここ3回の記事では「繰り返し公式」を用いてカオスを生成し描画する、ということをやってきました。今回で一区切りとします。繰り返し公式を用いてカオスを生成することについての詳しい説明は「チョーサ氏とゴルビツキー氏のカオス」の描画のあとにしたいと思います。

チョーサ氏とゴルビツキー氏の対称性をもつカオス

チョーサ氏とゴルビツキー氏は複素数の繰り返し公式を用いて、対称性のあるカオスを生成する公式を導きました。その公式は以下のようなものです。

複素数  z=x+iy (x,yは実数) とし、mを自然数とする。



ここに、



である。これを用いて、漸化式



を考える。この公式には大雑把に言ってm個の軸を持つ対称なカオスが見られる。
ここではm=3の場合について、カオスを描いてみる。





であるから、実数部xと虚数部yに分けて考えると漸化式は次式のようになる。





ここに、



である。



この公式のm=3の場合の繰り返し計算を行うc言語のプログラム(chosa.c)のソースは以下のようになります。


/*チョーサ氏とゴルビツキー氏の対称性をもつカオスを描く*/
#include<stdio.h>
#include<math.h>

int main(void){

double x0=0.1, y0=0.1, x1, y1,A ;//変数(x0とy0はともに0ではいけない)
double a, b,c,d ;//パラメータ
int i ;
FILE *fp;

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

}else{

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

for(i=0;i<50000;i++){
A=a*(x0*x0+y0*y0)+b*x0*(x0*x0-3.0*y0*y0)+c;
x1=A*x0+d*(x0*x0-y0*y0);
y1=A*y0-2.0*d*x0*y0; //数式部分

fprintf(fp,"%lf %lf\n",x1,y1);//現在のx,yを出力

x0=x1;
y0=y1;
}
}
fclose(fp);


return 0;
}



このプログラムを実行したあと、gnuplotで以下のように入力します。

$ gnuplot
gnuplot> set size square
gnuplot> set xrange [-2.0:2.0]
gnuplot> set yrenge [-2.0:2.0]
gnuplot> set term jpeg
gnuplot> set output "chosa.jpeg"
gnuplot> plot "chosa.dat" with points pt 5 ps 0 lt 3

xとyの範囲をxrange[-2.0:2.0]などと指定しているのは、出来上がった図が綺麗に対称に見えるようにするためです。適宜変更して調整します。(ちなみに最後の行の [lt 3]とある部分の数字を変えることで、図の色を決めることができます。環境によって色々違う模様)

以上よりできあがった図はこちら。クリックで拡大できます。

chosa1.jpegchosa2.jpeg
a=-1.0,b=0.05,c=2.275,d=-0.5      a=1.0,b=0.0,c=-2.25,d=0.2

chosa3.jpegchosa4.jpeg
a=1.0,b=0.0,c=-1.9,d=0.4        a=-1.0,b=0.1,c=1.6,d=-0.8

どの図もおおむね3つの軸に対称な形をした美しいカオスとなっています。


さて、ここ3回の記事で繰り返し公式によるカオスを見てきました。生成されたカオスを図に表してみるといずれもとても美しく、ただ法則性を指定してパソコンに描かせただけのようにも思えます。そもそもカオスとはいったい何なのか?これについて、私が参考とした「CによるカオスCG サイエンス社 川上博/上田哲史」の説明を引用することにします。

カオスは、過去およそ20年くらいの間(この本が書かれたのは1994年)、様々な分野で興味がもたれ研究されてきた。しかしまだまだ未知の部分も多い。そもそも「カオス」という言葉自身、まだ定まった定義ができていないのが現状である。

●カオスとは
この本では、コンピュータにより繰り返し公式を用いて点列を生成し、次の2つの点が確認できるとき、得られた点列をカオスとみなした。

(1) 周期倍分岐列の連鎖があり、その極限を過ぎた値にパラメータが選ばれていること。

(2) 点列は平面の有界領域を落ち着くことなく彷徨していること。
この性質を持つ点列の集合は、現在もっとも自然にカオスとして受け入れられている。欲をいえば、

(3) 点列の運動が不安定であること。すなわち、運動が初期値の微小な変化に伴って鋭敏に変化すること。
を確認しなけらばいけない。これにはエノン写像の実験ですこし述べた(ブログでは紹介していない)リヤプノフ指数を計算するとよい。本文では、(2)の性質を観察することで(3)の性質をもつとみなしてごまかしてしまった。

点列の動きがランダムに見えることは、カオスの中にうまっている周期点に出入りする2本の曲線の絡み具合によって生じている。数値計算の中に含まれる打ち切り誤差や丸め誤差は、この絡みによって、特定の方向に助長される。これが初期値の近い2つの点列を長時間後には別々に行動させる原因となっている。



枠内で青字で示した用語の説明をします。

周期倍分岐
過去の記事:ロジスティック写像の分岐図で示した図があります。これは横軸にパラメータaの値をとっているのですが、aの増加にともないxの収束先の周期点の周期が1,2,4,8・・・のように増加していきます(図で言うと、枝が1,2,4本・・・と増えて行っている)。このように周期点の周期が倍々に増えていく現象を周期倍分岐と呼びます。図ではa=3.5699以降周期が無限大となり、カオスが生まれています。

bunkizu.jpeg
ロジスティック写像の分岐図


リヤプノフ指数
力学系においてごく接近した軌道が離れていく度合いを表す量のことです。(Wikipedia参照)
詳しい式は少々複雑なので、今後の記事で取り上げることにします。この指数が負の時、軌道は安定な固定点または周期点をとり、0のとき準周期点、正のときカオス点列をとります。

周期点
ある点が時間とともに2点を行ったり来たりする軌道をとる場合、これを安定な2周期点と呼びます。周期的な運動がとる点のことです。運動がカオス状態の場合、周期点が無限にあることになります。


難しい内容となってきました。結論的にカオスって何の役に立つの?と思うかもしれません。しかし、カオスは自然のランダムさを生み出しているのが実は簡潔な法則だったということを示す、画期的な発見なのです。

今回の記事はここまでです。次回以降、カオス理論の用語の説明、あるいはカオスに携わった人の紹介など考えています。



2012/03/05

池田氏のアトラクタ

さて、久々の記事となりましたが前回の続きで「池田氏のアトラクタ」を描いていきたいと思います。

池田氏のアトラクタ

以下の写像のカオス的アトラクタを計算し、描画する。









上記の繰り返し計算を行うC言語のプログラムソース(ikedaatracter.c)はこちら。

/*池田氏のアトラクタを描く*/
#include<stdio.h>
#include<math.h>

int main(void){

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

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

}else{

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

for(i=0;i<50000;i++){
theta0=x0*x0+y0*y0-2.0;
x1=a+b*(x0*cos(theta0)-y0*sin(theta0));
y1=b*(x0*sin(theta0)+y0*cos(theta0));//数式部分

fprintf(fp,"%lf %lf\n",x1+y1,x1-y1);//現在のx,yを出力(描画のため引き伸ばし)

x0=x1;
y0=y1;
}
}
fclose(fp);
return 0;
}



端末上で上記のプログラムを実行し、いつもと同じようにgnuplotで描画します。
$ gnuplot
gnuplot> set size square
gnuplot> set term jpeg
gnuplot> set output "ikedaatracter.jpeg"
gnuplot> plot "ikedaatracter.dat" with points pt 5 ps 0 lt 3


すると、以下のような渦巻き模様が描けます。

ikedaatracter0.jpeg



パラメータの値を変えていくと、若干渦巻きの密度が違った図が描けます。

ikedaatracter.jpegikedaatracter1.jpeg
a=4.0, b=0.4               a=4.0, b=0.8


このアトラクタを生み出す上に示した式は、レーザー光学のモデルとして池田研介氏によって研究された写像です。
ヤコビ行列式はJ(x,y)=b^2となるので、|b^2|<1の場合縮小写像となります。また、θがx^2+y^2の関数であれば、ヤコビ行列式は一定b^2となるので、θに他の関数を選んでも形の違ったカオスを描くことができます。興味があれば、是非他のθを使った図を描いてみてください。

縮小写像について

写像によって、移った先にできる図形がもとの図形の面積より小さくなるような写像。


本当はもう少し詳しく記述したかったのですが、それは次回の「チョーサ氏とゴルビツキー氏のカオス」の記事で、またはその次の記事で書こうと思います。
2012/01/20

2枚の翼とカオスとなっているアトラクタ

今回と次回で3つのカオスCGを描きたいと思います。

1つ目は「2枚の翼」です。これは過去記事:神話の翼を描く【カオスでCG】で紹介したカオスCGと同系統の図です。鳥が羽ばたいているようにも見える不思議な図形です。今回はこちらを描きます。

2つ目と3つ目は池田研介氏のアトラクタとチョーサ氏とゴルビツキー氏のアトラクタによるカオスCGです。前者は渦を巻くような図形が、後者は対称的な図形が現れ、とても綺麗です。次回はこちらを描きます。

では、この3つのカオスを描いていきます。参考にしたのは、「CによるカオスCG」 川上博・上田哲史共著 サイエンス社 です。




2枚翼の写像

以下の繰り返し公式にて2枚羽根のカオスを求めます。




ただし、a,b,Rは定数で、nは0以上の整数となります。

上記の繰り返し計算を行うC言語のソース(2maitsubasa.c)はこちら。

/*2枚の翼を描く*/
#include<stdio.h>
#include<math.h>

int main(void){

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

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

}else{

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

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

x1=y0+a*x0+5.0*x0*x0/(1+x0*x0)+1.0-0.2*exp(-y0*y0);
y1=-b*x0+R*pow(-1.0,i); //数式部分
if(i%2==0){ //周期2なので2周ごとにプロット
fprintf(fp,"%lf %lf\n",x1+y1,x1-y1); //現在のx,yを出力(描画のため引き伸ばし)
}
x0=x1;
y0=y1; //現在のx,yを次回の入力値にしてループ
}
}
fclose(fp);
return 0;
}

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

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



2maitsubasa.jpeg
a=0.02,b=0.98,R=5.0として描画したときの2枚翼


私はこの図形を見たときにまるで鶴のようだ、と思いました。この記事を読んでくださっている方はどう思うでしょうか。

ソースの中にポツンと「周期2なので2周ごとにプロット」というコメント文があるのに気づかれた方もいらっしゃると思います。これは、上で示した2つ目の数式の中にR(-1)^nであるところの強制項があるためです。強制振動の微分方程式での外力みたいなものでしょうか。

ともかく、この強制項が2の倍数ごとに同じものになりますので、(n=0,2,4,...でR、n=1,3,5,...で-R)これをそのままプロットすると翼が2つ重なったものになってしまうのです。そのため、nが偶数または奇数のときだけをプロットすれば重ならない図を描くことができます。(ソース中ではnの役割をiで担っており、iが偶数の時プロットするようにしています。

重なって描かれた場合どうなるかを下に示します。

2maitsubasa2.jpeg
周期2による2つの2枚翼の重なり


先に示した図はこの図でいうところの下側の翼にあたります。上側の翼はnが奇数のときだけプロットした場合に現れるものです。2羽の鶴が交差した瞬間のようにも見え、これはこれで面白い図だと思います。

図の描き方ばかりの説明になってしまいましたが、これの一体どこがカオスなんだ?ということについて簡単に説明したいと思います。


キーワードは「アトラクタ」です。

ここに描いた図形は「ストレンジ・アトラクタ」と呼ばれるもののひとつです。アトラクタとは、平面上で言えば、周りの点を引き寄せる(attract)もののことです。アトラクタによる様々な平面図形は、総じて周りの点を引き寄せているような構造をしています。

さて、カオスの定義の一つとして「初期値に対する鋭敏な依存性」が挙げられます。すなわち、最初の点をごくわずかにずらしただけでも、将来的には決定的な差異が生じるという性質です。

このストレンジ・アトラクタも、パラメータa,b,Rの値をほんの僅かに変えるだけで、50000回反復後にはかなり異なる図形が描かれます。初期値に対する鋭敏な依存性、があるというわけですね。カオスです。



次回、冒頭で述べましたように、池田研介氏のアトラクタとチョーサ氏とゴルビツキー氏のアトラクタによるカオスCGを描きたいと思います。


2012/01/18

ロジスティック写像の分岐図

関数族において、十分時間が経過した後の写像の繰り返しから求めた分岐をパラメータの関数として表示したものを、分岐図(bifurcation diagram)という。分岐図の例として、以下の例を示す。




ロジスティック写像の分岐図

関数族は
   
で与えられる。

以下にこの計算を繰り返していくC言語のプログラム(logistic.c)を示す。

#include

int main(void){
double a,x=0.8;
int i,j;
for(a=2.0;a<4;a+=0.002){
for(i=0;i<100;i++){
x=a*x*(1-x);
}
for(i=0;i<100;i++){
x=a*x*(1-x);
printf("%f %f\n",a,x);
}
}

return 0;
}


こちらをubuntu端末上にて
$gcc logistic.c
$./a.out ->log.dat
$gnuplot

gnuplot> plot "log.dat" pt 5 lt 0 ps 0

と打ち込み表示された図をここに載せる。

bunkizu.jpeg

この図には、分岐の始まるところの周囲を含めた部分を拡大していくと、もとの図とまったく同じような形が見えてくるフラクタルな構造が見られる。
他にも様々な特徴が見られるのだが、今後の記事で紹介していきたい。
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。