カオス #34 相関次元

今宵は相関次元法というのを用いてフラクタルな図形の次元を求めてみます．
相関次元法では，以下で定義される相関積分を求めます．



H(x)は以下で定義されるHeaviside関数です．


ただし，P₁, P₂, …, P_Nは，
P_i=f(P_i-1)
という関係にあります．
また，|P_j-P_i|は，m次元空間内でのP_iと，P_jの距離で，以下に定義します．



rを変化させ，log₂(ε)と，log₂(C(r))をプロットして行き，グラフの傾きを読み取ることにより，相関次元を求める事が出来ます．

図で式の意味を示すとこんな感じ．


Fig.1

黒丸 - P_j
白丸 - {P₁, P₂, …, P_N}
円の半径 - ε

H(ε-|P_j-P_i|)は，P_jを中心とする半径εの円(球)にP_iが入っていれば1を，入っていなかったら0を返すというわけです．

■計算例

・Hénon写像

x_t+1=1-ax_t²+by_t
y_t+1=x_t
a=1.4
b=0.3


Fig.2

logε∈[-6,-4]について最小自乗法により計算した相関次元
D=1.23753

畑政義写像

z_n+1=f₁(z_n) or f₂(z_n)

f₁(z)=az+bz~
f₂(z)=c(z-1)+d(z~-1)+1
(z=x+yi => z~=x-yi)

・以前，次元が約1.736であるといったフラクタル図形
a=0
b=0.45+0.5i
c=0
d=0.45-0.5i


Fig.3

D=1.74996

・Koch曲線
Koch曲線のハウスドルフ次元はlog4/log3≒1.26186ですが果たして…．
a=0
b=0.5+0.3i
c=0
d=0.5-0.3i
※以前紹介したパラメータと違いますが，それも含めて多分これらのパラメータは正確で無いです．


Fig.4

D=1.27587

・C曲線
C曲線のハウスドルフ次元は実は2です．
a=0
b=0.5+0.5i
c=0
d=0.5+0.5i


Fig.5

D=1.94509

と，こんな感じです．

カオス #33 Star Julia

Star Juliaとはクリスマスツリーネタで使いましたが，参考文献にあるIFSです．
説明が手抜きなのは，リンク先にある漸化式が若干謎な為です．
実際には下のようなソースで動きます．

double xa=-4*x+Math.cos(m*w);
double ya=-4*y+Math.sin(m*w);
double r=Math.sqrt(xa*xa+ya*ya);
w=Math.atan2(ya,xa);
double x1,y1;
if(mt.nextBoolean(0.5)){
 x1=-Math.cos(n*w)-Math.sqrt(r)*Math.cos(w/2)/2;
 y1=-Math.sin(n*w)+Math.sqrt(r)*Math.sin(w/2)/2;
}else{
 x1=-Math.cos(n*w)+Math.sqrt(r)*Math.cos(w/2)/2;
 y1=-Math.sin(n*w)-Math.sqrt(r)*Math.sin(w/2)/2;
}
x=x1;
y=y1;

例:
某氏より要望があった為，大きいサイズで．

m=1, n=1


Fig.1

回転させて，重ねると雪になります．


Fig.2

m=2, n=2


Fig.3

こちらは星．

折角なので壁紙サイズ．


Fig.4


Fig.5

参考文献
Paul Bourke．“Star Julia”．＜http://local.wasp.uwa.edu.au/~pbourke/fractals/starjulia/＞．Paul Bourke - Personal Pages．（参照2010年1月13日）

カオス #32 Clifford Attractor

再々びClifford Attractor．そのうちまたいつかもしかしたらいつぞやに作ります．


Fig.1 a=1.1, b=-1.1, c=1.0, d=1.4


Fig.2 a=0.7, b=-1.8, c=1.6, d=0.9


Fig.3 a=0.7, b=-2.0, c=1.6, d=0.9

カオス #31 Langton's ant

「Langtonのアリ」です．
Christopher Langtonが考案しました．概要は以下．

格子状の平面を構成する．各マスは黒(=0)に初期化．ある1つのマスにアリを配置．アリは各ステップで上下左右いずれかのマスに移動できる．アリの移動規則は以下に従う．

・黒いマスにアリがいた場合
90°右に方向転換，そのマスの色を反転させ，1マス前進．
・白いマスにアリがいた場合
90°左に方向転換，そのマスの色を反転させ，1マス前進．

