gpt4 book ai didi

python - Scipy ODE 时间步长倒退

转载 作者:行者123 更新时间:2023-11-28 18:08:47 24 4
gpt4 key购买 nike

我在 Stackoverflow 上四处寻找,但找不到任何可以回答我问题的内容。

问题设置:

我正在尝试使用 scipy.integrate.ode 求解刚性 ODE 系统。我已将代码缩减为最小的工作示例:

import scipy as sp
from scipy import integrate
import matplotlib.pylab as plt
spiketrain =[0]
syn_inst = 0

def synapse(t, t0):
tau_1 = 5.3
tau_2 = 0.05
tau_rise = (tau_1 * tau_2) / (tau_1 - tau_2)
B = ((tau_2 / tau_1) ** (tau_rise / tau_1) - (tau_2 / tau_1) ** (tau_rise / tau_2)) ** -1
return B*(sp.exp(-(t - t0) / tau_1) - sp.exp(-(t - t0) / tau_2)) #the culprit

def alpha_m(v, vt):
return -0.32*(v - vt -13)/(sp.exp(-1*(v-vt-13)/4)-1)

def beta_m(v, vt):
return 0.28 * (v - vt - 40) / (sp.exp((v- vt - 40) / 5) - 1)

def alpha_h(v, vt):
return 0.128 * sp.exp(-1 * (v - vt - 17) / 18)

def beta_h(v, vt):
return 4 / (sp.exp(-1 * (v - vt - 40) / 5) + 1)

def alpha_n(v, vt):
return -0.032*(v - vt - 15)/(sp.exp(-1*(v-vt-15)/5) - 1)

def beta_n(v, vt):
return 0.5* sp.exp(-1*(v-vt-10)/40)

def inputspike(t):
if int(t) in a :
spiketrain.append(t)

def f(t,X):
V = X[0]
m = X[1]
h = X[2]
n = X[3]

inputspike(t)
g_syn = synapse(t, spiketrain[-1])
syn = 0.5* g_syn * (V - 0)
global syn_inst
syn_inst = g_syn

dydt = sp.zeros([1, len(X)])[0]
dydt[0] = - 50*m**3*h*(V-60) - 10*n**4*(V+100) - syn - 0.1*(V + 70)
dydt[1] = alpha_m(V, -45)*(1-m) - beta_m(V, -45)*m
dydt[2] = alpha_h(V, -45)*(1-h) - beta_h(V, -45)*h
dydt[3] = alpha_n(V, -45)*(1-n) - beta_n(V, -45)*n
return dydt

t_start = 0.0
t_end = 2000
dt = 0.1

num_steps = int(sp.floor((t_end - t_start) / dt) + 1)

a = sp.zeros([1,int(t_end/100)])[0]
a[0] = 500 #so the model settles
sp.random.seed(0)
for i in range(1, len(a)):
a[i] = a[i-1] + int(round(sp.random.exponential(0.1)*1000, 0))

r = integrate.ode(f).set_integrator('vode', nsteps = num_steps,
method='bdf')
X_start = [-70, 0, 1,0]
r.set_initial_value(X_start, t_start)

t = sp.zeros(num_steps)
syn = sp.zeros(num_steps)
X = sp.zeros((len(X_start),num_steps))
X[:,0] = X_start
syn[0] = 0
t[0] = t_start
k = 1

while r.successful() and k < num_steps:
r.integrate(r.t + dt)
# Store the results to plot later
t[k] = r.t
syn[k] = syn_inst
X[:,k] = r.y
k += 1

plt.plot(t,syn)
plt.show()

问题:

我发现当我实际运行代码时,求解器中的时间 t 似乎来回移动,这导致 spiketrain[-1] 大于t,并且值 syn 变得非常负,严重扰乱了我的模拟(如果运行代码,您可以在图中看到负值)。

我猜它与求解器中的可变时间步长有关,所以我想知道是否可以将时间限制为仅向前(正向)传播。

谢谢

最佳答案

