grbl1.1+Arduino CNCシールドV3.5+bCNCを使用中。
BluetoothモジュールおよびbCNCのPendant機能でスマホからもワイヤレス操作可能。
その他、電子工作・プログラミング、機械学習などもやっています。
MacとUbuntuを使用。

CNCマシン全般について:
国内レーザー加工機と中国製レーザー加工機の比較
中国製レーザーダイオードについて
CNCミリングマシンとCNCルーターマシンいろいろ
その他:
利用例や付加機能など:
CNCルーター関係:



*CNCマシンの制作記録は2016/04/10〜の投稿に書いてあります。


2019年12月22日日曜日

クォータニオン(四元数) / Quaternion / 回転制御(その1)

今回はクォータニオンのアルゴリズムについてです(Python3.6/Jupyter notebook)。
クォータニオンは主に3DCG内のオブジェクトの回転制御で使われています。通常のXYZ軸ごとの角度制御よりも便利らしい。
しかしながら、二次元を複素平面で表現したように、三次元空間を四元数(一つの実数と三つの虚数i、j、k)によって表現しようとしている部分も興味深い。

クォータニオンは、

Q = w + x*i + y*j + z*k    (i2=j2=k2=ijk=-1)

で表現され、wは実数で残りの三つはXYZ軸に対応した複素数のようなもの。
Pythonにおいてはijkはないので、4つに分解して、

Q = [w, x, y, z]

という配列で表現しています。


回転制御:
通常回転制御する場合は、XYZ軸を何度ずつ動かせばいいのかと考えますが、クォータニオンの場合は回転軸一つと角度一つだけで制御します。

まず青いA点が単位球面上にあり、赤い回転軸(単位ベクトル:unit vector)を45度回すとBにたどり着くという状況があります。
この時点で、

・Aのベクトル(XYZ座標)
・回転軸(unit vector)のベクトル(XYZ座標)
・回転角(45度)

はすでに与えられており、この操作をするとBのベクトル(XYZ座標値)が得られるということになります。
言い換えれば、この操作方法によってAがどこに向かうのかを知りたいのであって、AからBにいくための操作方法を知りたいということではないので、きっかけが異なるとこの部分で少し戸惑います(AとBの二点間の回転移動を逆算して単位ベクトルと回転角度を求める内容は次回)。
しかしながら、赤い回転軸があれば単位球面上を最短距離でBに辿りつくことができます。わざわざ、X軸に何度、Z軸に何度、Y軸に何度と角度を三つに分解しなくても一つの角度だけで実現できるメリットがあり回転ツールとしては便利です。ただし3点はすべて単位曲面上にあるのでスカラー(半径)は1で座標値がそのままベクトルになります。

これをクォータニオンで計算するには、Aのベクトル[x, y, z]をクォータニオン式に変換し(wを付け足すだけ、w=0)、

A = [w, x, y, z]

とし、回転軸(unit vector)はベクトルのまま、

U = [ux, uy, uz]

として、クォータニオンの回転制御の公式である

w, x, y, z = Q*A*Q-

にいれれば、到着点であるBのベクトル/座標値[x, y, z]がわかります。
ただ、この公式の中身(QとQ-)が少し複雑です。

Q  = [cos(θ/2),  ux*sin(θ/2),  uy*sin(θ/2),  uy*sin(θ/2)]
Q- = [cos(θ/2), -ux*sin(θ/2), -uy*sin(θ/2), -uy*sin(θ/2)]

上式のθに回転角、回転軸:単位ベクトルのux, uy, uz値を代入。尚、Q-はQの共役クォータニオンで式中にマイナスがついており、Q*Q-=1という関係となっています。
さらに、クォータニオンの掛け算では各XYZにijkの虚数が含まれているため、実部wと虚部ijkの項に分けて計算します。

q1*q2 = [[q1w*q2w - q1x*q2x - q1y*q2y - q1z*q2z],
         [q1w*q2x + q1x*q2w + q1y*q2z - q1z*q2y],
         [q1w*q2y - q1x*q2z + q1y*q2w + q1z*q2x],
         [q1w*q2z + q1x*q2y - q1y*q2x + q1z*q2w]]
      = [w, x, y, z]

4段ありますが、上からw, x, y, z(実部, 虚部3つ)に対応しています。掛け算の結果として、[w, x, y, z]が得られるので、このうち[x, y, z]だけを取り出せばベクトル(XYZ座標値)になります。
また、ベクトル、内積、外積を使えばもう少し式が短くなります。

q1*q2:  w = q1w * q2w - (v1・v2)
        x, y, z = q1w * v2 + q2w * v1 + (v1×v2)

1段目が実部wに対応しており、2段目の虚部の式はXYZに対応した三つの値を返すので、そのままベクトルとして使えます。尚、(v1・v2)の部分はドット積(numpy.dot(v1, v2))、また(v1× v2)の部分はクロス積(numpy.cross(v1, v2))。ちなみに、v1=[q1x, q1y, q1z]はq1=[q1w, q1x, q1y, q1z]のXYZ成分です。
掛け算は複雑そうですが、コーディングする場合は一度関数にしておいて代入するだけなので楽です。
角度は一つしか登場してこないのでアニメーションなどつくるときには便利そうです。


