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=ジャストゼロが得られなかった場合である。現時点では練習問題にそのような問題がみられないので問題とはならないが、そのような場合どうすればいいのだろうか。

渦巻き

　昨日のBZ反応のモデルを、もう少し大きい領域で実行した。(200*200)

　中央下部が渦を巻いているように見える。

　静止画だといまいち。渦の線を延々と追えたりすると「渦巻きだ！」となれるんだが、どうも格子の上でやっても渦巻きらしい渦巻きはつくれないような気もしてきた。釈然としない。ただし昨日は円形に線を辿れたので、それに比べれば渦に近いような気がする。
　何度も実行してみたが、似たような「渦巻いてそうな何か」は何度か見られた。

これは中央上部が渦巻きっぽい？

　なぜ領域を大きくしたら見え始めたのかがわからない。実は小さい領域でやってた時も発生していたが、拡大すぎるため渦巻きだと認識できなかった可能性がある。他の波紋までの距離が近いので渦っぽく見えないということもありそう。あとは一度に発生する波紋が多いので、単純に遭遇率が上昇しただけというのもありえる。
　おまけ。2個目のgifのこの辺が...(はぁと)

ベロウソフ・ジャボチンスキー反応のモデルで遊びたかった

　※注意
　　この記事には強い点滅を伴うgifが貼ってあります。




　毎度ながら「遊ぶ」まで至らずに終わりました。ただしつい先日学んだ関数内関数を使えたというだけで個人的には満足なのです。
　なお関数内関数を使う必要があったかどうかについては考えてはいけない。
　これを参考にした。
---------------------------------------------------------------------------------
#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def defs(k1, k2, g, n, nei):
    def BZreaction(plate):
        padded=zeropad(plate)
        output=np.zeros(padded.shape)
        i=nei
        while i<padded.shape[0]-nei:
            j=nei
            while j<padded.shape[0]-nei:
                if padded[i, j]==0:
                    b=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1]>=n)
                    a=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1]>0)-b
                    output[i, j]=a/k1+b/k2
                elif padded[i, j]>=n:
                    output[i, j]=0
                else:
                    s=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1])
                    b=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1]>=n)
                    a=np.sum(padded[i-nei:i+nei+1, j-nei:j+nei+1]>0)-b
                    output[i, j]=s/a+b+1+g
                j+=1
            i+=1
       
        return output[nei:-1*nei, nei:-1*nei]

    def zeropad(plate):
        padded=np.zeros((plate.shape[0]+2*nei, plate.shape[0]+2*nei))
        padded[nei:plate.shape[0]+nei, nei:plate.shape[1]+nei]=plate
        return padded
   
    return BZreaction, zeropad

#setting
g=10
k1=2
k2=1
n=20
l=50
plate=np.random.rand(l*l).reshape(l, l)*n
nei=1

BZreaction, zeropad=defs(k1, k2, g, n, nei)

fig=plt.figure()

ims=[]
for a in range(500):
 plate=BZreaction(plate)
 im=plt.imshow(plate, animated=True)
 ims.append([im])

ani=animation.ArtistAnimation(fig, ims, interval=50, blit=True, repeat_delay=1000)
ani.save('sample.gif', writer='imagemagick')
---------------------------------------------------------------------------------

　いい感じの波模様が出てきてくれたのはいいんだけど、合っているかはわからんね。

　3つのパラメータ(g、k1、k2)を色々と動かしてみたんだが、どうにもどれを動かすとどうなるのかさっぱり掴めなかった。

g=10, k1=5, k2=5

たまたま2つ近いところで回りだしたから維持できてるって感じ。実際同じパラメータで何回かやっても消えちゃうことが多い。

g=10, k1=1, k2=1

kが小さいからstate=0からすぐに離れるし、gも大きいからすぐに上限達している感じ。

g=5, k1=5, k2=5

クソ眩しい。

g=1, k1=1, k2=1

眩しいけどgが小さいので絵面がやわらかい印象

g=10, k1=10, k2=10

多くのセルが上限を一斉に超えるがために、みんな揃ってstate=0になってしまった感じ。

g=10, k1=5, k2=5

泡が割れるときに一瞬膜が見えるみたいなやつに似てる。いっぱい試した中で一番好き。

　初っ端で波紋が出てきたのが奇跡かってくらい波紋が出てきてくれないので悲しい。

