如何防止 odeint 给我解决方案

Question

我正在尝试使用 Scipy 的 odeint 使用不同的初始条件来求解和绘制 ODE。这是在下面的代码中完成的。请注意，对于三个初始条件（2、4 和 6），解决方案消失了，然后对于这 3 个解决方案，图表开始看起来很奇怪（在图中，最值得注意的是 ic/N0 = 6，它对应于绿色曲线，但您也可以在底部的边缘看到一些蓝色和橙色）。我该如何解决这个问题，以便对于那些最终会消失的解决方案，我只是得到这些曲线，之后没有奇怪的行为？显然，这样做的一种方法是根据解决方案从正变为负的时间停止绘制曲线，但我想知道是否有更优雅的方法来做到这一点。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
ic = [2,4,6,8,10,12,14,16,18,20]
def de(t, u):
    return u*(1-u/12)-4*np.heaviside(-(t-5), 1)

plt.xlim([0, 10])
plt.ylim([0, 20])
for N0 in ic:
    N1 = odeint(de, N0, np.linspace(0, 10, 10000), tfirst=True)
    plt.plot(np.linspace(0,10,10000), N1)

编辑： 这是解决此问题的一种不太优雅的方法：

ic = [8,10,12,14,16,18,20]
ic_to_zero = [2,4,6]

plt.xlim([0, 10])
plt.ylim([0, 20])
for N0 in ic:
    N1 = odeint(de, N0, np.linspace(0, 10, 10000), tfirst=True)
    plt.plot(np.linspace(0,10,10000), N1)
for N0 in ic_to_zero[0:1]:
    N1 = odeint(de, N0, np.linspace(0, 2, 10000), tfirst=True)
    plt.plot(np.linspace(0,2,10000), N1)
for N0 in ic_to_zero[1:2]:
    N1 = odeint(de, N0, np.linspace(0, 2, 10000), tfirst=True)
    plt.plot(np.linspace(0,2,10000), N1)
for N0 in ic_to_zero[2:3]:
    N1 = odeint(de, N0, np.linspace(0, 4, 10000), tfirst=True)
    plt.plot(np.linspace(0,4,10000), N1)

编辑：我第一次尝试和asked 讨论如何使用solve_ivp 解决这个问题，但事情变得更加复杂

Answer 1

由于您对负值范围不感兴趣，所以让我们在 u 的某个负值处切断 ODE 函数。一般来说，右侧有界或线性增长的 ODE 始终存在，并且在数值上也是良性的。

ic = [2,4,6,8,10,12,14,16,18,20]
def de(t, u):
    u = max(-10,u) # replace f(t,u) with f(t,-10) for u<-10
    return u*(1-u/12)-4*np.heaviside(-(t-5), 1)

plt.xlim([0, 10])
plt.ylim([0, 20])
t_plot = np.linspace(0, 10, 10000)
for N0 in ic:
    N1 = odeint(de, N0, t_plot, tfirst=True)
    plt.plot(t_plot, N1)

此结果在情节中没有任何其他变化

Answer 2

我相信在这种情况下Lutz Lehmann's solution 更好。但是由于我已经对其进行了编码，因此我发布了我的解决方案，该解决方案在积分曲线超出范围后也会停止。也许这对有一天来这个帖子的人有用。

ic = [2,4,6,8,10,12,14,16,18,20]

def de(t, u):
    return u*(1-u/12)-4*np.heaviside(-(t-5), 1)

for N0 in ic:
    # supressing ODEintWarning
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        N1 = odeint(de, N0, np.linspace(0, 10, 10000), tfirst=True).reshape(-1)
    cond = (0 > N1) | (N1 > 20)
    stop_at = len(N1)-1 if np.all(~cond) else np.argmax(cond)
    plt.plot(np.linspace(0,10,10000)[:stop_at], N1[:stop_at])