フィッツヒュー・南雲 (FitzHugh-Nagumo) 方程式をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算。
今回は神経を通る電気信号をモデル化したフィッツヒュー・南雲方程式
dv/dt = v - v3/3 -w +I
dw/dt = ε(v + a - b*w)
をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算してみた。
結果はこんな感じ。スパイクが見える。
パラメータはここと同じにしました。
https://www.comsol.com/blogs/understand-the-dynamics-of-the-fitzhugh-nagumo-model-with-an-app/
ソースコードはこちら:
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
a=0.7
b=0.8
eps=0.08
I=1
def FNmodel(t, u): #odeintのときとt,xの並びが逆
v=u[0]
w=u[1]
v_dot = (v*(1.0-(v**2)/3)-w+I)
w_dot = eps*(v+a-b*w)
return [v_dot, w_dot]
t0=0
tmax = 200
N=5000
u0=[0.0,0.0]
solver=ode(FNmodel)
solver.set_integrator('dop853')
solver.set_initial_value(u0,t0) #なぜか関数と並びが逆
t=np.linspace(0, tmax, N)
sol= np.empty((N, 2))
sol[0] = u0
k=1
while solver.successful() and solver.t < tmax:
solver.integrate(t[k])
sol[k] = solver.y
k+= 1
# Plot
fig = plt.figure(figsize=(10,24))
ax1=fig.add_subplot(3, 1, 1)
ax2=fig.add_subplot(3, 1, 2)
ax3=fig.add_subplot(3, 1, 3)
ax1.set_xlabel("Time")
ax1.set_ylabel("v")
ax1.grid(True)
ax2.set_xlabel("Time")
ax2.set_ylabel("w")
ax2.grid(True)
ax3.set_xlabel("v")
ax3.set_ylabel("w")
ax3.grid(True)
ax1.plot(t,sol[:,0],c='Red')
ax2.plot(t,sol[:,1],c='Blue')
ax3.plot(sol[:,0],sol[:,1],c='Green')
plt.show()
最近のコメント