ベクトルの正規化:
ベクトル[x, y, z]に実部wを加えればクォータニオンの式になりますが、先ほどの計算は全て正規化されていないとだめなので元のベクトルをスカラー値で割って正規化しておきます。

A = [x, y, z]

の場合、

√(x2 + y2 + z2) = 1

にならなければいけないので(ノルム:1)、

A_normalized = [x, y, z] / √(x2 + y2 + z2)

となります。実際Pythonではリストの割り算は一気にできないので、

A_normalized = np.array([x, y, z]) / ((x**2 + y**2 + z**2))**0.5
             = np.array(A) / np.linalg.norm(A)

になります。
同様にクォータニオンの場合は、

√(w2 + x2 + y2 + z2) = 1

でなければならないので、

Q = np.array([w, x, y, z]) / ((w**2 + x**2 + y**2 + z**2))**0.5

になります。
スカラー値で割って正規化して計算し、結果のベクトルを求め、それをスカラー倍すれば元のスケールに戻ります。


計算手順:
計算手順を再度まとめると、

A = [ax, ay, az] (これから動かそうとする点)

An = A / |A|     (正規化)

U = [ux, uy, uz] (単位ベクトル:回転軸)

Un = U / |U|     (正規化が必要であればしておく)

t = θ           (回転角)

以上が予め決めておくこと。
以下の式に角度と単位ベクトル(回転軸ベクトル)を代入しておきます。

Q  = [cos(t/2),  ux*sin(t/2),  uy*sin(t/2),  uy*sin(t/2)]
Q- = [cos(t/2), -ux*sin(t/2), -uy*sin(t/2), -uy*sin(t/2)]

そして、回転の公式で求めます。

w, x, y, z = Q*A*Q-  (QにAを掛けて、その結果にQ-を掛ける)

上式クォータニオン同志の掛け算は、以下の式を利用。

q1*q2 = [[q1w*q2w - q1x*q2x - q1y*q2y - q1z*q2z],
         [q1w*q2x + q1x*q2w + q1y*q2z - q1z*q2y],
         [q1w*q2y - q1x*q2z + q1y*q2w + q1z*q2x],
         [q1w*q2z + q1x*q2y - q1y*q2x + q1z*q2w]]
      = [w, x, y, z]

結果として[w, x, y, z]が求まるので、XYZ成分だけを取り出してベクトルにします。

[x, y, z]           (回転後のベクトル:単位球面上でのXYZ座標)

正規化される前のベクトルと同等のスカラー値に戻したい場合は、元のスカラー値を掛け直します。

[x, y, z] * |A|     (スカラー倍で元のスケールにする)

任意のベクトルA(青)を任意のunit vector(赤)を回転軸として45度動かしています。結果はベクトルB(緑)になります。

手順:
・まず、initial vector Aを正規化してnormalized vector Anにします。
・normalized vector Anをunit vector(回転軸)で45度回転させます。
・回転後はnormalized vector Bnが得られます。
・normalized vector BnにAのスカラー値をかけるとinitialized vector Bになります。
・よって、AからBに移動したことになります。

以下がPython(Jupyter Notebook)コード:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
A = np.array([1, -1, 1])             # initial vector A (arbitary vector)
As = np.linalg.norm(A)               # scalar of A: |A|
An = A / As                          # normalization: should be (sum(x**2 + y**2 + z**2))**0.5 = 1.0

U  = np.array([-0.1, -0.9, 0.2])     # unit vector (not normalized)
Us = np.linalg.norm(U)               # scalar of U: |U|
UV = U / Us                          # unit vector (normalized)

# quaternion multiplication function
def QbyQ(q1, q2): # q = [w, x, y, z]
    w = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
    i = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2]
    j = q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1]
    k = q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0]
    return np.array([w, i, j, k])

# quaternion rotation function
def QPQc(P, UV, t): # P: quaternion, UV: unit vector, t: rotation angle theta
    q =  [np.cos(t/2),  UV[0] * np.sin(t/2),  UV[1] * np.sin(t/2),  UV[2] * np.sin(t/2)]
    qc = [np.cos(t/2), -UV[0] * np.sin(t/2), -UV[1] * np.sin(t/2), -UV[2] * np.sin(t/2)]
    qP = QbyQ(q, P)
    result = QbyQ(qP, qc)
    return result

