Pythonでロジスティック写像(やったのは7/16)

非線形工学という講義を取っているんだが，内容が非常に高度でテスト前でもがき苦しんでる．でも，内容自体は非常に興味深いものが多く珍しく珍しくテスト対策にやる気が起きる教科でもある．

講義も終盤に差し掛かり，非線形の花形ともいえるカオス現象についての解説に入った．その中でも，導入に用いられたカオスの一つである，ロジスティック写像が非常に興味深かった．ロジスティック方程式は，もともと生物の個体の数を数理モデルに落とし込むために導入された方程式であり，連続な場合では以下の微分方程式で表される．

$\displaystyle \frac{dx}{dt}=\alpha x(1-x)$

上式において $x$ は個体数である．個体数は環境収容の限界である $x=1$ に近づくにつれて，その増加率は落ち込んでいく．また， $\alpha$ は内的自然増加率と言われる定数である．連続時間下ではこの式は非線形微分方程式であるが，変数分離法を用いて容易に解くことができ，

$\displaystyle x=\frac{1}{1+(x_0^{-1}-1)e^{-\alpha t}}$

のようになる．この解はS字を描き，

出典:第18回ロジスティック回帰：機械学習はじめよう｜gihyo.jp … 技術評論社

のようになる．この関数はロジスティック関数や，シグモイド関数などと呼ばれ，機械学習分野でも， $x_0=0.5$ とした，点 $(0,0.5)$ で点対称な関数を，ロジスティック回帰などのアルゴリズムで用いられる．

しかし，生物の増加を考えたとき離散時間で考えた方が世代の変遷が同時に行われると考えた下では妥当性が高い．そこで，先ほどの微分方程式を差分方程式へと変換する

$x\left[k+1 \right]=\alpha x\left[k \right](1- x\left[k \right])$

この方程式は連続時間内と大きく異なり， $\alpha$ の値によって個体数の変動は非常に複雑となる． $\alpha$ が小さなうちは連続時間と同様，一定の値に収束するが， $\alpha\ge 3$ 以上となると一定値に収束せず，2又，4又．．．と別れていき， $\alpha=4$ に近づくにしたがって全くカオスな値を取るようになることが知られている．そこで離散時間におけるロジスティック写像の挙動をpythonでシミュレーションさせてみる．

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.lines as lines
import matplotlib.animation as anim
import os
from tqdm import tqdm

使うパッケージは上のヤツ．イツメンって感じ．離散時間系の力学では，差分方程式を蜘蛛の巣図法と言われる手法で図的に表現することが可能である．

クモの巣図法 - Wikipedia

pythonで表すと以下のようになる

def calc_lines(func, init, num=300):
    vertices=[]
    nn = []
    yy = []
    
    x=init
    y=0
    vertices.append([x,y])
    nn.append(0)
    yy.append(x)
    
    for n in range(1, num):
        y = func(x)
        vertices.append([x, y])
        nn.append(n)
        yy.append(y)
        x=y
        vertices.append([x, y])
    return np.array(vertices), np.array(nn), np.array(yy)

x,yを入れ替えさえすればかなり容易に実装できる．あとはシミュレーションを行えばよい．

aaa = []
yyy = []
init = 0.3
os.makedirs("temp2",exist_ok=True)
for idx,a in enumerate(tqdm(np.arange(2.95,4.001,0.001))):
    fig = plt.figure(figsize=(12,8))
    fig.suptitle("$a=${:.4f}".format(a))
    ax3 = plt.subplot2grid((2,3), (0,1), rowspan=2,colspan=2)
    ax3.set_xlim(2.95, 4.0)
    ax3.set_ylim(0.0, 1.0)
    ax1 = plt.subplot2grid((2,3), (0,0))
    ax2 = plt.subplot2grid((2,3), (1,0))
    
    logistic = lambda x: a*x*(1-x)
    tt = np.arange(0.0, 1.0, 0.01)
    ax1.plot(tt, logistic(tt))
    ax1.plot(tt, tt)
    vs, ns, ys = calc_lines(logistic, init)
    ax1.plot(vs[:,0],vs[:,1],c='r',alpha=0.4,linestyle='-.')
    ax2.plot(ns, ys, alpha=0.7,marker='o')
    aa = np.ones((100,))*a
    aaa.append(aa)
    yyy.append(ys[-100:])
    ax3.scatter(aaa,yyy,c='black',s=1)
    plt.savefig("./temp2/{:04d}.png".format(idx))
    plt.close()
    del fig, ax1,ax2,ax3