Solve_ivp From Scipy Does Not Integrate The Whole Range Of Tspan

I'm trying to use solve_ivp from scipy in Python to solve an IVP. I specified the tspan argument of solve_ivp to be (0,10), as shown below. However, for some reason, the solutions

Solution 1:

What is going wrong

You should always look at what the solver returns. In this case it gives

message: 'Required step size is less than spacing between numbers.'

Think of the process of solving your initial value problem with scipy.integrate.solve_ivp as repeatedly estimating a direction and then going a small step in that direction. The above error means that the solutions to your equation change so fast that taking the minimal step size that machine accuracy allows is too far. But your equation is simple enough that at least for t =< 5 where 4*np.heaviside(-(t-5), 1) always gives 4 it can be solved exactly/symbolically. I will explain more for t > 5 later.

Symbolic Solution

Sympy can solve your differential equation. While you can provide it an initial value it would have taken much longer to solve it once for each of your initial values. So instead I told it to give me all solutions and then I calculated the parameters C1 for your initial value separately.

import numpy as np
import matplotlib.pyplot as plt
from sympy import *

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

f = symbols("f", cls=Function)
t = symbols("t")
eq = Eq(f(t).diff(t),f(t)*(1-f(t)/12)-4)
base_sol = dsolve(eq)
c1s = [solve(base_sol.args[1].subs({t:0})-ic) for ic in ics]
# Apparently sympy is unhappy that numpy does not supply a cotangent. # So I do that manually.
sols = [lambdify(t, base_sol.args[1].subs({symbols('C1'):C1[0]}), 
modules=['numpy', {'cot':lambda x:1/np.tan(x)}]) for C1 in c1s]
t = np.linspace(0, 5, 10000)

for sol in sols:
    y = sol(t)
    mask = (y > -5) & (y < 20)
    plt.plot(t[mask], y[mask])

At first glance the picture looks odd. Especially the blue and orange straight line part. This is just due to the values lying outside the masked range so matplotlib connects them directly. What is actually happening is a sudden jump. That jumped tipped off the numeric ode solver earlier. You can see it even more clearly when you make sympy print the first solution.

The tangent is known to have a jump at pi/4 and if you solve the argument of the tangent above you get 2.47241377386575. Which is probably where your plotting stopped.

Now what about t>5?

Unfortunately your equation is not continuous in t=5. One approach would be to solve the equation for t>5 separately for the initial values given by following the solutions of the first equation. But that is an other question for an other day.