Scipy ODE时间步长倒退 [英] Scipy ODE time steps going backward

问题描述

我已经查看了Stackoverflow，但是找不到任何可以回答我问题的东西.

I've looked around on Stackoverflow, but could not find anything that would answer my question.

问题设置:

我正在尝试使用scipy.integrate.ode解决由刚性ODE组成的系统.我将代码简化为最小的工作示例:

I am trying to solve a system of stiff ODEs using scipy.integrate.ode. I've reduced the code to the minimal working example:

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变得非常负，并且大大搞乱了我的仿真(如果运行代码，您可以在图中看到负值).

I find that when I actually run the code, time t in the solver appears to go back and forth, which results in spiketrain[-1] being greater than t, and the value syn becoming very negative and significantly messing up my simulations (you can see the negative values in the plot if the code is run).

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

I am guessing it has something to do with variable time steps in the solver, so I was wondering if it is possible to restrict time to only forward (positive) propagation.

谢谢

推荐答案

求解器实际上会来回移动，而且我认为也是因为时间步长可变.但是我认为困难来自于f(t, X)的结果不仅是t和X的函数，而且是对该函数先前的调用，这不是一个好主意.

The solver do actually go back and forth, and I think also because of the variable time stepping. But I think the difficulty comes from that the result of f(t, X) is not only a function of t and X but of the previous call made to this function, which is not a good idea.

您的代码可以通过替换以下内容来实现:

Your code works by replacing:

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(请参见下面的代码中的注释).

And by setting an "old event" for the "settle time" with a = np.insert(a, 0, -1e4). This is needed to always have a last_spike_date defined (see the comment in the code below).

这是您代码的修改版本:

Here is a modified version of your code:

我修改了最后一次尖峰发生的时间(使用这次Numpy函数搜索排序，以便可以对函数进行矢量化处理).我还修改了创建数组a的方式.这不是我的领域，所以也许我误解了意图.

I modified how the time of the last spike if found (using this time the Numpy function searchsorted so that the function can be vectorized). I also modified the way the array a is created. This is not my field, so maybe I misunderstood the intent.

我使用了 solve_ivp 代替ode，但仍带有 BDF求解器(但是，它与Fortran中的ode中的实现方式不同).

I used solve_ivp instead of ode but still with a BDF solver (However it's not the same implementation as in ode which is in Fortran).

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');

给出:

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

note: There is an events parameters in solve_ivp which could be useful.

这篇关于Scipy ODE时间步长倒退的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！