在 ODE 中使用事件的时间答案

【问题标题】：Using the timing of the event in the ODE在 ODE 中使用事件的时间
【发布时间】：2019-02-15 03:49:47
【问题描述】：

（这是与Scipy ODE time steps going backward相关的后续问题）

我有一个方程组，我试图用 scipy 的solve_ivp 求解。这是一个最小的工作代码：

import numpy as np
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 event(t,X):
    return X[0] + 20
event.terminal = False
event.direction = +1

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

    last_inputspike = inputspike[inputspike.searchsorted(t, side='right') - 1 ]
    last_t_event = -100 #Not sure what to put here

    g_syn_in = synapse(t, last_inputspike)
    g_syn_spike = synapse(t, last_t_event)
    syn = 0.5 * g_syn_in * (V - 0) + 0.2 * g_syn_spike * (V + 70)

    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)
inputspike = np.cumsum( np.random.exponential(beta, size=nbr_spike) ) + first_spike_date
inputspike = np.insert(inputspike, 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, events=event)
print(sol.message)

我想检测何时出现尖峰（定义为 V > 20），并让尖峰的时间通过更改 g_syn_spike 影响 ODE 中的 syn，其方式类似于随机输入影响它。

本质上，我想知道这是否可能以及如何在求解器的给定迭代中访问sol.t_events 的最后一个值？

【问题讨论】：

看起来您的 f 已经将 t 和 V 作为参数，并且您为 g_syn_spike 分配了一个具有范围内的值。你不能在 f 里面放一个“if event(t, X)...”吗？

标签： python scipy ode

【解决方案1】：

我一直在寻找一种方法来模拟不同的微分方程的连续系统中的离散事件。建模此类不连续性并不是微不足道的，而两个（最近）包装在那里可以帮助您合作：

arsimulo - https://pypi.org/project/Assimulo/

simpy - https://pypi.org/project/simpy/

（而不是SAMPY，这仅适用于离散系统）

我希望这有用，万一你已经发现了一个不同的解决方案，我想听到它的内容

【讨论】：