何故、谷川柑菜は報われなかったのか

(追記)なぜ文字サイズがバラバラになるのかいまだに原因がわかりません。

　ストロガッツの5章で、様々なタイプのロミオとジュリエットについての恋愛問題を解いた。このついでに、よりややこしい恋愛問題について考察を試みようと思った。具体的には、3人でもつれあう恋愛を扱おうとした。
　そこで複雑な恋愛を呈する作品何かあったかな～と思いを巡らしたが、普段あまり恋愛ものを観ない・読まないために、「あの夏で待ってる」まで遡ることとなった。あの夏では6人が恋愛の連鎖をしているが、その中でも最も報われなかった柑菜(異論は認める)を取り上げることとした。
　柑菜(Kanna: K)、海人(Kaito: U)、イチカ(Ichika: I)として、それぞれの好意を次のように表すこととした。

KU: 柑菜 to 海人
UK: 海人 to 柑菜
IU: イチカ to 海人
UI: 海人 to イチカ
(カイトがUなのは、Kがカンナと被ったために「うみんちゅ」の頭文字を使った)

　それぞれの好意の特徴を次のように設定した。異論は認める。

・かんな to かいと (KU)' = IU + UI
　イチカ-海人 間の恋愛が熟すことにより燃え上がる略奪要素を有する。KUとUKの効果は、それまでに進展が無かったことを考えるとゼロなのかなと。

・かいと to かんな (UK)' = 0.5*KU - UI
　基本的に柑菜に興味はないのだが、柑菜からそれとなく好意が出ていると若干その気になってくるヘタレ。でもイチカ先輩が好きになっちゃうと柑菜への好意は低下するヘタレ。

・いちか to かいと (IU)' = -0.5*KU - UK + IU + UI
　自分から海人への好意、海人から自分への好意が増すほどどんどん好きになる積極的お姉さん。ただし柑菜には若干の同情心がある。

・かいと to いちか (UI)' = IU + UI

　イチカへの好意に柑菜は関与しないヘタレ。

　開始時点は海人がイチカを気になりだしたタイミングとし、この時柑菜はすでに海人へ好意を持っているとし、初期値はこんな感じにした。

KU0=2
UK0=0
IU0=0
UI0=1

-----------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt

KU0=2
UK0=0
IU0=0
UI0=1

KUI=np.array([KU0, UK0, IU0, UI0]).reshape(4, 1)

A=np.array([
    [0, 0, 1, 1],
    [0.5, 0, 0, -1],
    [-0.5, -1, 1, 1],
    [0, 0, 1, 1]
])

dt=0.001
Tmax=2
T=np.arange(0, Tmax, dt)

for t in T:
    x=KUI[:, -1]
   
    #Runge-Kutta method
    k1=np.dot(A, x)
    temp=x+k1*dt/2
    k2=np.dot(A, temp)
    temp=x+k2*dt/2
    k3=np.dot(A, temp)
    temp=x+k3*dt
    k4=np.dot(A, temp)
    x=x+(k1+2*k2+2*k3+k4)*dt/6
   
    KUI=np.c_[KUI, x.reshape(4, 1)]
    t+=dt

plt.plot(np.arange(0, Tmax+dt, dt), KUI[0], label='Kanna to Kaito')
plt.plot(np.arange(0, Tmax+dt, dt), KUI[1], label='Kaito to Kanna')
plt.plot(np.arange(0, Tmax+dt, dt), KUI[2], label='Ichika to Kaito')
plt.plot(np.arange(0, Tmax+dt, dt), KUI[3], label='Kaito to Ichika')
plt.legend()
plt.title('Anonatsu')
plt.show()

-----------------------------------------------------------------------

　振られるのは避けられなかった。。。

　柑菜がもし「好きになればなるほど好きになる」タイプだったとしたらどうなったのだろうか。 (KU)' =0.5*KU + 0.5*UK + IU + UI としてRedo。

　どの夏も待ってない！あるのは冬のみ！氷河期到来！

　不憫なのでもうちょっと強くして (KU)' =0.9*KU + 0.9*UK + IU + UI としてみた。

　これくらい行くと逆転できたか。

