So far, I have been using scipy.integrate.odeint as my "workhorse" ODE solver. My current problem requires that I solve a large system (up to ~5000) ODEs.

Here's the general workflow of the odeint program, that I have noticed:

  1. get initial values and function to calculate Jacobian (or RHS) of the ODE system at $t = t_0$
  2. make estimate of $t_n = t_{n-1} + \delta t$, using standard numerical quadrature methods (requires 1 call to Jacobian function)
  3. tackle non-linearity by "converging" solution at $t_1$ (might require many calls, say $N_n$, to Jacobian function

So, to proceed at each step $\delta t$, we need at least $N_n + 1$ calls to the Jacobian. If calculating the Jacobian requires some effort, then this can be computationally expensive.

One way to reduce this cost that I have been thinking about is by considering some coefficients/parameters as staying constant (almost like a "quasi-equilibrium" assumption over some small time scale?) over $\delta t$. That way, while the ODE system is solved over $\delta t$, the repeated calls to the Jacobian function will not be as expensive, since it won't have to calculate certain things, until the next time step.

My naive implementation of such a simplification used a method similar to the one described here: scipy.integrate.odeint: how can odeint access a parameter set that is evolving independently of it?

Basically, in order to solve the ODE system over a the interval $[0, 10]$ (as an example), I would:

  1. consider first the interval $[0, 1]$, and pre-calculate the coefficients/parameters to be used for this time step
  2. pass the pre-calculated values from step 1. as constant arguments to the Jacobian function for solving the ODE system over $[0, 1]$,
  3. solve the ODE system over $[0, 1]$ but discard all results except for the data obtained for $t = 1$
  4. use the data obtained for $t = 1$ to pre-calculate the coefficients/parameters to be used over the time-step $[1, 2]$
  5. similar to steps 2-3, call odeint again
  6. rinse and repeat

When my supervisor found out what I was doing, similar to the answer for the question I linked above, he warned that this implementation might be both inefficient (since I am breaking the internal flow of odeint) and inaccurate (since I am sort of artificially introducing discontinuities constantly).

How could I go about implementing the "quasi-equilibrium" simplification correctly, if it is possible?


1 Answer 1


I recommend that you check out Rosenbrock or Rosenbrock-W methods. These are multi-stage (like Runge-Kutta) methods that compute the Jacobian only once per step. The methods have more order conditions than diagonally-implicit Runge-Kutta methods (due to the lagged and/or inexact Jacobian), but are otherwise similar. There is a chapter in Hairer and Wanner's second volume that explains the order conditions and discusses some classes of methods.

You can access the PETSc implementation, which implements many Rosenbrock and Rosenbrock-W methods, via petsc4py. See the Oregonator example for usage.

  • $\begingroup$ Thanks so much! Looks like I have got some reading to do :) $\endgroup$
    – bzm3r
    Oct 6, 2014 at 23:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.