この規則でアリはかなりカオスな挙動を示します．
しかし，10000程歩いた後，真っ直ぐな(斜め方向の)「道(highway)」を作り始めます．この道は，LangtonのアリのAttractorかもね，らしいです．

Case 1:
格子:1024x1024
アリの数:1
初期条件:
位置:(512, 512)
方向:上


Fig.1 step=0
白い点がアリです．


Fig.2 step=1000
小さい範囲でゴニョゴニョと．赤いマスは黒を反転させたということです．


Fig.3 step=5000
範囲が大きくなりましたが，直進する気配がありません．


Fig.4 step=10000
10000step経ちました．


Fig.5 step=10500
こんな，


Fig.6 step=11000
感じで，


Fig.7 step=11500
のんびりと，


Fig.8 step=12000
のびていきます．


Fig.9 step=500000
長時間経過後．

Case 2:
格子:512x512
アリの数:3
初期条件:
位置:(128, 128), (256, 256), (384, 384)
方向:上, 上, 上


Fig.10 step=0
3アリ配置したナリ．


Fig.11 step=10000
ここまでは一緒．どのアリが反転したかによって色を変えたナリヨ．


Fig.12 step=50000
カオス．


Fig.13 step=500000
サイケ．

参考文献

Wikipediaの執筆者たち．“ラングトンのアリ”．＜http://ja.wikipedia.org/wiki/%E3%83%A9%E3%83%B3%E3%82%B0%E3%83%88%E3%83%B3%E3%81%AE%E3%82%A2%E3%83%AA＞．Wikipedia．（参照2009年12月23日）

カオス#30 リアプノフ指数

リアプノフ指数(Lyapunov exponent)

リアプノフ指数とは，力学系において接近した軌道が離れていく度合いを表す量を表します．
1次元写像に関して，この値は非常に簡単に求めることが出来ます．
具体的な式は以下の通り．



例:

テント写像

テント写像T(x)は以下で定義されます．



パラメータaを0から2まで変化させたときのリアプノフ指数を計算してみます．



Fig.1

0≦a<1 ⇒ λ>0
1≦a≦2 ⇒ λ<0
となっている所がポイントです．

ここで，テント写像の分岐図を見てみましょう．


Fig.2

分岐図からは，a>1のとき，カオスになっていることが分かります．
これより，(強引だけど)λ>0の時，写像はカオス性を持つと予想出来ます．

ロジスティック写像

ロジスティック写像L(x)=ax(1-x)において，aを変化させた時のリアプノフ指数を求めてみます．


Fig.3 0≦a≦4
3.5≦a≦4が特にゴチャゴチャしているので，拡大してみます．



Fig.4 3.5≦a≦4

所々λ<0となっています．

さて，ロジスティック写像の分岐図はFig.5です．


Fig.5

これらを比較してみると，λ<0の時，分岐図はカオスではありません，ある周期の周期軌道に落ち込んでいます．

そんなわけで，リアプノフ指数は写像がカオス性を持つかどうかを判断するのに非常に役立つといえるでしょう．

カオス #29 クリスマスツリーの製作

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

クリスマスが今年もやって来る
悲しかった出来事を 消し去るように
さあ パジャマを脱いだら 出かけよう
少しずつ白くなる 街路樹を駆け抜けて

クリスマスがもうじきやって来る
うれしさをかくせない 犬や猫まで
もみの木に灯る明かりをうけて
いつもよりやさしそうな パパの目が笑ってる

クリスマスは誰にもやって来る
もしひとりぼっちでも 淋しがらずに
心に住むサンタに呼びかけて
幼い頃の夢を 思い出してごらんよ

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

すてきなホリデイ/竹内まりや

やってきましたクリスマス．
今年は「カオスで作るクリスマスツリー」です．

■木
ツリーは以前作ったシダのパラメータを変えることで表現してみます．

初期座標を(x₀=0, y₀=0)とし，以下に示す4つの座標変換の内一つを選択し，反復適用していきます．Pは選択される確率です．

x_n+1=-0.08x_n
y_n+1=0.72y_n-1
P₁=1%
x_n+1=-0.16x_n-0.32y_n
y_n+1=0.4x_n-0.18y_n+1.6
P₂=7%
x_n+1=0.16x_n+0.32y_n
y_n+1=0.4x_n-0.18y_n+1.6
P₃=7%
x_n+1=0.86x_n
y_n+1=0.91y_n+0.65
P₄=85%


