FitzHugh-Nagumo modelで遊ぶという強い意志

昨日のおさらい。
頭から書いていったコードを整理する。
------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt

def floor(vector):
    shape=vector.shape
    vector=vector.reshape(-1)
    for i in range(vector.size):
        vector[i]=np.int(vector[i])
 
    vector.reshape(shape)
    return vector

def i_input(u, i):
    con=C[:, i]
    delay=D[:, i]
    temp=np.zeros(N)+u.shape[1]-1
    temp=floor(temp-delay)
    con[temp<0]=0
    temp[temp<0]=0
    pulse=np.zeros(N)
    for j in range(N):
        pulse[j]=u[j, temp[j]]
 
    iinput=np.dot(con, pulse)
    return iinput

def d_u(u, v):
    base=g(u[:, -1], v[:, -1])
    for i in range(u.shape[0]):
        base[i]+=i_input(u, i)
 
    return base

def g(u, v):
    return tau*(v+gamma*u-u**3/3)

def d_v(u, v):
    return h(u[:, -1], v[:, -1])

def h(u, v):
    return -(u-a-b*v)/tau

a=0.7
b=0.8
c=10
tau=1
gamma=1
N=2
C=np.ones((N, N))
D=np.ones((N, N))
for i in range(N):
    C[i, i]=0
    D[i, i]=0

t=0
T=20
dt=0.005

u=np.array([5,3]).reshape(-1, 1)
v=np.array([5,3]).reshape(-1, 1)

t+=dt
while t<T:
    du=d_u(u, v)
    dv=d_v(u, v)
    temp_u=u[:, -1]+du*dt
    temp_v=v[:, -1]+dv*dt
    u=np.c_[u, temp_u]
    v=np.c_[v, temp_v]
    t+=dt

fig=plt.figure()
ax=fig.gca()
ax.plot(np.arange(0, u.shape[1], 1)*dt, u[0, :], c='g', label='u1')
ax.plot(np.arange(0, v.shape[1], 1)*dt, v[0, :], c='b', label='v1')
ax.plot(np.arange(0, u.shape[1], 1)*dt, u[1, :], c='m', label='u2')
ax.plot(np.arange(0, v.shape[1], 1)*dt, v[1, :], c='r', label='v2')
ax.legend()
plt.show()
fig2=plt.figure()
ax2=fig2.gca()
ax2.plot(u[0, :], v[0, :], c='g', label='node1')
ax2.plot(u[1, :], v[1, :], c='m', label='node2')
ax2.legend()
plt.show()
------------------------------------------------------
（なおmatplotlibの使い方がいまだによくわかっていないので、やみくもにやっていい感じで表示してくれたものをとりあえずそのまま利用しているため、ラストがきたない）

bを半分にしたら振動はするけどu-v平面での回転が逆では、という話をしていた。
なおuとvがマイナスに突入することにも疑問を抱いていたが、それは問題ない気がしてきた。というのも振動子1つの場合の図でも突入している故。考えてみれば偶数乗の項が無いので、むしろゼロを挟む感じで振動するのが自然なのかもしれない。

で、間違えてuとv逆にしてたりしない？と思い式を確認したところ、参照してたレビューの式が間違ってました。これは悲しい。

修正したものが以下。
------------------------------------------------------
def d_u(u, v):
    base=g(u[:, -1], v[:, -1])
    for i in range(u.shape[0]):
        base[i]+=i_input(u, i)　#これが加算か減算かは要確認
 
    return base

def g(u, v):
    return tau*(-v+gamma*u-u**3/3)　#vの前にマイナス

def d_v(u, v):
    return h(u[:, -1], v[:, -1])

def h(u, v):
    return -(-u-a+b*v)/tau　#uの前にマイナス、b*vをプラスへ
------------------------------------------------------
なおd_uで外部からの入力に関しては加減のどちらが正しいのかわからん。感覚的には加算なのでプラスに変更しましたがこれも確認が必要。

パラメタを
a=0.7, b=0.2, c=10, tau=1, gamma=1, N=2
C=[[0, 1], [1, 0]]
D=[[0, 1], [1, 0]]

T=20
dt=0.005
u=[5;3]
v=[5;3]

とした実行結果。

素敵や。

さて今後は
・ノードを増やす
・ノイズを加える
・C, Dに乱数を入れてみる
といったことに着手しようと思うが、この後ろに続けるとぐちゃぐちゃするのでいったんここまで。