我正在尝试用 Python 求解微分方程。
在这两个系统微分方程中，如果第一个变量的值 (v) 超过阈值 (30)，应将其重置为另一个值 (-65)。下面我贴上我的代码。问题是第一个变量的值在达到 30 后保持不变，不会重置为 -65。这些方程描述了单个神经元的动力学。方程取自此website http://www.mjrlab.org/2014/05/08/tutorial-how-to-write-a-spiking-neural-network-simulation-from-scratch-in-python/和这个PDF file http://www.mjrlab.org/wp-content/uploads/2014/05/network_python_tutorial2013.pdf.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
from scipy.integrate import odeint
plt.close('all')
a = 0.02
b = 0.2
c = -65
d = 8
i = 0
p = [a,b,c,d,i]
def fun(u,tspan,*p):
du = [0,0]
if u[0] < 30: #Checking if the threshold has been reached
du[0] = (0.04*u[0] + 5)*u[0] + 150 - u[1] - p[4]
du[1] = p[0]*(p[1]*u[0]-u[1])
else:
u[0] = p[2] #reset to -65
u[1] = u[1] + p[3]
return du
p = tuple(p)
y0 = [0,0]
tspan = np.linspace(0,100,1000)
sol = odeint(fun, y0, tspan, args=p)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plt.plot(tspan,sol[:,0],'k',linewidth = 5)
plt.plot(tspan,sol[:,1],'r',linewidth = 5)
myleg = plt.legend(['v','u'],\
loc='upper right',prop = {'size':28,'weight':'bold'}, bbox_to_anchor=(1,0.9))
解决方案如下所示：
Here is the correct solution by Julia, here u1 represent v:
这是Julia code:
using DifferentialEquations
using Plots
a = 0.02
b = 0.2
c = -65
d = 8
i = 0
p = [a,b,c,d,i]
function fun(du,u,p,t)
if u[1] <30
du[1] = (0.04*u[1] + 5)*u[1] + 150 - u[2] - p[5]
du[2] = p[1]*(p[2]*u[1]-u[2])
else
u[1] = p[3]
u[2] = u[2] + p[4]
end
end
u0 = [0.0;0.0]
tspan = (0.0,100)
prob = ODEProblem(fun,u0,tspan,p)
tic()
sol = solve(prob,reltol = 1e-8)
toc()
plot(sol)
这使用事件并在每次不连续性之后单独集成。
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
a = 0.02
b = 0.2
c = -65
d = 8
i = 0
p = [a,b,c,d,i]
# Define event function and make it a terminal event
def event(t, u):
return u[0] - 30
event.terminal = True
# Define differential equation
def fun(t, u):
du = [(0.04*u[0] + 5)*u[0] + 150 - u[1] - p[4],
p[0]*(p[1]*u[0]-u[1])]
return du
u = [0,0]
ts = []
ys = []
t = 0
tend = 100
while True:
sol = solve_ivp(fun, (t, tend), u, events=event)
ts.append(sol.t)
ys.append(sol.y)
if sol.status == 1: # Event was hit
# New start time for integration
t = sol.t[-1]
# Reset initial state
u = sol.y[:, -1].copy()
u[0] = p[2] #reset to -65
u[1] = u[1] + p[3]
else:
break
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
# We have to stitch together the separate simulation results for plotting
ax.plot(np.concatenate(ts), np.concatenate(ys, axis=1).T)
myleg = plt.legend(['v','u'])
看来你的方法工作得很好solve_ivp.
Warning我认为在朱莉娅和solve_ivp，处理这种事情的正确方法是使用事件。我相信下面的方法依赖于一个实现细节，即传递给函数的状态向量与内部状态向量是同一对象，这允许我们就地修改它。如果是副本，这个方法就行不通了。此外，这种方法不能保证求解器采取足够小的步长，从而能够踩到达到极限的正确点。使用事件将使这一点更加正确并可推广到其他微分方程，这些微分方程在不连续性之前可能具有较低的梯度。
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import FormatStrFormatter
from scipy.integrate import solve_ivp
plt.close('all')
a = 0.02
b = 0.2
c = -65
d = 8
i = 0
p = [a,b,c,d,i]
def fun(t, u):
du = [0,0]
if u[0] < 30: #Checking if the threshold has been reached
du[0] = (0.04*u[0] + 5)*u[0] + 150 - u[1] - p[4]
du[1] = p[0]*(p[1]*u[0]-u[1])
else:
u[0] = p[2] #reset to -65
u[1] = u[1] + p[3]
return du
y0 = [0,0]
tspan = (0,100)
sol = solve_ivp(fun, tspan, y0)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plt.plot(sol.t,sol.y[0, :],'k',linewidth = 5)
plt.plot(sol.t,sol.y[1, :],'r',linewidth = 5)
myleg = plt.legend(['v','u'],loc='upper right',prop = {'size':28,'weight':'bold'}, bbox_to_anchor=(1,0.9))
Result