# How to deal with solving coupled ODE systems where variables are updated multiple times within each timestep?

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
break
else:
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)

• 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.
– Community Bot
Mar 10 at 22:06
• 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. Mar 11 at 8:47
• 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. Mar 11 at 13:08
• 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? Mar 11 at 17:35
• 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 Mar 11 at 17:44

• 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. Mar 12 at 4:29