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)
