我试图用Python解决微分方程 . 在这两个系统微分方程中,如果第一个变量( v
)的值大于阈值(30),则应将其重置为另一个值(-65) . 下面我把我的代码 . 问题是第一个变量达到30后的值保持不变,不会重置为-65 . 这些方程描述了单个神经元的动力学 . 方程式取自website和PDF file .
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))
解决方案如下:
这是 Julia
的正确解决方案,这里 u1
代表 v
:
这是 Julia
代码:
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)
1 回答
推荐的解决方案
这使用事件并在每次不连续后单独集成 .
最小变化“解决方案”
看起来好像你的方法适用于
solve_ivp
.Warning 我认为在Julia和_2494657中,处理这种事情的正确方法是使用事件 . 我相信下面的方法依赖于一个实现细节,即传递给函数的状态向量是与内部状态向量相同的对象,这允许我们在适当的位置修改它 . 如果是副本,这种方法是行不通的 . 此外,在这种方法中无法保证求解器采取足够小的步骤以使得达到极限的正确点将被踩踏 . 使用事件将使其更加正确并且可以推广到其他微分方程,这些微分方程可能在不连续之前具有较低的梯度 .
结果