QA = np.insert(An, 0, 0)    # quaternion of An: [0, x, y, z]
theta = np.deg2rad(45)      # rotation angle
QB = QPQc(QA, UV, theta)    # quaternion rotation for An with UV and theta
Bn = QB[1:]                 # vectorize the quaternion result
B = Bn * As                 # new vector as same scalar as initial vector A
上記のプロット:
def plot_base(elev=25, azim=-70):
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    ax.view_init(elev=elev, azim=azim)
    ax.set(xlim=(-1, 1), ylim=(-1 ,1), zlim=(-1, 1))
    ax.set(xlabel='X', ylabel='Y', zlabel='Z')
    ax.xaxis.pane.fill = ax.yaxis.pane.fill = ax.zaxis.pane.fill = False

    t = np.linspace(0, 2*np.pi, 128+1)
    alpha = 0.7
    ax.plot(np.cos(t), np.sin(t), [0]*len(t), linestyle=':', c='red', alpha=alpha)
    ax.plot(np.cos(t), [0]*len(t), np.sin(t), linestyle=':', c='red', alpha=alpha)
    ax.plot([0]*len(t), np.cos(t), np.sin(t), linestyle=':', c='red', alpha=alpha)
    ax.plot([-1, 1], [0, 0], [0, 0], linestyle=':', c='red', alpha=alpha)
    ax.plot([0, 0], [-1, 1], [0, 0], linestyle=':', c='red', alpha=alpha)
    ax.plot([0, 0], [0, 0], [-1, 1], linestyle=':', c='red', alpha=alpha)
    return ax

ax = plot_base()
ax.plot([0, A[0]], [0, A[1]], [0, A[2]], c='b', alpha=0.4) # vector A
ax.scatter(A[0], A[1], A[2], c='b', alpha=0.4)
ax.text(A[0], A[1], A[2], s='  initial vector A', va='top')
ax.plot([0, An[0]], [0, An[1]], [0, An[2]], c='b')         # vector A normalized
ax.scatter(An[0], An[1], An[2], c='b')
ax.text(An[0], An[1], An[2], s='  normalized vector An', va='top')

ax.plot([0, UV[0]], [0, UV[1]], [0, UV[2]], c='r')         # unit vector
ax.scatter(UV[0], UV[1], UV[2], c='r')
ax.text(UV[0], UV[1], UV[2], s='  unit vector', va='top')

ax.plot([0, B[0]], [0, B[1]], [0, B[2]], c='g', alpha=0.4) # vector B
ax.text(B[0], B[1], B[2], s='initialized vector B  ', ha='right', va='top')
ax.scatter(B[0], B[1], B[2], c='g', alpha=0.4)
ax.plot([0, Bn[0]], [0, Bn[1]], [0, Bn[2]], c='g')         # vector B normalized
ax.scatter(Bn[0], Bn[1], Bn[2], c='g')
ax.text(Bn[0], Bn[1], Bn[2], s='normalized vector Bn  ', ha='right', va='top')

Arc = [QPQc(QA, UV, t) for t in np.linspace(0, theta, 33)] # arc between A and B
arw, arx, ary, arz = np.array(Arc).T
ax.plot(arx, ary, arz, linestyle='--', c='b', alpha=0.6)
ax.plot(arx*As, ary*As, arz*As, linestyle='--', c='b', alpha=0.4)

まとめ:
クォータニオンの回転制御の公式に従えば簡単に操作できることがわかりました。それよりも3Dグラフに補助線を引いたりするコーディングのほうが大変でした。回転操作の計算内訳は少々複雑ですが、一旦ファンクション化してしまえば、あとは代入するだけなので手順はシンプルです。回転制御だけでなく、四元数そのものも気になります。
次回は、少し疑問に思っていた二点間の回転移動を逆算してみようと思います。

関連:
クォータニオン / 回転制御(その2):二点間から単位ベクトルと回転角度を求める
クォータニオン / Quaternion:立方体の回転操作(その3)

四元数
四元数
posted with amazlet at 19.12.22
今野 紀雄
森北出版
売り上げランキング: 67,351

2019年12月15日日曜日

ロジスティック・マップ/分岐図 / Logistic Map/Bifurcation Diagram

今回はロジスティック・マップのアルゴリズムについてです(Python3.6/Jupyter Notebook)。
方程式はシンプルで、

xn+1 = a * xn * (1 -  xn)

前回のローレンツアトラクタのように、係数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がグラフとともに表示されます。
*Web上ではグラフが静止画となってしまうため単なる画像を掲載しています。Jupyter Notebook上、もしくはbinder上でコードをランさせてください。このコードGistURLはこちら


これはa=3.84のとき、3箇所の高さに点がならんでいます。
このように、ある一定の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あたりでまた二分して4つに経路が増え、それ以降は複雑。


ループ回数を上げて分岐図を点描表示する:
ということから、方程式のループ回数を上げ、ループ回数の多いグラフをある一定数点描表示させることにします。

コード:
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万くらいにあげれば直線が曲線になるのかもしれません。


まとめ:
ロジスティック・マップ/分岐図もカオス理論のはじめによく登場してきますが、前回試したローレンツアトラクタのように直接数式で描かれる形ではないので要注意。それと、線で描かれたグラフではなく、点の分布図であり、その中にパターンを見つけるという感じで、やや抽象的。
グラフの右のほうはカオスに見えまずが、分岐図と呼ばれているだけあって、無数に分岐した点群が重なり合って周期性のある模様が浮かび上がっています。形やパターンというのは連続した線で直接描かれているばかりではなく、このように無数の点が重なることで潜在的な形を浮かび上がらせることもありうるという点で興味深い。

人気の投稿