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