Scipy odeint giving lsoda warning

Question

I am totally new to coding and I want to solve these 5 differential equations numerically. I took a python template and applied it to my case. Here's the simplified version of what I wrote:

import numpy as np
from math import *
from matplotlib import rc, font_manager
import matplotlib.pyplot as plt
from scipy.integrate import odeint

#Constants and parameters
alpha=1/137.
k=1.e-9     
T=40.    
V= 6.e-6
r = 6.9673e12
u = 1.51856e7

#defining dy/dt's
def f(y, t):
       A = y[0]
       B = y[1]
       C = y[2]
       D = y[3]
       E = y[4]
       # the model equations
       f0 = 1.519e21*(-2*k/T*(k - (alpha/pi)*(B+V))*A) 
       f1 = (3*B**2 + 3*C**2 + 6*B*C + 2*pi**2*B*T + pi**2*T**2)**-1*(-f0*alpha/(3*pi**3) - 2*r*(B**3 + 3*B*C**2 + pi**2*T**2*B) - u*(D**3 - E**3))
       f2 = u*(D**3 - E**3)/(3*C**2)
       f3 = -u*(D**3 - E**3)/(3*D**2)
       f4 = u*(D**3 - E**3)/(3*E**2) + r*(B**3 + 3*B*C**2 + pi**2*T**2*B)/(3*E**2)
       return [f0, f1, f2, f3, f4]

# initial conditions
A0 = 2.e13
B0 = 0. 
C0 = 50.           
D0 = 50.
E0 = C0/2.

y0 = [A0, B0, C0, D0, E0]       # initial condition vector
t  = np.linspace(1e-15, 1e-10, 1000000)   # time grid

# solve the DEs
soln = odeint(f, y0, t, mxstep = 5000)
A = soln[:, 0]
B = soln[:, 1]
C = soln[:, 2]
D = soln[:, 3]
E = soln[:, 4]

y2 = [A[-1], B[-1], C[-1], D[-1], E[-1]]  
t2  = np.linspace(1.e-10, 1.e-5, 1000000)  
soln2 = odeint(f, y2, t2, mxstep = 5000)
A2 = soln2[:, 0]
B2 = soln2[:, 1]
C2 = soln2[:, 2]
D2 = soln2[:, 3]
E2 = soln2[:, 4]

y3 = [A2[-1], B2[-1], C2[-1], D2[-1], E2[-1]]     
t3  = np.linspace(1.e-5, 1e1, 1000000)
soln3 = odeint(f, y3, t3)
A3 = soln3[:, 0]
B3 = soln3[:, 1]
C3 = soln3[:, 2]
D3 = soln3[:, 3]
E3 = soln3[:, 4]

#Plot
rc('text', usetex=True)
plt.subplot(2, 1, 1)
plt.semilogx(t, B, 'k')
plt.semilogx(t2, B2, 'k')
plt.semilogx(t3, B3, 'k')

plt.subplot(2, 1, 2)
plt.loglog(t, A, 'k')
plt.loglog(t2, A2, 'k')
plt.loglog(t3, A3, 'k')

plt.show()

I get the following error:

lsoda--  warning..internal t (=r1) and h (=r2) are         
   such that in the machine, t + h = t on the next step 
   (h = step size). solver will continue anyway         
  In above,  R1 =  0.2135341098625E-06   R2 =  0.1236845248713E-22

For some set of parameters, upon playing with mxstep in odeint (also tried hmin and hmax but didn't notice any difference), although the error persists, my graphs look good and are not impacted, but most of the times they are. Sometimes the error I get asks me to run with the odeint option full_output=1 and in doing so I obtain:

A2 = soln2[:, 0]
TypeError: tuple indices must be integers, not tuple

I don't get what this means when searching for it.

I would like to understand where the problem lies and how to solve it. Is odeint even suitable for what I'm trying to do?

Answer 1

In that lsoda warning, t refers to the current time value and h refers to the current integration step size. The step size has become so close to zero that the current time plus the step size evaluates as equal to the current time due to rounding error (i.e. r1 + r2 == r1). This sort of problem usually occurs when the ODEs you are integrating are badly behaved.

On my machine the warning message only seems to occur when computing soln2. Here I've plotted the values of each of the parameters in the region where the warnings are occurring (note that I've switched to linear axes for clarity). I've marked the time step where the lsoda error first appeared (at r1 = 0.2135341098625E-06) with a red line:

It's no coincidence that the appearance of the warning message coincides with the 'kink' in E!

I suspect that what's happening is that the kink represents a singularity in the gradient of E, which is forcing the solver to take smaller and smaller steps until the step size reaches the limits of floating point precision. In fact, you can see another inflection point in D which coincides with a 'wobble' in B, probably caused by the same phenomenon.

For some general advice on how to deal with singularities when integrating ODEs, take a look at section 5.1.2 here (or any decent textbook on ODEs).

---

You were getting an error with full_output=True because in this case odeint returns a tuple containing the solution and a dict containing additional output information. Try unpacking the tuple like this:

soln, info = odeint(..., full_output=True)