Scipy odeint gives lsoda warning

Question

I am completely 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 is a 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 a certain set of parameters when playing with mxstep in odeint (I also tried hmin and hmax , but did not notice any difference), although the error persists, my graphs look good and do not affect, but in most cases they are. Sometimes the error that I cause runs with the option odeint full_output=1 , and in doing so I get:

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

I don’t understand what this means when searching.

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 this warning, lsoda t refers to the current time value, and h refers to the current size of the integration step. The step size has become so close to zero that the current time plus the step size is estimated to be equal to the current time due to rounding errors (i.e. r1 + r2 == r1 ). This problem usually occurs when the ODEs that you integrate behave badly.

On my machine, only a warning appears only when calculating soln2 . Here I built the values ​​of each of the parameters in the area where the warnings occur (note that I switched to linear axes for clarity). I noted a time step when the lsoda error first appeared (in r1 = 0.2135341098625E-06 ) with a red line:

It is no coincidence that the appearance of a warning message coincides with a “kink” in E!

I suspect that it happens that the bend is a singularity in the gradient of E, which forces the solver to take smaller and smaller steps until the step size reaches the limits of precision with a floating point. In fact, you can see another inflection point in D that coincides with the “swing” in B, probably caused by the same phenomenon.

For some general tips on how to handle features when integrating ODEs, see section 5.1.2 here (or any decent ODE tutorial).

---

You got an error with full_output=True , because in this case, odeint returns a tuple containing the solution, and a dict containing additional output. Try unpacking the tuple as follows:

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