如何在微分方程（SciPy）中使用if语句？ [英] How to use if statement in a differential equation (SciPy)?
问题描述
我正在尝试使用Python解微分方程。
在这两个系统的微分方程中，如果第一个变量（ v ）的值大于阈值（30），则应将其重置为另一个值（-65 ）。下面我放我的代码。问题在于，第一个变量的值达到30后保持不变，并且不会重置为-65。这些方程式描述了单个神经元的动力学。方程式取自
这是朱莉娅的正确解决方案，此处 u1 表示 v ：
这是 Julia 代码：
使用微分方程
使用图表
a = 0.02
b = 0.2
c = -65
d = 8
i = 0
p = [a，b，c，d，i]
函数fun（du，u，p，t）
如果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]）
否则
u [1] = p [3]
u [2] = u [2] + p [4]
结束
结束
u0 = [0.0; 0.0]
tspan =（0.0,100）
概率= ODEProblem（fun ，u0，tspan，p）
tic（）
sol = solve（prob，reltol = 1e-8）
toc（）
plot（sol）
推荐的解决方案
这将使用事件并集成sep
import matplotlib.pyplot as plt
从scipy导入numpy为np
。集成importsolve_ivp
a = 0.02
b = 0.2
c = -65
d = 8
i = 0 0
p = [a， b，c，d，i]
＃定义事件函数并将其设为终端事件
def event（t，u）：
返回u [0]-30
event.terminal =真
＃定义微分方程
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]）]
返回du
u = [0,0]
ts = []
ys = []
t = 0
往往= 100
而True：
sol = resolve_ivp（有趣，（t，趋势），u，事件=事件）
ts.append（sol.t）
ys.append（sol.y）$ b如果sol.status == 1，则为$ b：＃发生事件
＃集成的新开始时间
t = sol.t [-1]
＃重置初始状态
u = sol .y [:, -1] .copy（）
u [0] = p [2] #reset t o -65
u [1] = u [1] + p [3]
其他：
休息
图= plt.figure（）
ax = fig.add_subplot（1，1，1）
＃我们必须将绘制
的单独模拟结果拼接在一起ax.plot（np.concatenate（ts），np.concatenate（ys，axis） = 1）.T）
myleg = plt.legend（['v'，'u']）
最小更改解决方案
似乎您的方法对于 solve_ivp 来说很好用。
警告我认为在Julia和 solve_ivp 中，正确的处理方式事情是使用事件。我相信下面的方法依赖于实现细节，即传递给函数的状态向量与内部状态向量是同一对象，这使我们可以就地对其进行修改。如果是副本，则此方法无效。此外，这种方法不能保证求解器会采取足够小的步长，以使达到极限的正确点被踩到。使用事件将使此方法更正确并且更易于推广到其他在不连续性之前可能具有较低梯度的微分方程。
import matplotlib.pyplot as plt
从numpy导入为np
从matplotlib.ticker导入FormatStrFormatter
从scipy.integrate import resolve_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]
如果u [0]< 30：＃检查是否已达到阈值
du [0] =（0.04 * u [0] + 5）* u [0] + 150-u [1]-p [4]
du [1] = p [0] *（p [1] * u [0] -u [1]）$ b $ b else：
u [0] = p [2]＃重置为-65
u [1] = u [1] + p [3]
return du
y0 = [0,0]
tspan =（0,100）
sol = resolve_ivp（fun，tspan，y0）
图= plt.figure（）
轴= fig.add_subplot（1、1、1）
plt.plot（sol.t，sol.y [0，：]，'k'，线宽= 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））
结果
I am trying to solve a differential equation with Python.
In this two system differential equation if the value of first variable (v) is more than a threshold (30) it should be reset to another value (-65). Below I put my code. The problem is that the value of first variable after reaching 30 remains constant and won't reset to -65. These equations describe the dynamics of a single neuron. The equations are taken from this website and this 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))
The solution looks like:
Here is the correct solution by Julia, here u1 represent v:
This is the 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)
Recommended solution
This uses events and integrates separately after each discontinuity.
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'])
Minimum change "solution"
It appears as though your approach works just fine with solve_ivp.
Warning I think in both Julia and solve_ivp, the correct way to handle this kind of thing is to use events. I believe the approach below relies on an implementation detail, which is that the state vector passed to the function is the same object as the internal state vector, which allows us to modify it in place. If it were a copy, this approach wouldn't work. In addition, there is no guarantee in this approach that the solver is taking small enough steps that the correct point where the limit is reached will be stepped on. Using events will make this more correct and generalisable to other differential equations which perhaps have lower gradients before the discontinuity.
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
这篇关于如何在微分方程（SciPy）中使用if语句？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
