I have a system of 4 coupled ODEs that model an astrophysical system. The dynamic range of the state variables (masses) and the independent time variable are huge – the masses are expected to go from ~0 grams in the early universe to ~10^43 grams ~10 billion years later. So far I have been manually using a forward Euler integration loop to “brute-force” solve for the evolution of the state variables – my results match theoretical calculations, except at very early times where there can be oscillations and the masses can become negative (though I force them to be >= 0).
I have not thought deeply about the stability/error of my ODEs at very early times. However, recently I have become interested in porting my code to standard ODE solvers just to see what happens. I started by trying Python’s scipy.integrate.solve_ivp but all of the methods fail because “'Required step size is less than spacing between numbers.” – they are basically hanging at early times.
So today I decided to take a shot at diffeqpy and Julia (first time). Unfortunately all of the diffeqpy auto and implicit solvers I tried similarly failed at early times with “Warning: Instability detected. Aborting”. So now I’m wondering – is it a coincidence that my fixed-step forward Euler integration succeeds in getting over the initial hump and matching theoretical expectations at late times, or is my problem just not suited for higher order predictor-corrector methods due to the large dynamic range of variables? Are there ODE systems where stiffness early on basically requires using lower order methods and if so how do I verify my “manual” Euler integration? Note that I need to be able to solve my ODE system in ~1 sec so even if I were to switch to a general solver, it would need to not get stuck early on (I’m not interested in the exact values of the state variables at very early times anyway).
I tried to create a minimal working example below (purely in Julia) that gives a flavor of what I’m working with. In this simple case I’m able to show that the Julia solver and my brute-force Euler method agree decently. However my actual ODEs (in Python) are more complicated and are fed interpolated values of external parameters as a function of time (these can be kind of noisy which I suspect may also be contributing to the instability at early times). I am still trying to think of ways to “break” this minimal example so that it mimics my real ODE system, but in the meantime I’d appreciate any suggestions for how to overcome early stiffness in ODEs with general solvers and if I should just stick with my brute-force forward Euler approach.
# load needed packages
using Random
using Distributions
using Interpolations
using DifferentialEquations
using Plots
# define units
u_s = 3.15e7 # seconds in 1 year
u_Msun = 1.989e33 # grams in 1 solar mass
# initialize a list of 500 times at which we are given input parameters
times = LinRange(0.1e9*u_s,10e9*u_s,500) # 0.1-10 Gyr since big bang in seconds
# assume at every time we are given the amount of new mass flowing into the system
# for simplicity generate 500 random log-normal values
d = LogNormal() # default mu=0, sigma=1
mdot_grow_array = rand(d,500) * u_Msun / u_s # Msun/yr to grams/sec
# do a linear interpolation so the ODE solver has a smooth function to get Mdot_grow at any time t
interp_grow = LinearInterpolation(times,mdot_grow_array,extrapolation_bc=Flat())
# assume that some fraction of mass also gets converted into another form at every time t
# for simplicity make these fractions random uniform numbers between 0 and 0.1
d2 = Uniform(0.0,0.1)
frac_convert_array = rand(d2,500) / (1e9*u_s) # conversion fraction of M1->M2 per second (over 1 Gyr)
# and again interpolate the conversion fractions so the ODE solver has a smooth function
interp_convert = LinearInterpolation(times,frac_convert_array,extrapolation_bc=Flat())
# set up the integrator function following the ODE Example in the DifferentialEquations documentation
function integrator!(du,u,p,t)
# the differential equations are
# dM1/dt = mdot_grow - frac_convert*M1 (mass of M1 grows by new accretion but decreases due to conversion to M2)
# dM2/dt = frac_convert*M1 (mass of M2 steadily grows due to conversion from M1 to M2)
# get the interpolated values of mdot_grow and frac_convert at current time t
mdot_grow_now = interp_grow(t)
frac_convert_now = interp_convert(t)
# compute mdot_convert as frac_convert * M1
mdot_convert_now = frac_convert_now * u[1]
# return/update the derivatives for u[1] (aka dM1/dt) and u[2] (aka dM2/dt)
du[1] = mdot_grow_now - mdot_convert_now
du[2] = mdot_convert_now
end
# finally set up initial conditions, timespan and run the ODE solver
u0 = [0.0;0.0] # initially we have zero mass for both M1 and M2
tspan = (times[1],times[end])
prob = ODEProblem(integrator!,u0,tspan)
sol = solve(prob)
#### Brute-force forward Euler integration
# initialize arrays of zeros for values of M1 and M2 that we will update throughout integration
M1_vals = zeros(size(times));
M2_vals = zeros(size(times));
# begin the Euler integration loop
for tnum in range(1,size(times)[1]-1)
print("Working on time ",tnum,"\n")
# get current and next time, and the timestep in between in seconds
tnow, tnext = times[tnum], times[tnum+1]
dt = tnext-tnow # in seconds
# get current values of state variables
M1_now = M1_vals[tnum]
M2_now = M2_vals[tnum]
# get the interpolated value of new mass flowing in at this time and add it to M2_now
mdot_grow_now = interp_grow(tnow)
M1_now = M1_now + mdot_grow_now * dt
# compute conversion rate from M1 to M2 (depends on the updated M1 value)
frac_convert_now = interp_convert(tnow)
mdot_convert_now = frac_convert_now * M1_now
# subtract the converted mass from M1 and add it to M2
M1_now = M1_now - mdot_convert_now * dt
M2_now = M2_now + mdot_convert_now * dt
# append the final values of the state variables at this timestep to our lists
M1_vals[tnum+1] = M1_now
M2_vals[tnum+1] = M2_now
end
plot(sol)
plot!(times,M1_vals,linecolor=:green,linestyle=:dash,label="Euler M1(t)")
plot!(times,M2_vals,linecolor=:black,linestyle=:dash,label="Euler M2(t)")
