<p>正如在另一个答案中指出的,每次调用<code>odeint</code>中的ODE函数时,时间可能不会严格增加,特别是对于僵硬的问题</p>
<p>处理ode函数中这种不连续性的最可靠的方法是使用事件来查找示例中<code>(A * x - B * x * y) - 0.6</code>的零的位置。对于不连续解,使用终端事件精确地在零处停止计算,然后更改ode函数。在<code>solve_ivp</code>中,可以使用<code>events</code>参数执行此操作。参见<a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html" rel="nofollow noreferrer">solve ivp documentation</a>,特别是与炮弹轨迹相关的示例<code>odeint</code>不支持事件,并且<code>solve_ivp</code>有一个可用的LSODA方法,该方法调用与<code>odeint</code>相同的Fortran库</p>
<p>下面是一个简短的示例,但是您可能希望在解决<code>sol2</code>之前额外检查<code>sol1</code>是否到达了终端事件</p>
<pre><code>from scipy.integrate import solve_ivp
tend = 10
def discontinuity_zero(t, y):
return y[0] - 10
discontinuity_zero.terminal = True
def ode_func1(t, y):
return y
def ode_func2 (t, y):
return -y**2
sol1 = solve_ivp(ode_func1, t_span=[0, tend], y0=[1], events=discontinuity_zero, rtol=1e-8)
t1 = sol1.t[-1]
y1 = [sol1.y[0, -1]]
print(f'time={t1} y={y1} discontinuity_zero={discontinuity_zero(t1, y1)}')
sol2 = solve_ivp(ode_func2, t_span=[t1, tend], y0=y1, rtol=1e-8)
plt.plot(sol1.t, sol1.y[0,:])
plt.plot(sol2.t, sol2.y[0,:])
plt.show()
</code></pre>
<p>这将打印以下内容,其中不连续时间精确到7位</p>
<pre><code>time=2.302584885712467 y=[10.000000000000002] discontinuity_zero=1.7763568394002505e-15
</code></pre>