如何防止 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

我相信在这种情况下更好。但是因为我已经对它进行了编码，所以我 post 我的解决方案也会在积分曲线超出范围后停止它。也许这对有一天来到这个话题的人有用。

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])

如何防止 odeint 给我解决方案失败

How to prevent odeint from giving me solutions that blow up

python

scipy

differential-equations