Number of equations and precision of SciPy's integrate.odeint()

Question

Is there any reason why SciPy's integrate.odeint() should become less precise when the number of equations increase? I'm trying to solve these two sets of differential equations:

$\frac{dy_1}{dx} = x^2$ — (1)

$\frac{dy_2}{dx} = x^3$

and,

$\frac{dy_1}{dx} = x^2$ — (2)

Since the two equations in (1) are uncoupled, $y_1$ from integrating (1) and (2) should be the same. I'm implementing the two in the following snippet:

import numpy as np
from scipy.integrate import odeint

def f(y, x):
    return [x**2, x**3]

def g(y, x):
    return x**2

a = odeint(f, [0.0, 0.0], np.arange(0, 5, 0.0001))
b = odeint(g, 0.0, np.arange(0, 5, 0.0001))

print a[-1][0], b[-1][0], abs(a[-1][0] - b[-1][0])

Running the above code gives me:

41.6641667161 41.6641667354 1.93298319573e-08

The difference seems to be very insignificant here. But in cases where the number of equations becomes large (for e.g. in Lyapunov exponent calculations), it appears to cause significant differences.

What could be going on here?

---

UPDATE: Like @horchler explained, using numpy.linspace instead of numpy.arange did increase the accuracy of the final values, but the difference between two answers is of the same order:

41.6666666618 41.6666666811 1.93298248519e-08

Answer 1

To begin, it might be a better idea to use linspace in this case if you care about the final value or just let the integrator choose the intermediate steps. That will give you accurate final values, but the difference in precision between the two cases will remain.

There are a few interrelated reasons why this is occurring. Just because two equations are uncoupled does not mean that odeint solves them separately. The solver has to choose one time step on each iteration and this will depend on the behavior of all of the equations. In the case of your simpler scalar quadratic equation, the solver may be able to take much larger time steps, whereas the cubic equation in the other case grows much more quickly and may cause the solver to reduce step sizes or even switch between modes. The time steps chosen depend on (amongst other things) the absolute and relative tolerances which you can adjust.

Answer 2

scipy.integrate.odeint is really a wrapper to LSODA, part of ODEPACK, written in the 1980s. It uses some heuristics to switch between Adams (nonstiff) and BDF (stiff) methods, depending on diagnostics used to determine the stiffness of the ODE system. Since BDF methods are implicit, presumably the Adams methods used are Adams-Moulton methods (which are also implicit, as opposed to Adams-Bashforth methods, which are explicit).

Without a more concrete example comparing two large systems, the errors you see could be caused by:

• stiffness/ill-conditioned Jacobian matrix in your right-hand side
• roundoff from doing substantially more floating point operations (say, in the limit as your larger system becomes infinitely larger than the smaller system)
• error tolerances, depending on how you set them for any "new" variables