方程式はシンプルで、
xn+1 = a * xn * (1 - xn)
前回のローレンツアトラクタのように、係数aと初期値xを与えて何回かループさせれば軌道を描くと思いましたが、描画するには少し工夫が必要だったのでメモ。
前回のローレンツアトラクタのように、係数aと初期値xを与えて何回かループさせれば軌道を描くと思いましたが、描画するには少し工夫が必要だったのでメモ。
横軸を係数aの変化（0.0〜4.0）、縦軸をxの値（0.0〜1.0）とし、左からグラフが徐々に二分されていき、右のほうではカオスになっているようです。二分岐の繰り返しによって、拡大すると同じようなパターンが見えるフラクタルの性質も持ち合わせています。しかしながらカオス特有の初期鋭敏性はないようです。
係数aと方程式のループ数の関係：
そのままこの方程式に値を入れて描く前に、この方程式から得られる値の特性を見ておきます。
まずは、係数aを固定し（範囲：0.0〜4.0）、初期値x（範囲：0.0〜1.0）を代入し何回かループさせてみます。
コード：
a = 3.07 x = 0.01 X = [] iterations = list(range(100)) for _ in iterations: x = a * x * (1 - x) X.append(x) plt.figure(figsize=(8, 4)) plt.axis([0, 100, 0, 1]) plt.xlabel('iterations') plt.ylabel('X') plt.scatter(iterations, X, c='r', s=3) plt.plot(iterations, X, lw=0.5)
結果のグラフ：
a=3.07、x=0.01で100回方程式をループさせると、このようなグラフになります。横軸はループさせた回数、縦軸はループ回数に応じたxの値。
これから分かることは、立ち上がりの数ループ（5ループ程度）を除けば、あるaの値のときには、何回ループさせても値は2箇所に集約されるということ。
また、a=3.53にすると、以下のように4箇所のx値を行き来します。
同様にいくつかaの値を変えてみると、
これはa=3.63、この場合左端の立ち上がりを除けば6箇所の高さに点がならんでいます。
ipywidgets：
パラメータaを動かすと挙動が変わるため、ipywidgetsのスライダーをつけてインタラクティブに操作できるようにしてみました（Jupyter Notebook用）。
from ipywidgets import interact @interact(a=(0, 4, 0.01)) def logistic(a): N = 100 x = 0.01 X = [] iterations = list(range(N)) for i in iterations: x = a * x * (1 - x) X.append(x) plt.figure(figsize=(8, 4)) plt.axis([0, N, 0, 1]) plt.xlabel('iterations') plt.ylabel('X') plt.scatter(iterations, X, c='r', s=3, label=a) plt.plot(iterations, X, lw=0.5)
ipywidgetsを使うと、以下のようにスライダーなどのUIがグラフとともに表示されます。
これはa=3.84のとき、３箇所の高さに点がならんでいます。
このように、ある一定のaの値のときには、x値に周期性や規則性が現れるようです。またこのような規則性を得るには方程式をある程度の回数ループさせる必要があるということです。
ロジスティック・マップ/分岐図の下準備：
aを0.0〜4.0の間で変化させるとxの値に規則性が生まれることから、次に横軸をa、縦軸にxとしたグラフを描いてみます。x値の規則性が得られるまでにある程度のループ数も必要なことから、ループ数についても調べてみます。
コード：
import numpy as np import matplotlib.pyplot as plt A_steps = 1000 A = np.arange(0.0, 4.0, 4.0/A_steps) x = 0.01 x = A * x * (1 - x) fig = plt.figure(figsize=(12, 6)) ax = fig.add_subplot(111) ax.axis([0, 1, 0, 4]) ax.set_xlabel('A', fontsize=14) ax.set_ylabel('X', fontsize=14) ax.plot(A, x)
グラフ：
横軸aは範囲0.0〜4.0、分解能を1000にしておきます。xの初期値は0.01（0にならない程度で）。
方程式のループ回数が1回だと、このようにわずかながら右上がりの直線になります。あまり変化はありません。
方程式のループ回数を増やす：
x値に規則性がでてくるには、ある程度ループさせる必要があるのでループ回数を増やしてみます。
コード：
A_steps = 1000 A = np.arange(0.0, 4.0, 4.0/A_steps) x = 0.01 X_iter = 10 for i in range(X_iter): x = A * x * (1 - x) fig = plt.figure(figsize=(12, 6)) ax = fig.add_subplot(111) ax.set_xlabel('A', fontsize=14) ax.set_ylabel('X', fontsize=14) ax.plot(A, x)
方程式のループ数を10回に増やしてみました。
グラフ：
aの値（0.0〜4.0、1000ステップ、刻み幅0.04）に応じたxの値（10回ループ後の値）です。
1回ループと比べると、かなり複雑な曲線になっています。右に行くほど上下の振幅が激しい。
方程式のループ回数ごとのグラフの表示：
次にループ数を20回まで増やし、1ループから20ループまでの20種類の線を表示してみます。
コード：
A_steps = 1000 A = np.arange(0.0, 4.0, 4.0/A_steps) x = 0.01 X_iter = 20 fig = plt.figure(figsize=(12, 6)) ax = fig.add_subplot(111) ax.set_xlabel('A', fontsize=14) ax.set_ylabel('X', fontsize=14) for i in range(X_iter): x = A * x * (1 - x) ax.plot(A, x, lw=1, label=i) plt.legend()
グラフ：
左に見える番号がグラフにおける各ループ数（イテレーションの番号）です。グラフ下の青い線はループ1回のときの線で、ループ数が増えるごとに右肩上がりとなり、さらには折れ曲がった複雑な線になっていきます。
ループ回数が多いグラフは、横軸Aの1.0から2.7くらいまで右上がりの曲線に収束しているようにも見えます。
方程式のループ回数の多いグラフだけを表示する：
ループ回数の少ないほうのグラフはまだ規則性が安定していないため、ループ回数の多いグラフを5本だけを表示させます。つまり、ループ回数16、17、18、19、20回の5種類のグラフの表示で、ループ回数1から15回までのグラフは非表示。
コード：
A_steps = 1000 A = np.arange(0.0, 4.0, 4.0/A_steps) x = 0.01 X_iter = 20 X_lines = 5 fig = plt.figure(figsize=(12, 6)) ax = fig.add_subplot(111) ax.set_xlabel('A', fontsize=14) ax.set_ylabel('X', fontsize=14) for i in range(X_iter): x = A * x * (1 - x) if i >= X_iter - X_lines: ax.plot(A, x, lw=1, label=i) plt.legend()
グラフ：
このようにループ回数が少ないグラフを除けば、最初に掲載したロジスティックマップ/分岐図に近いグラフとなります。
横軸の0.0から1.0まではほぼ水平、1.0から2.7あたりまでは右上がりの一本、そして二分されて、3.26付近で上下が入れ替わり、3.5あたりでまた二分して４つに経路が増え、それ以降は複雑。
ループ回数を上げて分岐図を点描表示する：
ということから、方程式のループ回数を上げ、ループ回数の多いグラフをある一定数点描表示させることにします。
コード：
A_steps = 1000 A_min = 0 A_max = 4.0 A = np.arange(A_min, A_max, (A_max-A_min)/A_steps) x = 0.01 X_iter = 400 X_lines = 200 X_min = 0 X_max = 1 fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111) ax.axis([A_min, A_max, X_min, X_max]) ax.set_xlabel('A', fontsize=14) ax.set_ylabel('X', fontsize=14) for i in range(X_iter): x = A * x * (1 - x) if i >= X_iter - X_lines: ax.scatter(A, x, color='black', marker='o', s=0.1, alpha=0.6)
ループ数を400回にし、そのうち後半の200回分を表示（点描）させることにしました。
これでロジスティック・マップ/分岐図のアルゴリズムはほぼ完成しました。フラクタルな性質も持ち合わせているため、部分的にグラフを拡大できるように改良しようと思います。
直感的にこのグラフの仕組みがわかりにくいので、方程式で得られる連続する一本の線だけでグラフを描くと以下のようになります。
これは400回ループ後のxの値です。
横軸Aはnp.array()なので、一回のループで一本のグラフができあがります。それを400回ループするので400本のグラフができあがります。一つ前のグラフ（分岐図）では400回ループの内201回から400回ループさせた結果（200種類のグラフ）を重ね合わせて点描していることになります。おそらく、規則性が安定するのは最初の数十ループ程度だと思うので、400回の内後半300回を表示させても問題ないかもしれません。
試しに、400回ループの後半200回分の線だけでなく全て（1〜400回分）を表示させると以下のようになります。
規則性がでていない立ち上がりの状態も含むため余計な線が見えます。特に横軸Aの値1.0から3.0までの右上がりの線は、複数の線が収束して重なっており太くなっているのがわかります。特に分岐点付近では線の太さが目立っています。このことから、線を細くシャープに見せるには、ループ前半の立ち上がりの線をできるだけ含めない方がいいのかもしれません。
通常は、横軸の変化に応じて縦軸に結果が反映されますが、横軸Aによって縦軸xが決まるというだけでなく、方程式を何回ループさせるかでも縦軸xの値が変わるのでわかりづらい。二次元グラフなのに、実はAだけでなくループ数Nという変数もあるので、忠実に計算過程を表すには三次元グラフになるはず（この辺は次回への課題）。
部分拡大可能にする（フラクタル確認用）：
ipywidgetsを使ってインタラクティブに部分拡大できるようにしてみます。またグラフの画像を撮影するボタンもつけておきます。ipywidgetsには、interact、interactive、interactive_outputがありますが、細かな設定が可能なinteractive_outputを使うことにしました。
コード：
import numpy as np import matplotlib.pyplot as plt from ipywidgets import interactive_output, VBox, HBox, Label, Button from ipywidgets import FloatSlider, FloatRangeSlider, Layout from IPython.display import display import datetime %matplotlib inline def logistic_bifurcation(A_range, X_range, Dot_size): A_steps = 1000 A = np.arange(A_range[0], A_range[1], (A_range[1]-A_range[0])/A_steps) x = 0.01 X_iter = 400 X_lines = 200 X_min = 0.0 X_max = 1.0 A = np.arange(A_range[0], A_range[1], (A_range[1]-A_range[0])/A_steps) fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111) ax.set_xlim((A_range[0], A_range[1])) ax.set_ylim((X_range[0], X_range[1])) plt.xlabel('A: step size at ' + str(A_steps), fontsize=14) plt.ylabel('X: every ' + str(X_iter) + ' iterations', fontsize=14) for i in range(X_iter): x = A * x * (1 - x) if i >= X_iter - X_lines: ax.scatter(A, x, color='black', marker='o', s=Dot_size, alpha=0.6) A_range = FloatRangeSlider(value=[2.8, 4.0], min=0.0, max=4.0, step=0.01, layout=Layout(width='60%')) X_range = FloatRangeSlider(value=[0.0, 1.0], min=0.0, max=1.0, step=0.01, layout=Layout(width='60%')) Dot_size = FloatSlider(value=0.2, min=0.01, max=1.0, step=0.01, layout=Layout(width='40%')) BT = Button(value=False, description="save image") S1 = HBox([Label('RANGE A:'), A_range]) S2 = HBox([Label('RANGE X:'), X_range]) S3 = HBox([Label('DOT SIZE:'), Dot_size, BT]) UI = VBox([S1, S2, S3]) W = interactive_output(logistic_bifurcation, {'A_range':A_range, 'X_range':X_range, 'Dot_size':Dot_size}) display(UI, W) def button_pressed(b): logistic_bifurcation(A_range.value, X_range.value, Dot_size.value) d = str(datetime.datetime.now()).split('.')[0] plt.savefig('logistic' + d + '.png') plt.close() print('image saved at ' + d) BT.on_click(button_pressed)
ipywidgetsのButtonはvalueを持たないため少し扱いづらいのですが、ボタン用にbutton_pressed()ファンクションをつくり、ボタンが押されたら一度logistic_bifurcation()を現在のvalueを引数として発動させます。そして画像保存をし、再度もう一つ描画されないようにplt.close()を追加しておきます。
UIを含めて外観は以下。
横軸Aの範囲設定スライダー、縦軸Xの範囲設定スライダー、点描の点の大きさ調整スライダー、画像保存ボタンがついています。
＊問題点としては、スライダーを動かしても画像更新に時間がかかるので多少待たないといけない。ipywidgetsのスライダーの刻み幅が0.01までなので、それほど拡大できない（直接数値入力で範囲設定したほうがいいかもしれない）。スライダーが動かしにくい時はキーボードの矢印キーで微調整するといい。ちなみに、print()で変数値などを出力させようとすると、更新ごとに画面がチラつくのでつかわないほうがいい。
A：3.52-3.62、X：0.80-0.90の範囲を拡大してみると、分岐したあとにさらに分岐点が見えます。さらに拡大しても同じパターンの繰り返し（フラクタル）になっているようです。
また縦に空隙がありますが、この間隔も一定の割合で配置されておりフラクタルになっているようです。
プログラム上、横軸は常に分解能1000あるので、拡大しても画像が荒れるということはありませんが、0.01刻みのスライダーだと操作に限界があるので、この部分は数値入力に改良したほうがよさそうです。
数値手入力で拡大：
ipywidgetsをつかわず、そのまま変数に数値を手入力して拡大してみました。
範囲はA:3.5695-3.5703、X:0.8921-0.8927です。横：0.0008、縦：0.0005の画面となり、元の画像（A:2.8-4.0）の1500倍拡大したことになります。元画像中央上にみえる分岐点（左から3段目の分岐点付近）。
描画している線の数が足りていないためか、点分布というよりは曲線が見えてしまっています。400ループのうちの後半200本の線を描画しているので、もっとループ数をあげ、表示する線の数も増やした方がよさそうです。しかし、これはこれで細部の仕組みがわかって面白い。
Aの範囲を-2.0から4.0で表示：
これまでAの値は0.0から4.0でしたが、Aが負の場合-2.0までであれば表示できるようです。ちなみにXは今まで通り0.0から1.0の範囲です。
Aが-2.0から4.0で、Xが0.0から1.0なので、縦横比1:6で表示してみました。
グラフ左側では、Xの値が上下にまだありそうなので、A=-2.0のときの最小値と最大値を求めてみます。
コード：
a = -2.0 x = 0.001 x_min = np.inf x_max = -np.inf for _ in range(1000): x = a * x * (1 - x) x_min = min(x_min, x) x_max = max(x_max, x) print(x_min, x_max)
とりあえず、1000回ループ後のxの値を求めてみました。
結果は、
最小値：-0.49999979181325677
最大値：1.4999991672531139
なので、Xは-0.5から1.5の間にx存在していそうです。
ということから、Xの範囲を拡張して再度グラフを表示すると、
こんな感じで全体像が見えます。
グラフの左端、つまりA<-2.0になると、Xの最大値が無限になってしまうので、これが限界。
同様にグラフの右端、A>4.0になると、Xの最小値がマイナス無限となり、これも限界。
つまり、A:-2.0〜4.0、X:-0.5〜1.5という範囲となります。
単純なアルゴリズム：
横軸Aの変化と連続ループ計算による方法でコーディングもしてみました。先ほどの方法では、Aの刻み幅ごとに400ループさせたx値をプロットしていましたが、今回は刻み幅ごとではなく、全体を通して連続でループさせます。つまり、Aが上限4.0に近づくほどループ数が増えていくことになります。
コード：
import numpy as np import matplotlib.pyplot as plt A_steps = 200000 A_min = 1.0 A_max = 4.0 A = np.arange(A_min, A_max, (A_max-A_min)/A_steps) x = 0.01 X_min = 0.0 X_max = 1.0 X_iter = 400 X_lines = 200 fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111) ax.axis([A_min, A_max, X_min, X_max]) ax.set_xlabel('A', fontsize=14) ax.set_ylabel('X', fontsize=14) X = [] for a in A: x = a * x * (1 - x) X.append(x) ax.scatter(A, X, color='black', marker='o', s=0.1, alpha=0.6)
横軸Aの分解能を200000（20万）まで増やし、全体を通して200000回方程式を繰り返してx値を更新していきます。ただし、Aの初期値がA<0.4817...だとずっと0になってしまうので、とりあえず1.0を入れておきます。おそらく、最初のほうはループ数が足りないためか、Aが0.0から1.0までは安定しない形になってしまいます。Aの分解能200000に応じたx値が得られるので、点の個数も200000となります。
グラフ：
結果は似ていますが、分岐点においてやや縦に直線が現れます。
分解能を50000に減らしてみると、以下のように分岐点の縦線がやや長くなります。点の数も50000に減ったので全体的に色の薄い画像になります。
このことから、分解能を増やせば分岐点の直線は短くなって曲線を帯びて、ループ数を使ったアルゴリズムの結果に近づいていくのではないでしょうか。
ちなみに分解能を2000000（200万）にして試してみました。CPU計算で17秒。
200万回でも分岐点に縦の直線がまだ少し現れています。1000万くらいにあげれば直線が曲線になるのかもしれません。
まとめ：
ロジスティック・マップ/分岐図もカオス理論のはじめによく登場してきますが、前回試したローレンツアトラクタのように直接数式で描かれる形ではないので要注意。それと、線で描かれたグラフではなく、点の分布図であり、その中にパターンを見つけるという感じで、やや抽象的。
グラフの右のほうはカオスに見えまずが、分岐図と呼ばれているだけあって、無数に分岐した点群が重なり合って周期性のある模様が浮かび上がっています。形やパターンというのは連続した線で直接描かれているばかりではなく、このように無数の点が重なることで潜在的な形を浮かび上がらせることもありうるという点で興味深い。
