I'm solving a system of coupled ODEs using Euler integration for simplicity. To make this concrete, please see the (extremely simplified) minimal working example below in Python. Imagine we have a box that initially has 0 mass and 0 energy. At each timestep, we add some mass to the box from the outside. Some fraction of the box mass at each timestep is converted into energy -- this means we need to decrease the mass and increase the energy variables. But then later in the same timestep we also allow some fraction of the energy to be converted back into mass -- so we update the mass and energy variables again. At the end we record only the final values of mass and energy at the end of each timestep but clearly there was "evolution" within the timestep (in the example below, the mass variable changes 3 times and energy changes 2 times in each timestep).

What is the proper way to think about this "time-substep" evolution? When is it important to keep track of the intermediate values of variables? It's not even clear to me what "sub-times" you would assign to each of the intermediate variable values in between the "big" timestep.

This has to be an extremely common numerical issue -- are there well-defined strategies for dealing with it (and to know when one should or shouldn't care)? I know that something like this is an issue when solving Maxwell's E&M equations and that people typically use "leapfrog integration" -- and maybe there is also a connection to higher-order integration schemes like RK4. But I'm having a hard time understanding the bigger context and would love clarification or references to pedagogical literature or code examples.

import numpy as np

# initialize 100 times
# assume dimensionless code units throughout
times = np.linspace(0,100,100)

# at every timestep, assume some mass gets added into a box from outside
mass_growth_rates = np.abs(np.random.randn(100)) 

# also at each timestep, some fraction of the available mass will get converted into another form (call it energy) 
conversion_fractions = np.random.uniform(0,1,100) 

# finally at each timestep some fraction of that energy will be returned back to the mass reservoir
return_fractions = np.random.uniform(0,1,100) 

# initialize two lists that will record the box mass and energy at each time
mass_array, energy_array = np.array([0.0]), np.array([0.0]) # initial conditions: both are 0

# simple euler integration loop to keep track of mass and energy vs time
for tnum in np.arange(len(times)):
    # compute timestep between current and next time
    if tnum == len(times)-1: # we are on the final time, no more steps to take 
        dt = times[tnum+1] - times[tnum]
    # get the mass and energy at the end of the previous timestep
    mass_now = mass_array[tnum]
    energy_now = energy_array[tnum]

    # add the new mass 
    mass_now = mass_now + mass_growth_rates[tnum]*dt
    # subtract a fraction of mass and add it to energy
    mass_now = mass_now - conversion_fractions[tnum]*mass_now*dt
    energy_now = energy_now + conversion_fractions[tnum]*mass_now*dt
    # subtract some fraction of energy and return it to the mass component
    energy_now = energy_now - return_fractions[tnum]*energy_now*dt
    mass_now = mass_now + return_fractions[tnum]*energy_now*dt
    # append the final values of mass and energy at this time to our list
    mass_array = np.append(mass_array, mass_now)
    energy_array = np.append(energy_array, energy_now)
  • $\begingroup$ Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. $\endgroup$
    – Community Bot
    Commented Mar 10, 2022 at 22:06
  • 2
    $\begingroup$ Why don’t you make your mass and energy are dynamical variables themselves and have their conversions governed by the differential equation, either with given fluctuating rates or as an SDE? Right now, you have a discrete-time model for your mass–energy conversion and try to shoehorn this into a continuous-time solving algorithm (presumably for continuous-time components of your model that we do not see). I suggest that you first think about whether your really need this hybrid approach. If yes, please provide details on the nature and relative timescales of the two components. $\endgroup$
    – Wrzlprmft
    Commented Mar 11, 2022 at 8:47
  • $\begingroup$ So you are solving the ODE system $\dot m(t) = g(t)-c(t)m(t)+r(t)e(t)$, $\dot e(t)=c(t)m(t)-t(r)e(t)$. This you can easily code as right-hand-side function and plug into any ODE solver. $\endgroup$ Commented Mar 11, 2022 at 13:08
  • $\begingroup$ Thanks both. @LutzLehmann it's super helpful to put it into that succinct form. But now what if c(t) and r(t) were not just given like in the above example but instead depended on the value of mass_now or energy_now. Then it becomes ambiguous what value of m(t) and e(t) to use in your simplified notation since c(t) becomes c(t,m) and r(t,m) becomes r(t,m). That was the crux of my question -- doesn't this necessitate "operator splitting" as suggested in the answer below, and would solver packages (say scipy solve_ivp) do this under the hood? $\endgroup$ Commented Mar 11, 2022 at 17:35
  • $\begingroup$ To clarify, it becomes ambiguous because m(t) and e(t) need to be updated multiple times within a single timestep and hence take on multiple values within a single timestep. In some sense this is also imposing some kind of "order of operations" for carrying out the various terms on the RHS of your mdot(t) and edot(t) differential equations, and I don't know if there is a formal way to refer to this requirement for solving ODEs in the literature $\endgroup$ Commented Mar 11, 2022 at 17:44

1 Answer 1


It would have been helpful to see the ODE you're actually trying to solve, but the way I interpret what you write is that the right hand side of your ODE system consists of a number of terms corresponding to different effects, and that you update variables one after the other for these terms. This is generally called "operator splitting" and there is a large amount of literature on the topic.

I will note that operator splitting is not the best way to go about this. The better approach is to start by writing down the complete ODE and then handing the right hand side (and possibly its derivatives) to one of the many well-developed ODE integrator packages.

  • 1
    $\begingroup$ Thanks so much. This is EXACTLY what I was hoping for -- now I know that what I'm asking about is formally called "operator splitting". I agree that I can just use an ODE integrator package like scipy's solve_ivp but wouldn't there be certain coupled ODE systems for which which operator splitting is effectively required? For example, what if the conversion rate in my simple code above was not just given as a function of time but also depended on the value of "mass_now"? $\endgroup$ Commented Mar 11, 2022 at 17:18
  • $\begingroup$ That's just another way of saying that your ODE's right hand side is not just an $f(t)$, but an $f(x(t),t)$. The latter is typically the case for ODEs. $\endgroup$ Commented Mar 12, 2022 at 4:29

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.