scipy odeint: excess work done on this call depending on initial values even with analytically solvable ODE

Question

I am trying to solve a differential equation in the form: dx/dt = funct(x) using scipy odeint.

However, for some initial values, I get a "ODEintWarning: Excess work done on this call", even if the ODE itself seems to be solvable analytically. Why does the numerical ODE solver give me a warning, and how can I make it work instead of jumping to x > 3e^7?

Here is a minimal (non) working example:

x0 =0.49 # works, but changing to 0.51 gives the ODEint warning
t = np.linspace(0,10,100)
def funct(x,t):
    return 2*x**2-x
x, out = integrate.odeint(funct, x0, t, full_output = 1)
plt.plot(t,x)

I am not too good at math, but I think the analytical answer should be: x = 1/(2 + C * exp(t)) with C being any constant.

If so, I think using initial value of x0=0.51 should also be numerically solvable because C exists (C = -0.04).

The closest question I could find was: scipy odeint: excess work done on this call and very sensitive to initial value However, the ODE in the previous question blows up analytically, whereas mine does not.

Edit: My real function and x are vectors and matrices about networks, not an analytically solvable one.

Answer 1

If I am not mistaken, for $x>0.5$, your example ODE yields $dx/dt>0$ and thus $x$ will increase and its time derivative as well. Therefore the solution diverges, even faster than a simple exponential after some time (when $x^2 \gg x$). If you start at $x=0.5$, the solution is a constant. For $x<0.5$, the solution converges to 0.

In the case of a divergence, the error estimate that is used to dynamically adapt the time step will require time steps smaller and smaller, and at one point lower than $10^{-15}$, or at least such that $t+dt=t$ at machine precision. The solver then assumes that there is no point in trying to pursue the integration and notifies you. That behaviour is absolutely normal.

Answer 2

To be more precise in the influence of the initial value, the ODE $\dot x=2x^2-x$ is Bernoulli. With $u=x^{-1}$ one gets $$ \dot u=-x^{-2}\dot x=u-2\implies u=2+(u_0-2)e^t\implies x=\frac{x_0}{2x_0+(1-2x_0)e^t}. $$ The denominator has a pole for $x_0>\frac12$ at $t_{pole}=\ln(\frac{x_0}{x_0-\frac12})$. For $x_0=0.51$ this gives $t_{pole}=\ln(51)=3.931825...$ This singularity provides an insurmountable obstacle for the integrator. There could be more informative error messages, but odeint just reports the failure of the step-size controller.