カスプカタストロフィー面

　式が出てきた時点でコードし始めていたため、「カスプカタストロフィー面」って名前がついていることを知ったのはこの記事を上げた後であった......。
　まずピッチフォーク分岐(dx/dt=rx-x^3)の可視化に取り組んだ。(r, x*)平面にドットを打ちまくることで線分のように見せて解決を得た。

import numpy as np
import matplotlib.pyplot as plt

def f(x, b):
    return b*np.tanh(x)-x

B=np.arange(0, 2, 0.001)

for b in B:
    x=np.arange(-5, 5, 0.1)
    x_d=f(x, b)

    i=1
    while i+1<x.shape[0]:
        if x_d[i]*x_d[i+1]<=0:
            tempx=x[i]+(x[i+1]-x[i])*x_d[i]/(x_d[i]-x_d[i+1])
            if x_d[i]>x_d[i+1]:
                plt.plot(b, tempx, marker='o', color='red', markersize=2)
            else:
                plt.plot(b, tempx, marker='o', color='blue', markersize=2)
        i+=1

plt.xlabel('r')
plt.ylabel('x*')
plt.grid()
plt.show()

　ドットが離れてしまった分岐部に若干の隙間が見えているが、これ以上どうする気もない。rを細かく刻めば埋まるんだろうが。

　次に不完全分岐(dx/dt=h+rx-x^3)の可視化をしようとした。同様に(r, h, x*)空間に大量のドットを打つ。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def f(x, r, h):
    return h+r*x-x**3

X=np.arange(-5, 5, 0.05)

R=np.arange(-10, 10, 0.05)
H=np.arange(-10, 10, 0.05)

fig=plt.figure()
ax=fig.add_subplot(111, projection='3d')

for r in R:
    for h in H:
        for x in X:
            temp=f(x, r, h)
            if np.abs(temp)<1e-12:
                ax.scatter3D(r, h, x)

plt.show()

　こうなってしまう。
　ピッチフォーク分岐の可視化では「平面に点を打って線分」としたので、不完全分岐の可視化では「空間に線分を引いて平面」としようと考え、そのまま実装。rを定めた平面上に分岐図を描くことを繰り返す。

def f(x, r, h):
    return h+r*x-x**3
x=np.arange(-2, 2, 0.01)
R=np.arange(-5, 5, 1)
H=np.arange(-10, 10, 1)
fig=plt.figure()
ax=fig.add_subplot(111, projection='3d')
for r in R:
    temph=[]
    tempx=[]
    for h in H:
        x_d=f(x, r, h)
        i=1
        while i+1<x.shape[0]:
            if x_d[i]*x_d[i+1]<0:
                temph.append(h)
                tempx.append(x[i]+(x[i+1]-x[i])*x_d[i]/(x_d[i]-x_d[i+1]))
            i+=1
    ax.plot(np.ones(len(tempx))*r, temph, tempx)
plt.grid()
plt.show()

　線を引く順番が違う！
　今はrを固定してhが小さい方から線を引いているので、r軸方向から見てピッチフォーク分岐的なことが起きるr>1でギザギザしてしまう。なのでxを昇順にソートするのと同時にhの順番も入れ替えるようなことをする。同時に解像度も上げる。

def f(x, r, h):
    return h+r*x-x**3
X=np.arange(-5, 5, 0.01)
R=np.arange(-5, 5, 0.2)
H=np.arange(-10, 10, 0.2)
fig=plt.figure()
ax=fig.add_subplot(111, projection='3d')
for r in R:
    temph=[]
    tempx=[]
    for h in H:
        x_d=f(X, r, h)
        i=0
        while i+1<x_d.shape[0]:
            if x_d[i]*x_d[i+1]<0:
                temph.append(h)
                tempx.append(X[i]+(X[i+1]-X[i])*x_d[i]/(x_d[i]-x_d[i+1]))
            i+=1
   
    i=len(tempx)-1
    while i>0:
        j=0
        while j<i:
            if tempx[j]<tempx[j+1]:
                temp=tempx[j]
                tempx[j]=tempx[j+1]
                tempx[j+1]=temp
               
                temp=temph[j]
                temph[j]=temph[j+1]
                temph[j+1]=temp
            j+=1
        i-=1
    ax.plot(np.ones(len(tempx))*r, temph, tempx, color='b')
plt.grid()
plt.show()

　いい感じ。解像度をあげすぎるとベタ塗りになってしまって形が分かりにくかったので、線分同士の間に少し隙間があるくらいのほうが見やすかった。
　最初はメッシュプロットで出そうとしたんだけど、こういう折り重なったものをメッシュプロットに入れられるのかを試すのすらめんどくさそうだったので断念した。
　これに安定/不安定の判定して色付けするのはちょっと今の俺には無理そう。安定な固定点に赤いドットを上書きする程度ならできそうだけど。


追記
r=-2, 0, 2の平面で分岐図をドット打ちで確認した

#r=-2, 0, 2に関して平面で確認
def f(x, r, h):
    return h+r*x-x**3

X=np.arange(-5, 5, 0.05)
R=np.array([-1, 0, 1])
H=np.arange(-10, 10, 0.05)

plt.subplot(131)
plt.grid()
plt.title('r=-1')
r=R[0]
for h in H:
    x=np.arange(-2, 2, 0.01)
    x_d=f(x, r, h)

    i=1
    while i+1<x.shape[0]:
        if x_d[i]*x_d[i+1]<=0:
            tempx=x[i]+(x[i+1]-x[i])*x_d[i]/(x_d[i]-x_d[i+1])
            if x_d[i]>x_d[i+1]:
                plt.plot(h, tempx, marker='o', color='red', markersize=2)
            else:
                plt.plot(h, tempx, marker='o', color='blue', markersize=2)
        i+=1

plt.subplot(132)
plt.grid()
plt.title('r=0')
r=R[1]
for h in H:
    x=np.arange(-2, 2, 0.01)
    x_d=f(x, r, h)

    i=1
    while i+1<x.shape[0]:
        if x_d[i]*x_d[i+1]<=0:
            tempx=x[i]+(x[i+1]-x[i])*x_d[i]/(x_d[i]-x_d[i+1])
            if x_d[i]>x_d[i+1]:
                plt.plot(h, tempx, marker='o', color='red', markersize=2)
            else:
                plt.plot(h, tempx, marker='o', color='blue', markersize=2)
        i+=1
       
plt.subplot(133)
plt.grid()
plt.title('r=1')
r=R[2]
for h in H:
    x=np.arange(-2, 2, 0.01)
    x_d=f(x, r, h)

    i=1
    while i+1<x.shape[0]:
        if x_d[i]*x_d[i+1]<=0:
            tempx=x[i]+(x[i+1]-x[i])*x_d[i]/(x_d[i]-x_d[i+1])
            if x_d[i]>x_d[i+1]:
                plt.plot(h, tempx, marker='o', color='red', markersize=2)
            else:
                plt.plot(h, tempx, marker='o', color='blue', markersize=2)
        i+=1
       
plt.show()

　つまりさっきの(r, h, x*)空間の分岐平面(?)で"よじれている"部分が、不安定な固定点になっているということですな。