　空き時間でちょこちょこいじって、最終的には渦巻きを作ってThe 9thを聞くのが目標。

　ところで文字サイズがそろわないのは何故なのか。

感染伝播モデル

　待ち時間がたくさんできたので、数学やらネットワーク解析やらの基礎学習が捗る。が、数学疲れるのでまったり感染伝播モデルを実装してみた。

　2次元正方格子の各交点に確率qで木が生えている。中央で火を発生させた時、前後左右の4方向に木があれば燃え移る。燃え移った先の木の前後左右に木があればまた燃え移る。周りに木がなくなれば燃え移りは止まる。
　参考書に忠実に従っているわけではないのだが、とりあえず上記設定で、燃えた木の割合の平均θ(q) が急上昇するような相転移がみられるかどうかを調べる。

------------------------------------------------------------
#-*- coding:utf-8 -*-
#argv[1]: edge
#argv[2]: number of each simulation
#argv[3]: step size of q

import sys
import numpy as np
import matplotlib.pyplot as plt

def zeropad(data):
    x, y=data.shape
    data=np.c_[np.zeros(x), data, np.zeros(x)]
    data=np.r_[np.zeros(y+2).reshape(1, -1), data, np.zeros(y+2).reshape(1, -1)]
    return data

def spread(data, x, y):
    if data[x-1, y]==1:
        data[x-1, y]=2
        data=spread(data, x-1, y)

    if data[x+1, y]==1:
        data[x+1, y]=2
        data=spread(data, x+1, y)

    if data[x, y-1]==1:
        data[x, y-1]=2
        data=spread(data, x, y-1)

    if data[x, y+1]==1:
        data[x, y+1]=2
        data=spread(data, x, y+1)

    return data
               
N=np.int(sys.argv[1])
count=np.int(sys.argv[2])
q_list=np.arange(0, 1, np.float(sys.argv[3]))

if (N%2)!=1:
    print('N must be an odd integer.')
    quit()

theta=[]
for q in q_list:
    results=[]
    for i in range(count):
        data=np.random.choice(2, N**2, p=[1-q, q]).reshape(N, N)
        data=zeropad(data)
        center=np.array(data.shape)/2
        data[center[0], center[1]]=2
        out=spread(data, N/2, N/2)
        out=(out==2)
        results.append(np.mean(out[1:N+1, 1:N+1]))

    theta.append(np.mean(results))

plt.plot(q_list, theta)
plt.show()
-------------------------------------------------------------
Input:  edge:37, count:50, step size: 0.05
Output:

 それっぽくなった。ちなみにedge数37だと1/2の確率で再帰潜り過ぎでボイコットされる。39は不可であった。

お次はこれを立法格子へ応用。やることは一緒。
------------------------------------------------------------
import sys
import numpy as np
import matplotlib.pyplot as plt

def zeropad3d(data):
    x, y, z=np.array(data.shape)
    temp=np.zeros((x+2, y+2, z+2))
    temp[1:x+1, 1:y+1, 1:z+1]=data
    return temp

def spread3d(data, x, y, z):
    if data[x-1, y, z]==1:
        data[x-1, y, z]=2
        data=spread3d(data, x-1, y, z)

    if data[x+1, y, z]==1:
        data[x+1, y, z]=2
        data=spread3d(data, x+1, y, z)

    if data[x, y-1, z]==1:
        data[x, y-1, z]=2
        data=spread3d(data, x, y-1, z)

    if data[x, y+1, z]==1:
        data[x, y+1, z]=2
        data=spread3d(data, x, y+1, z)

    if data[x, y, z-1]==1:
        data[x, y, z-1]=2
        data=spread3d(data, x, y, z-1)

    if data[x, y, z+1]==1:
        data[x, y, z+1]=2
        data=spread3d(data, x, y, z+1)

    return data
              
N=np.int(sys.argv[1])
count=np.int(sys.argv[2])
q_list=np.arange(0, 1, np.float(sys.argv[3]))

if (N%2)!=1:
    print('N must be an odd integer.')
    quit()

theta=[]
for q in q_list:
    results=[]
    for i in range(count):
        data=np.random.choice(2, N**3, p=[1-q, q]).reshape(N, N, N)
        data=zeropad3d(data)
        center=np.array(data.shape)/2
        data[center[0], center[1], center[2]]=2
        out=spread3d(data, center[0], center[1], center[2])
        out=(out==2)
        results.append(np.mean(out[1:N+1, 1:N+1, 1:N+1]))

    theta.append(np.mean(results))