求解器实际上来回移动,我认为也是因为可变时间步长。但我认为困难在于 f(t, X) 的结果不仅是 tX 的函数,而且是之前调用过这个函数,这不是一个好主意。

您的代码通过替换来工作:

inputspike(t)
g_syn = synapse(t, spiketrain[-1])

通过:

last_spike_date = np.max( a[a<t] )
g_syn = synapse(t, last_spike_date)

并通过使用 a = np.insert(a, 0, -1e4) 为“稳定时间”设置“旧事件”。这是始终定义 last_spike_date 所必需的(请参阅下面代码中的注释)。

这是您的代码的修改版本:

我修改了找到最后一个尖峰时间的方式(这次使用 Numpy 函数 searchsorted 以便可以对该函数进行矢量化)。我还修改了创建数组 a 的方式。这不是我的领域,所以我可能误解了意图。

我用了solve_ivp而不是 ode 但仍然带有 BDF solver (但是它与 Fortran 中的 ode 中的实现不同)。

import numpy as np  # rather than scipy 
import matplotlib.pylab as plt
from scipy.integrate import solve_ivp

def synapse(t, t0):
tau_1 = 5.3
tau_2 = 0.05
tau_rise = (tau_1 * tau_2) / (tau_1 - tau_2)
B = ((tau_2 / tau_1)**(tau_rise / tau_1) - (tau_2 / tau_1)**(tau_rise / tau_2)) ** -1
return B*(np.exp(-(t - t0) / tau_1) - np.exp(-(t - t0) / tau_2))

def alpha_m(v, vt):
return -0.32*(v - vt -13)/(np.exp(-1*(v-vt-13)/4)-1)

def beta_m(v, vt):
return 0.28 * (v - vt - 40) / (np.exp((v- vt - 40) / 5) - 1)

def alpha_h(v, vt):
return 0.128 * np.exp(-1 * (v - vt - 17) / 18)

def beta_h(v, vt):
return 4 / (np.exp(-1 * (v - vt - 40) / 5) + 1)

def alpha_n(v, vt):
return -0.032*(v - vt - 15)/(np.exp(-1*(v-vt-15)/5) - 1)

def beta_n(v, vt):
return 0.5* np.exp(-1*(v-vt-10)/40)

def f(t, X):
V = X[0]
m = X[1]
h = X[2]
n = X[3]

# Find the largest value in `a` before t:
last_spike_date = a[ a.searchsorted(t, side='right') - 1 ]

# Another simpler way to write this is:
# last_spike_date = np.max( a[a<t] )
# but didn't work with an array for t

g_syn = synapse(t, last_spike_date)
syn = 0.5 * g_syn * (V - 0)

dVdt = - 50*m**3*h*(V-60) - 10*n**4*(V+100) - syn - 0.1*(V + 70)
dmdt = alpha_m(V, -45)*(1-m) - beta_m(V, -45)*m
dhdt = alpha_h(V, -45)*(1-h) - beta_h(V, -45)*h
dndt = alpha_n(V, -45)*(1-n) - beta_n(V, -45)*n
return [dVdt, dmdt, dhdt, dndt]


# Define the spike events:
nbr_spike = 20
beta = 100
first_spike_date = 500

np.random.seed(0)
a = np.cumsum( np.random.exponential(beta, size=nbr_spike) ) + first_spike_date
a = np.insert(a, 0, -1e4) # set a very old spike at t=-1e4
# it is a hack in order to set a t0 for t<first_spike_date (model settle time)
# so that `synapse(t, t0)` can be called regardless of t
# synapse(t, -1e4) = 0 for t>0

# Solve:
t_start = 0.0
t_end = 2000

X_start = [-70, 0, 1,0]

sol = solve_ivp(f, [t_start, t_end], X_start, method='BDF', max_step=1, vectorized=True)
print(sol.message)

# Graph
V, m, h, n = sol.y
plt.plot(sol.t, V);
plt.xlabel('time'); plt.ylabel('V');

给出:

result for V

注意:solve_ivp 中有一个事件参数可能很有用。

关于python - Scipy ODE 时间步长倒退,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/51973039/

24 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com