n分岐

[ストロガッツ演習問題3.4.12]
　3分岐(超臨界ピッチフォーク分岐)は次の式だった。
　　　f(x)=dx/dt=rx-x^3
　4分岐でも同じことをすればいいんだろう。r<0で解が無く、r=0で4重解が出現し、r>0では散り散りに4つの解を持つような4次式を微分方程式とすれば良い。それらの解を±r^(1/2)、±(2r)^(1/2)とすれば、微分方程式は次の形になるはず。
　　　f(x)=(r-x^2)(2r-x^2)
　これを実装。

import numpy as np

import matplotlib.pyplot as plt

def f(x, r):
    return (r-x**2)*(2*r-x**2)

R=np.arange(-1, 5, 0.01)
x=np.arange(-5, 5, 0.01)

for r in R:
    x_d=f(x, r)

    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(r, tempx, marker='o', color='red', markersize=2)
            else:
                plt.plot(r, tempx, marker='o', color='blue', markersize=2)
        i+=1

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

plt.title('4-furcation')
plt.grid()
plt.show()

　ピッチフォーク分岐ではなくサドルノード分岐っぽくなっているが、問題文では「r<0で固定点を持たない」と書いてあるので、それでよいのだと思う。

　次は5分岐。これは4分岐にx=0を解に追加すればよい。この場合にはr<0でもx=0が固定点としてあるので、ピッチフォーク分岐になりそう。6分岐は4分岐にx=±(3r)^(1/2)を解に追加すればよい。結果としてnが偶数ならサドルノード分岐っぽく、奇数ならピッチフォーク分岐になりそう。

　そうして話は一般化へ。

(1)nが偶数：f(x)=Π [1≦k≦n/2] (kr-x^2)

(2)nが奇数：f(x)=x*Π [1≦k≦(n-1)/2] (kr-x^2)

def f(x, r, n):
    k=n%2
    n=int(n/2)
    temp=1
   
    for i in range(n):
        temp*=((i+1)*r-x**2)

    return temp*(x**k)

R=np.arange(-2, 2, 0.01)
x=np.arange(-2, 2, 0.01)
N=np.arange(5, 10, 1)

for n in N:
    for r in R:
        x_d=f(x, r, n)

        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(r, tempx, marker='o', color='red', markersize=2)
                else:
                    plt.plot(r, tempx, marker='o', color='blue', markersize=2)
            i+=1

    plt.xlabel('r')
    plt.ylabel('x')
    plt.title(str(n)+'-furcation')
    plt.grid()
    plt.show()

　nが奇数の時、関数fの最初でtempにxを入れるとなぜかエラーが吐かれる。returnする部分でxをかけたら大丈夫な理由もわかってないです。

　調子に乗ってrを増やす。

　ぐちゃぐちゃ。

　ところでnが奇数の場合に、n=4k+1とn=4k+3では少し性質が違いそう。n=4k+3ではx*=0の固定点の状態がr=0で反転するが、n=4k+1の時には反転しないという見方と、r>0の最も"外側"の状態とr<0における固定点x*=0の状態が同じになるという見方ができるのかな。r>0の固定点x*=0はどっちでも不安定であることを考えると後者の方が本質っぽそうだなと思うが、特に何か理解しているわけではないのでアテにならん。