如何在python中求解一组包含dirac delta函数（或逐步递增）的ode？

Question

我正在尝试用python可视化振荡器系统。只要其中一个变量没有积累（使用dirac delta函数），系统就可以工作。但是，当我试图包含它时，系统似乎没有任何变化。在下面你可以找到我所经历的过程和我的一些问题。在

我用python编写了一些代码，能够正确地可视化常规的振荡器行为。然而，当在其中一个变量中包含逐步递增时，系统不会对此作出反应。下面，您可以找到包含增量的版本。代码运行，但给出的输出就好像不存在增量一样。在这段代码下，你可以找到我在获得正确解决方案方面的问题。在

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

a4 = 5.68
aAPD = 9.09
deltaZ = 0.41
gamma = 1.33
negV = 1
n=500

def f(V):
    if V<=-60:
        return V/20 + 4
    elif V>=20:
        return V/20 - 1
    else:
        return -V/80 + 1/4

def alpha(V):
    if V<0:
        return a4
    else:
        return -aAPD

def g(q):
    return q/(q+1)

def step(V):
    global negV
    if V>0 and negV == 1:
        negV = 0
        return deltaZ
    elif V<=0 and negV == 0:
        negV = 1
        return 0
    else:
        return 0

# function that returns dz/dt
def model(z,t):
    V = z[0]
    y = z[1]
    q = z[2]
    dVdt = 25000*(y - f(V))
    dydt = alpha(V) - a4*g(q)
    dqdt = -gamma*g(q) + step(V)
    dzdt = [dVdt,dydt,dqdt]
    return dzdt

# initial condition
z0 = [0,0,0]

# time points
t = np.linspace(0,0.5,n)

# solve ODE
z = odeint(model,z0,t)

# plot results
plt.plot(t,z[:,0],label='Voltage')
plt.plot(t,z[:,1],label='y')
plt.plot(t,z[:,2],label='Z')
plt.ylabel('response')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

问题出在函数step(V)。振荡器在某一时刻与值V=0相交。此时，q的值应该以deltaZ递增。在

在python的odesolver中有几点我不理解或不知道如何访问正确的值。因此，这里有一个简短的问题清单：

• 为了找到V与0相交的点，我使用了一个全局值negV。但是，如果我能够包含一个沿着if V[i]>0 and V[i-1]<0: return deltaZ的语句，那就更好了。这些信息应该可以在某处找到，因为你最后画出了V的值，但是我不知道如何才能访问V之前的值
• q值的增加应该是瞬时的。因此，一种解决方案可能是先求解dqdt，然后将q值增量为deltaZ。但我不明白这是怎么做到的。有没有办法直接取V，y，q的值？稍后，我还想手动更改V的值，因此这将非常有用。在
• 另一个解决方法可能是在dqdt中保持deltaZ值，但是我应该知道微分方程求解的时间步长，这样我就可以用这个时间步长除以deltaZ=0.41，这样q的增量就是deltaZ=0.41。但是，当查看odeint信息时，我发现时间步长是可变的，所以这可能不是一个好的解决方案。。。在

我最好能解决这三个问题，但如果你能给出其中一个的答案，我也会得到很多帮助。在

Answer 1

我强烈建议在scipy中使用solve_ivp函数，使用一个事件来检测V=0交叉。见here。在

如果使用终端事件，则解决方案将精确地在V=0处停止。然后，您可以使用您的即时更改更新ODE集成器外部的解决方案，并将其用作第二次调用solve_ivp以完成集成的初始条件。我假设V=0的交叉只发生一次，但是您可以构建逻辑来处理事件的多个交叉。在