plt.plot(q_list, theta)
plt.show()
------------------------------------------------------------
Input: edge:9, count:50, step size:0.05
Output:

もうedge数9の時点でドキドキが足りない。しかし3*3*3のルービックキューブを27個積んであると思うと思ったよりでかい。

で、特に意味はないが4次元格子でも作ってみた。現実世界の時間に当てはめた場合にはそらもうSFよ。
実装はほんとにちょっと書き足すだけである。
------------------------------------------------------------
import sys
import numpy as np
import matplotlib.pyplot as plt

def zeropad4d(data):
    w, x, y, z=np.array(data.shape)
    temp=np.zeros((w+2, x+2, y+2, z+2))
    temp[1:w+1, 1:x+1, 1:y+1, 1:z+1]=data
    return temp

def spread4d(data, w, x, y, z):
    if data[w-1, x, y, z]==1:
        data[w-1, x, y, z]=2
        data=spread4d(data, w-1, x, y, z)

    if data[w+1, x, y, z]==1:
        data[w+1, x, y, z]=2
        data=spread4d(data, w+1, x, y, z)

    if data[w, x-1, y, z]==1:
        data[w, x-1, y, z]=2
        data=spread4d(data, w, x-1, y, z)

    if data[w, x+1, y, z]==1:
        data[w, x+1, y, z]=2
        data=spread4d(data, w, x+1, y, z)

    if data[w, x, y-1, z]==1:
        data[w, x, y-1, z]=2
        data=spread4d(data, w, x, y-1, z)

    if data[w, x, y+1, z]==1:
        data[w, x, y+1, z]=2
        data=spread4d(data, w, x, y+1, z)

    if data[w, x, y, z-1]==1:
        data[w, x, y, z-1]=2
        data=spread4d(data, w, x, y, z-1)

    if data[w, x, y, z+1]==1:
        data[w, x, y, z+1]=2
        data=spread4d(data, w, x, y, z+1)

    return data
              
N=np.int(sys.argv[1])
count=np.int(sys.argv[2])
q_list=np.arange(0, 1, np.float(sys.argv[3]))

if (N%2)!=1:
    print('N must be an odd integer.')
    quit()

theta=[]
for q in q_list:
    results=[]
    for i in range(count):
        data=np.random.choice(2, N**4, p=[1-q, q]).reshape(N, N, N, N)
        data=zeropad4d(data)
        center=np.array(data.shape)/2
        data[center[0], center[1], center[2], center[3]]=2
        out=spread4d(data, center[0], center[1], center[2], center[3])
        out=(out==2)
        results.append(np.mean(out[1:N+1, 1:N+1, 1:N+1, 1:N+1]))

    theta.append(np.mean(results))

plt.plot(q_list, theta)
plt.show()
------------------------------------------------------------
Input: edge:5,  count:50, step size:0.05
Output:

悲しみのedge数5。
5*5*5のルービックキューブ5個並んでるの想像するとそのしょーもなさが身にしみる。

Balloon-Windkessel model

　そもそもFitzHugh-Nagumo modelがちゃんと出来ているのかどうかがよくわかっていないのだが、それで出てきたデータをBalloon-Windkessel modelにぶちこんでみる。

まず前回までで作っていたFitzHugh-Nagumo modelをそのまま実行して出てきたのがこちら。

ここではu0を使うことにして、Balloon-Windkesselモデルを当てはめる。原文ママなのでコードなしでドン。

動きがいまいちわかりにくいので、最大最小が同じになるようにして表示。

ぐちゃぐちゃや・・・。
内部の値を非表示にして入出力のみドン。

　uが上昇するのに遅れてyが上昇している様子がわかる。yの低下が2段になっているが、遠くから見ればおおまかにuを反映していると言えるのかな。

　FitzHugh-NagumoのノリでBalloon-Windkesselも人名2つをつなげたものだと思っていたが、どうやらそうではないらしい。私見ではバルーンさんっていてもあまり違和感がない。日本語にすると風船さんだしいないか、いや、このご時世だと風船さんいそう。Windkesselも(聞きなれない単語なので)いそう。Windさんもいそう。風さんもどこかにいそう。Windkesselはググると空気室らしい。空気室さんはいなさそう、いや、このご時世(以下略