模造5(バタフライエフェクト)
このページはマス旅の一部です。
今回は「フラクタルとカオス」のシミュレーションです。
フラクタルとカオスは美しく不思議な画像と関係が深い。
フラクタルは自己相似のことで、自然物に近い美しさが光りますね。
カオスはローレンツが使った「バタフライ効果」のような初期値への敏感性からの混沌状態や、過渡状態から定常したアトラクタ状態へと落ち着く現象などを研究する領域です。
単純に「カオスやねん」というのは日常のコトバではよくある。
課題1:バーンスレイのシダの図を行列計算でかこう。
発生確率つきのアフィン変換(行列積+移動)f1、f2、f3、f4を用意して反復すると
シダの形がかける。このことを数学者のバーンスレイが示したことは有名だ。
k=0,(x_k,y_k)=(0,0)で初期化する。
関数fnの行列計算を
f_n(x,y)=F_n x+ An、発生確率p_nを使って、列記しょう。
F_1={{ 0, 0},{ 0, 0.16}}, A_1=t{0, 0 } ,p_1 =0.01 #茎
F_2={{ 0.85, 0.04},{-0.04, 0.85}}, A_2=t{0, 1.6} ,p_2 =0.85 #小さい葉
F_3={{ 0.2,-0.26},{ 0.23, 0.22}}, A_3=t{0, 1.6} ,p_3 =0.07 #大きい左葉
F_4={{-0.15, 0.28},{ 0.26, 0.24}}, A_4=t{0,0.44} ,p_4 =0.07 #大きい右葉
これを関数形にすると、
f_1={ 0, 0.16}
f_2={ 0.85x +0.04y ,-0.04x +0.85y + 1.6}
f_3={ 0.2x -0.26y , 0.23x +0.22y + 1.6}
f_4={-0.15x +0.28y , 0.26x +0.24y +0.44}
関数fは、リストfs、psに関数と確率を入れておき、
psの確率に応じてfsのどれかを選択するしくみだ。
k=k+1
(x_k,y_k)=f(x_k,y_k)
計算の結果はxが-2.2と2.7の間、yが0以上10以下になる。
<コード化>
「作成」ボタンのクリック時のスクリプトに貼り付けます。
Fs ={{{0, 0}, {0, 0.16}},{{0.85, 0.04},{-0.04, 0.85}},{{0.2, -0.26}, {0.23, 0.22}},{{-0.15, 0.28}, {0.26, 0.24}}}
As ={{0, 0},{0, 1.6},{0, 1.6},{0, 0.44}}
Ps = {0.05, 0.75, 0.1, 0.1}
step = 0
P = {0, 0}
Track = {{0, 0}}
text1 = "ステップ数 = " + step + " | 打点数 = " + Length(Track)
shida=Zip(Point(tr),tr,Track)
#点のサイズを3にします。
「更新」ボタンクリックのスクリプト
SetValue(step, step + 1)
num = RandomDiscrete({1, 2, 3, 4}, Ps)
F = Element(Fs, num)
A = Element(As, num)
Fk = F Vector(Point(P))+Vector(Point(A))
SetValue(Track, Append(Track, {x(Fk),y(Fk)} ))
SetValue(P, {x(Fk),y(Fk)} )
「リセット」ボタン
SetValue(step, 0)
SetValue(P, (0, 0))
SetValue(Track, {(0, 0)})
タイトル「打点数1000でシダが見える」
<振り返り>
処理によってデータ形式を変えました。
計算処理は行列で、
点データ保存はリストのリストで、
視覚化は点オブジェクト化で。
メモリの負担を減らすことで1000点くらい打てるようになります。
また、点Fkが動くことで、F*P+Aの動きがわかります。
先まで上ると下からまた上へと打点することを繰り返します。
それが、どのラインに行くかの確率を変えることでシダの形状が決まるのです。
打点数が万単位であれば、確率分布はバーンスレイと同じでいいです。
しかし、打点が100で茎1点、900で茎9点では、茎が見えず葉だけに
なります。それをさけるために、Ps = {0.05, 0.75, 0.1, 0.1}
としました。
すると、打点が900くらいでもかなり茎まわりが太くなり植物らしくなるでしょう。
スライダーを貼り付けて、最新情報に更新スクリプトを貼り付けても
更新が進みますが、打点数とスライダーのカウントがあいません。
これは、商品ではありません。だから、
アナログなおもちゃのように地味に打点することで「ロジックに関心を持つ」
という方向で止めておいてもよいでしょう。
打点数1000でシダが見える
気象モデルの研究者ローレンツがちょとした初期値のちょっとしたちがいでトーラスが変形していき、ずっこけた奇妙な定常状態(ストレンジアトラクタ)になることを発見した。だから、その図形を発見者にちなんでローレンツアトラクタと呼ぶようになったんだね。ローレンツアトラクタを発端に初期値敏感性というカオスの特徴が注目されましたが、もっと手軽に初期値敏感性を感じることができる素材がある。
課題2:ロジスティク方程式で初期値敏感性を感じよう。
x_{k+1}=a x_k(1-x_k)という数列で、aは正の実数
aの1, 2, 3, 3.449489, 4を境目として、数列xkの挙動が変わる。
xは1でない正数とするが、xを0.000001変えただけグラフの形状が変わる。
<コード化>
画面をクリックして、アプリの上位メニューの下部にある設定を選ぶ
設定>グラフィックスビユー>次元>
比率>x軸:y軸=100:1
次元>x最小-3,x最大51、y最小-0.2,y最大1.5
設定>一般>丸め>小数点以下10桁
「作成ボタン」クリック時に貼り付ける。
Anum = Slider(1,7,1)
As = {1.2, 2.8, 3.1, 3.5, 3.7, 3.9, 4.0}
a = Element(As,Anum)
x0 = Slider(0.001, 0.002, 0.000001)
f(x)= a x(1-x)
xx=Sequence(k,k,0,50)
yy=IterationList(f,x0,50)
logis=Zip((p,q),p,xx,q,yy)
PolyLine(logis)
#f,logisは非表示にする。
aが大きいとき、x0をわずかに増やしただけ数列xnの挙動が激変するのがわかるね。
これがカオスだ!
カオスの不気味さを際立たせるために、画面の背景を黒、グラフを黄色、スライダーをそれぞれ水色と明るいピンクにする。x0をアニメーションオンにすると、地震や心電図の波形のようなグラフが神経質に
変動するのが見える。
タイトルは「これがカオスだ!(ロジスティック方程式)
これがカオスだ!(ロジスティック方程式)
曲線的なイメージが続いたので
最後に直線的フラクタルを作って心を整えます。
5回のテーマを通して、シミュレーションそのものそうですが、更新のデザインパタンも
いろいろ学んだので、その総集編となる課題を用意しました。
まとめ課題:シェルピンスキーのギャスケットを数値ではなく「点」で表示しましょう。
まなんだことを思い出して取り組んでみよう。
模造3(状態の変化)では、01のリストで結果を表示していました。
今回は、これを点集合にします。
また、サイズを多くして、ギャスケットがきれいに見えるようにしたいですね。
1行の列数をLeに入れます。1をはさんで50個の0が並ぶ最初の行のサイズが、
Le=50+1+50=101です。
また、あらかじめ0が並ぶだけの行を作ってそれを上書きしてました。
最初の行Cを更新したら、行単位でCsに追加してCsを上書きすることで、
作成とリセットがカンタンになりますね。
Csは更新するたびに行リストが増える、リストのリストです。
Csの各行を読み込んだ01の状態リストと、1から順にx座標を入れただけのリストpxをかけます。
この2つのリストの同じ位置どうしの積をZipで求めpmulとましょう。
pmulはCsと同じ行数になりますね。
pmulの各行リストでは0が点のない位置で、0でない数は点のx座標になります。
そこで、pmulのp行q列が0でないときだけ、その値をx座標として読み取り、その行番号にマイナスをつけた座標点PntをPoint({pmulのp行q列の値、-p)})としましょう。
この点生成をp行を1から更新行位置cntまで、q列を1からLeまでSequenceを2重にして実行します。
これは通常言語でのfor文の入れ子と同じ効果があるね。
Cの次の行、つまり更新のロジックは模造3(状態の変化)の振り返りに書いたように、
LとRのデータを作り、排他的論理和と同値になる剰余計算でnextCを出せばよいですね。
<コード>
# 初期セル1をはさんで左右に50個の0を並べる
#「作成」ボタン
cnt=1
Le=101
C={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
Cs={{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}
Rule={0,1,0,1,1,0,1,0}
CL = Join({Last(C), Take(C, 1, Le - 1)})
CR = Join({Take(C, 2, Le), First(C)})
px=Sequence(k,k,1,Le)
pmul=Sequence(Zip(p*q,p,Element(Cs,k),q,px),k,1,cnt)
Pnt=Sequence(Sequence(If(Element(pmul,p,q)>0,Point({Element(pmul,p,q),-p}),?),p,1,cnt),q,1,Le)
#「更新」ボタン,
cnt=cnt+1
LRSum= CL+CR
nextC=Zip(Mod(c,2),c,LRSum)
SetValue(Cs,append(Cs,nextC))
SetValue(C,nextC)
#「リセット」ボタン
cnt=1
C={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
Cs={{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}
<振り返り>
可変な0のリストでやるにはPythonとかJSなどの言語が必要なります。
それはそれで、今度は描画のお約束に手間取るはずです。
どっちもどっちだね。
たとえば、1行に101個の点を打ちたいなら、
pythonで、ls=[0 for x in range(50)]
zeros=ls+[1]+ls
とするだけでコマンドラインにコピペ可能なデータが表示できるね。