如何避免Scipy的odeint求解常微分方程时得到发散的解

使用Scipy的odeint求解不同初始条件下的常微分方程并绘制结果曲线时，初始条件为2、4、6的三组解在衰减到0之后，曲线会出现异常表现。希望修复该问题，让最终会衰减到0的解只保留正常的曲线段，不需要手动为每个初始条件单独设置求解时间区间，且适配odeint接口，不使用solve_ivp。

原始代码如下：

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)
plt.show()

现有手动分段处理的方案不够灵活，无法适配任意初始条件的情况。

核心思路是在求解完成后自动截断解的有效区间：

• 统一生成全局时间数组，避免循环中重复计算
• 对每个初始条件的求解结果，找到解首次小于等于0的位置
• 仅绘制从起点到该位置的有效数据段，无需提前拆分初始条件组，也不需要手动设置每个组的求解区间

如果担心浮点数计算精度问题，可以将判断条件调整为N1 <= 1e-6，过滤接近0的极小数值的干扰。

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)

# 统一生成时间数组，避免循环中重复创建
t = np.linspace(0, 10, 10000)

plt.xlim([0, 10])
plt.ylim([0, 20])
for N0 in ic:
    # ravel把odeint返回的二维数组转为一维，方便处理
    N1 = odeint(de, N0, t, tfirst=True).ravel()
    # 找到第一个解<=0的索引，所有解都为正的话返回0
    first_neg_idx = np.argmax(N1 <= 0)
    if first_neg_idx > 0:
        # 存在负数值，截断到第一个负数值之前
        plt.plot(t[:first_neg_idx], N1[:first_neg_idx])
    else:
        # 所有解都为正，直接绘制全部
        plt.plot(t, N1)
plt.show()

该方案可以自动适配任意初始条件，不需要手动调整每个初始值对应的求解区间，逻辑简洁高效，完全基于odeint接口实现。

内容的提问来源于stack exchange，提问作者An Ignorant Wanderer