　結論：柑菜はツンデレをやめるべきだった。

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はどっちでも不安定であることを考えると後者の方が本質っぽそうだなと思うが、特に何か理解しているわけではないのでアテにならん。

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

　式が出てきた時点でコードし始めていたため、「カスプカタストロフィー面」って名前がついていることを知ったのはこの記事を上げた後であった......。
　まずピッチフォーク分岐(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*)空間の分岐平面(?)で"よじれている"部分が、不安定な固定点になっているということですな。

アトラクター

　"ストロガッツ 非線形ダイナミクスとカオス"を読み始めた。5月末までにこれの8章までを読むことを目標とする。誓います！（強気）ただしわからない部分をとばすのは可(えー

　備忘録もかねて、必要の有無にかかわらずコーディングできそうなところをjupyterにまとめていくことにした。そして先日githubを始動させたので、そこにアップロードすることでモチベーションを保っていくことにした。

https://github.com/nkknj/strogatz

　本日は2章の練習問題を解いていたのだが、その中でアトラクターの表示に手間取った話。

import numpy as np
import matplotlib.pyplot as plt

x=np.arange(-3.5, 1, 0.1)
x_d=np.exp(-x)*np.sin(x)
plt.plot(x, x_d, linestyle='-', color='b', label='x_d')

#x軸上に流れの方向へ三角形を描出
for i in range(x.shape[0]):
    if x_d[i]>0:
        plt.plot(x[i], 0, marker='>', color='dodgerblue')
    elif x_d[i]<0:
        plt.plot(x[i], 0, marker='<', color='dodgerblue')
    else:
        #流れがゼロの点
        if x_d[i+1]>0:
            #直後が正の時
            plt.plot(x[i], 0, color='red', marker='s', markersize=5)
        else:
            #直後が負の時
            plt.plot(x[i], 0, color='red', marker='o', markersize=5)

plt.legend()
plt.xlabel('x')
plt.grid()
plt.show()

　この問題の直前までは単純なxのn次式とか三角関数だったので、xを0.1おきとかpi/10おきとかで設定すればピタッとx_d=0になる点が出てきて、簡単にアトラクターだとわかった。しかし指数関数が出てくるとなかなかぴったりゼロとはいかない。今後のためにも他の方法を考えた方がよい。

　そんなわけで、2点間でゼロをまたいだ時に、その2点を結んだ線分とx軸の交点をアトラクターとしてマーカーを表示する感じにした。

import numpy as np
import matplotlib.pyplot as plt

x=np.arange(-3.5, 1, 0.1)
x_d=np.exp(-x)*np.sin(x)
plt.plot(x, x_d, linestyle='-', color='b', label='x_d')

#x軸上に流れの方向へ三角形を描出
for i in range(x.shape[0]):
    if x_d[i]>0:
        plt.plot(x[i], 0, marker='>', color='dodgerblue')
    elif x_d[i]<0:
        plt.plot(x[i], 0, marker='<', color='dodgerblue')
    else:
        #流れがゼロの点
        if x_d[i+1]>0:
            #直後が正の時
            plt.plot(x[i], 0, color='red', marker='s', markersize=5)
        else:
            #直後が負の時
            plt.plot(x[i], 0, color='red', marker='o', markersize=5)

i=0
while i+1<x.shape[0]:
    if x_d[i]*x_d[i+1]<0:
        if x_d[i]>0:
            plt.plot(x[i]+(x[i+1]-x[i])*(x_d[i]/(x_d[i]-x_d[i+1])), 0, color='red', marker='o', markersize=5)
        else:
            plt.plot(x[i]+(x[i+1]-x[i])*(-x_d[i]/(x_d[i+1]-x_d[i])), 0, color='red', marker='s', markersize=5)
   
    i+=1
   
plt.legend()
plt.xlabel('x')
plt.grid()
plt.show()

　単純にwhile文加えて、後ろから1つずつ確認する仕組み。ただしぴったりゼロになった場合の対処は、x軸に△表示するついでにお願いする形となっている。

　この方法で対処できないのは、x軸に接しているのにx_d=ジャストゼロが得られなかった場合である。現時点では練習問題にそのような問題がみられないので問題とはならないが、そのような場合どうすればいいのだろうか。