Fig.1

■星，雪
これらはStar juliaというフラクタルを用いてみます(詳細は後日)．

・星
m=2, n=2

Fig.2

・雪
m=1, n=1

Fig.3
このパラメータだと三角形が出来あがるので，60°ずらしたものを重ねています．

■クッキー
Gingerbread man(カオス #7 GingerbreadMan参照)が正に其のものです．


Fig.4

■飾り
いいネタが無いので，良く分からない飾りを幾つか．
グモウスキー・ミラ写像というものを使いました(詳細は後日)．


Fig.5
a=0
m=0.302


Fig.6
a=0
m=-0.36


Fig.7
a=0
m=-0.237

■仕上げ
これらをいい感じに組み合わせます．


Fig.8

完成．カオスとともにクリスマス・イブ…それは予測不可能の中にある秩序を見出す喜び…

カオス #28 Clifford Attractor

再びClifford Attractor．少しサイケな世界をお楽しみください．


Fig.1 a=1.3, b=1.7, c=0.5, d=1.4


Fig.2 a=-1.5, b=-3.0, c=1.5, d=1.5


Fig.3 a=1.4, b=-1.4, c=0.5, d=0.8

カオス #27 畑政義写像

以前，大量のアトラクタを紹介した畑政義写像でしたが，今回は，それらをアニメーション化した物を．
数々のアトラクタのパラメータを連続的に変化させることで，非常に面白い動きを見ることが出来ます．

カオス #26 馬蹄写像

楽しい写像の御提供．

馬蹄写像

馬蹄写像は下のような2変数写像です．固定された数式では無いので，定義域，値域をそれぞれ図示しています．


Fig.1 馬蹄写像
2つの半円D1,D2と正方形Sが定義域です．図では赤のチェックの領域です(一部隠れています)．
馬蹄写像を定義域全体に対して適用すると，定義域を折り曲げたような形になります．図では青のチェックの領域分です．この値域の形が馬蹄に似ているから馬蹄写像と呼ぶのです．

この写像を2回適用するとFig.2のようになります．


Fig.2

3回適用．


Fig.3

繰り返すと倍々に折り重なっていくのが分かります．また，D1，D2内に入った点は以後そこから出られなくなります(D1内の点はD1内へ，D2内の点はD2内へ写像するから)．そして最終的にD1，D2内に入る点は到る所無数にあるので，初期値の極僅かな差によって最終的な挙動が全く変わります．ロジスティック写像の2次元版みたいな感じです．

疑似コードはこんな感じ．上にも書いた通り，固定された数式ではないので，別にこのような定義でなくともかまいません．

horseshoe(x, y)
  if(-1/2≦x≦1/2 && -1/2≦y≦1/2) # Sに含まれる
    if(-1/2≦x≦-1/8)
      u=map(x, -1/2, -1/8, -5/8, 1/2)
      v=map(y, -1/2, 1/2, 1/8, 3/8)
    else if(1/8≦x≦1/2)
      u=map(x, 1/8, 1/2, 1/2, -5/8)
      v=map(y, -1/2, 1/2, -1/8, -3/8)
    else # -1/8≦x≦1/8 Uの字に曲がる所
      r=map(y, -1/2, 1/2, 1/8, 3/8)
      d=map(x, -1/8, 1/8, π/2, -π/2)
      u=r・cos(d)+1/2
      v=r・sin(d)
    return (u, v)
  else if(dist(-1/2, 0, x, y)≦1/2) # D1に含まれる
    u=map(x, -1, -1/2, -3/4, -5/8)
    v=map(y, -1/2, 1/2, 1/8, 3/8)
    return (u, v)
  else if(dist(1/2, 0, x, y)≦1/2) # D2に含まれる
    u=map(x, 1/2, 1, -5/8, -3/4)
    v=map(y, -1/2, 1/2, -1/8, -3/8)
    return (u, v)
  else # 定義域外
    return null

map(value, low1, high1, low2, high2)
  return (value-low1)(high2-low2)/(high1-low1)+low2;

dist(x1, y1, x2, y2)
  return sqrt(sq(x2-x1)+sq(y2-y1))

写像する様子を映像化するとこんな感じ(イメージです)．