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

クモの巣図法（2011/2/20）のつづき。

今まで使ってきたロジスティック写像は、だった。このグラフのクモの巣図法とフィードバックループを回したときの関数出力値の変化は下のグラフのようになっており、初期値0.2からはじめると最終的に0.5に収束していた。

上の図を描くPythonプログラム。matplotlibのsubplotを使うと2つのグラフを表示できて便利！

#coding:utf-8
import numpy as np
from pylab import *
import matplotlib.lines as lines

"""ロジスティック写像族"""

def drawGraph(func):
    """関数funcを描画"""
    xList = []
    yList = []
    for x in np.arange(-0.1, 1.0, 0.01):
        xList.append(x)
        yList.append(func(x))
    plot(xList, yList)

def drawCobweb(func, initial):
    """funcの軌道に対する出力の遷移とクモの巣図法を描画"""
    vertices = []
    nList = []
    yList = []

    # 初期座標
    x = initial
    y = 0
    vertices.append([x, y])
    nList.append(0)
    yList.append(x)

    for n in range(1, 100):
        # 垂直方向
        y = func(x)
        vertices.append([x, y])
        nList.append(n)
        yList.append(y)
        print n, y
        # 水平方向
        x = y
        vertices.append([x, y])
    vertices = np.array(vertices)
    plot(vertices[:,0], vertices[:,1], '--')

    # 2つめのグラフを描画
    subplot(212)
    plot(nList, yList)

if __name__ == "__main__":
    # 1つめのグラフを描画
    subplot(211)

    # y = xを描画
    drawGraph(lambda x: x)

    # y = ax(1-x)を描画
    a = 3.3
    drawGraph(lambda x: a * x * (1 - x))

    # 初期値を変えてクモの巣図法を描画
    initial = 0.9
    drawCobweb(lambda x: a * x * (1 - x), initial)

    subplot(211)
    axis([0.0, 1.0, 0.0, 1.0])
    xlabel("x")
    ylabel("y")

    subplot(212)
    axis([0, 100, 0, 1.0])
    xlabel("n")
    ylabel("y")

    show()

これを拡張してロジスティック写像の族 を考えるとすごく面白いことが起きる。まず、a = 3.3として初期値を0.2、0.5、0.9と変えてクモの巣図法を描いてみた（上のaとinitialを変えるだけ）。

初期値によって収束の速さは違うけど最終的に {0.4794, 0.8236} の2つの点の周期的な軌道をとることがわかる。次に、a = 3.5としてみた（初期値は0.2で固定）。

今度は、{0.3828, 0.8269, 0.5009, 0.8750} の4つの点の周期的な軌道をとることがわかる。次は、a = 4.0としてみた。

画面全体を覆いつくすような不規則な軌道をとる。これはすごい！

最後に、パラメータaを横軸にとってそのときの周期的な軌道が取る値（さっきの {0.3828, 0.8269, 0.5009, 0.8750} など）を縦軸にとったグラフを描いてみる。これが有名なロジスティック写像の分岐図。Pythonを使ってちょっと書いてみよう。カオス 第1巻 - 力学系入門によると下のような手順で描くとのこと。

1. a = 1を初めのパラメータ値として選ぶ
2. 初期値xを[0,1]の中からランダムに選ぶ
3. の下でのxの軌道を計算する
4. 初めの100回の反復を無視し、101回目の反復から軌道をプロットする
5. aを決められた増分だけ増やして、(1)に戻る

(4) で101回目からプロットするのは、軌道が収束するのを待っているのだろう。Pythonで描いてみた（Pythonの画像ライブラリPILが必要）。

# coding:utf-8
import random
from PIL import Image

# キャンバスを作成
width = 500
height = 500
canvas = Image.new("RGB", (width+1, height+1), (0, 0, 0))

# aが1.0から4.0の範囲を描画
a1 = 1.0
a2 = 4.0

for i in range(width):
    # aの増分はキャンパスサイズによって決まる
    a = a1 + (a2 - a1) * float(i) / width
    # [0, 1]からランダムに初期値を決める
    x = random.random()
    for j in range(1000):
        x = a * x * (1 - x)
        # 101回目からの軌道をプロット
        if j > 100:
            # キャンバスの座標は上下逆なので注意
            canvas.putpixel((i, height - int(x * height)), (255, 255, 255))

canvas.save("bifurcation.png")

結果は、

一番左端がa = 1.0、右端がa = 4.0を意味する。図1.7 (a)にあるようにa = [3.4, 4.0]の部分に拡大してみる（上のプログラムのa1とa2を変えるだけ）と

真ん中よりちょっと右寄りのところに大きな空白地帯があることに気づく。その部分（a = [3.82, 3.86]）をさらに拡大してみると・・・

でたらめな周期を取っていた軌道が3周期の軌道に落ち着く箇所がある。この部分は「カオスの窓」と呼ばれるそうだ。何かぴったりの名前だね。あと上の2枚の画像を比べるとわかるように拡大したカオスの窓の部分にも1枚目に示した全体と同じように倍分岐していく様子が見え、小さなカオスの窓もある。これは、部分と全体が自己相似するフラクタルの例になっている。この事実を初めて発見した人はさぞかし興奮しただろうなー。うらやましい。

分岐図については、ストレンジアトラクタ（2006/3/25）でも描いたことがある。前はJavaだったけど、今回はPythonを使ってみました。

参考リンク

• Recipe 577332: Bifurcation diagram of logistic map (Python)
• Bifurcation diagram - Wikipedia