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)