How to prevent odeint from giving me solutions that blow up

I’m trying to solve and plot ODE using Scipy’s odeint using different initial conditions. This is done in the code below. Note that for three of the initial conditions (2, 4, and 6), the solutions die out, then the graphs start looking weird for these 3 solutions (in the plot, it’s most notable for ic/N0 = 6, which corresponds to the green curve, but you can also see some blue and orange on the fringes at the bottom). How do I fix this so that for those solutions that do die out eventually, I just get these curves, without the weird behavior afterwards? Obviously, one way to do this would be to just stop plotting the curves based on when the solution goes from positive to negative, but I’m wondering if there’s a more elegant way of doing this.

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)

EDIT: Here’s a not very elegant way of solving this problem:

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)

EDIT: I first tried and asked about how to solve this using solve_ivp but things got more complicated

Answer

As you are not interested in the negative range, let’s cut off the ODE function at some negative value for u. In general ODE with a bounded or linearly growing right side exist for all times and are also numerically benign.

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)

This results without any other changes in the plot

Advertisement

Answer