Possible to reduce effort needed to solve non-linear ODEs by taking some coefficients/parameters as constant over small time intervals?

Question

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?

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.