How does odeint() from scypy python module work?

Question

I am a physics student interested in solving ODEs numerically. I usually write my own solvers in C using Runge–Kutta methods.

I recently learned Python, and I used SciPy’s odeint function to solve ODEs. But I am worried about how the function algorithm works, because it not take a step size argument. So, how can I learn how it works? How can I know what is the precision of their results?

I consulted this documentation, but it does not offer very much information, and I don’t really understand the optional arguments that they describe.

Answer 1

SciPy has three modules for integrating ODEs:

• scipy.integrate.odeint, which uses the LSODA algorithm.
• scipy.integrate.ode, which supports four different backends (LSODA, DoPri5, DoP853, and VODE).
• scipy.integrate.solve_ivp, which supports five backends (LSODA, DoPri5, Bogacki–Shampine, Radau, and a backwards-differentiation method).

Note that while asking about the first module, you used the documentation for the second.

All of these solvers adapt the step size. This means that they estimate the integration error, e.g., by comparing the result of two simple integrators (which are intelligently intertwined to save runtime). This estimate is then used to adapt the step size such that the estimated error is below a user-defined threshold specified by the arguments rtol and atol. These arguments allow you to control the precision.