ODE system solvable with forward Euler but not general solvers (early stiffness)

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

# 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


plot!(times,M1_vals,linecolor=:green,linestyle=:dash,label="Euler M1(t)")
plot!(times,M2_vals,linecolor=:black,linestyle=:dash,label="Euler M2(t)")

You are correct about your problem. If your parameters are noisy, you aren’t solving an ODE, you’re solving an SDE (stochastic differential equation). The reason more sophisticated solvers are failing is because they are looking at points that are very close together, and seeing almost infinite derivatives due to the noise.

Thanks Oscar_Smith! But in my example code above, I only draw one random number for each of my 500 sample times once at the very beginning just to initialize the growth rate and conversion rate time series. And then I fit a 1D linear interpolation to those so that my integrator function has a “smooth curve” for getting the externally-provided growth and conversion rates at any given time. It’s not like I’m adding (or providing) a completely new random value for the growth rate and conversion rate inside my integrator function at each time – it should be deterministic, no?

Moreover in my actual ODE system my various time series of input parameters are not just uncorrelated random numbers like in the minimal working example above – the input time series (growth rates, conversion rates, other rates) are generated by a known physical process and look to be correlated over some timescale and correlated amongst each other.

If I really am dealing with a system of stochastic differential equations, do you have any recommendations for how I can translate my existing ODEs to SDEs and how to integrate those? How different should I expect the values of the state variables to be vs. time with an SDE solver vs my existing forward Euler integration? The one re-assuring thing is that despite the simplicity of my manual forward Euler integration I do roughly recover the theoretical expectations for the state variables vs time with it (except at very early stiff times where the ODE solvers are also failing).

As an aside I have also been wondering if the huge order of magnitude is throwing off the general ODE solvers somehow, and if say log-transforming the state variables and derivatives would help. That is: solve for log(M1) and log(M2) and have the integrator function return dlog(M1)/dlog(t) and dlog(M2)/dlog(t) (or just dt in the denominators). That way the solvers aren’t working with numbers spanning 40 OOM but just the exponents – and presumably there would be no more risk of getting negative masses (which is a problem even for the Euler integration which I catch with if/else statements). Thoughts on if this would help prevent tiny timesteps and help early convergence?

The Euler integration only works because it’s not actually checking for correctness (i.e. error estimate) and it’s just letting anything go. The adaptive algorithms are failing because your f is only continuous. It’s not differentiable because of the linear interpolation. The highest order you can integrate a function is dependent on the regularity of the ODE being integrated. Here, the regularity is <1 (based on the Holder continuity), and so anything other than Euler is just extra work.

Normally what someone would do in this case is use something like a regression spline to smooth the data. DataInterpolations.jl has a bunch of tools for this kind of stuff.

Lowess is another approach people use, etc. Just some kind of smoothing.

No. In its current form, it’s easiest handled as a random ordinary differential equation (RODE). This is not necessarily easy to write as an SDE because you’d have to derive the construction from a Brownian process.


Wow @ChrisRackauckas thank you! So it turns out that my input parameters as a function of time did have sharp jumps in them (in addition to noise), which meant that the 2nd and higher order derivatives exploded with huge numbers. That was probably causing the ODE solvers to try to take very small timesteps and forcing them to abort. The fact that I was using linear interpolation was also contributing to the problem since as you said, linear interp is not differentiable at the input nodes.

Since my code is in python and I haven’t converted all of it to Julia, what I ended up doing in python is:

  1. take the log10 of any input time series that is positive everywhere and discontinuous (or large OOM)
  2. convolve each of the log10 input parameter time series by a gaussian (using scipy.ndimage.gaussian_filter1d with sigma~10)
  3. switch from scipy.interpolate.interp1d(kind=‘linear’) to scipy.interpolate.UnivariateSpline(k=3,s=0) so it’s a cubic interpolating (not smoothing) function which should be twice differentiable.

Those 3 changes allow scipy.integrate.solve_ivp() to successfully run with the explicit methods RK23, RK45 and DOP853 but none of the other implicit methods work. The explicit solutions agree with my manual Euler solution but are better behaved at early times (no oscillations or negative masses). However the explicit methods take a very long time (about 4-5 min) because they are still taking extremely small timesteps early on. I use the default absolute tolerance = 1e-6 and relative tolerance = 1e-3.

What’s even more troubling is that the diffeqpy Julia solvers fail with “Warning: Instability detected. Aborting” (I made sure to change f(t,y) to f(y,p,t) for diffeqpy instead of solve_ivp). I tried AutoTsit5(Rosenbrock23()), Tsit5() and even RK4(). I didn’t play with the ODE solver options.

Any idea what else could be going on and why this would work in scipy but not diffeqpy/Julia? I have a feeling it has something to do with my smoothing+interpolation algorithm still not being good enough – I still see wiggles (discontinuities) in the 2nd and 3rd (numerical) derivatives of my cubic interpolating splines for some reason – see attached figure. That is probably throwing off how the solvers compute their adaptive timesteps (leading to slow convergence in scipy and immediate abort in Julia, especially for the 5th order solvers since my 4th/5th derivatives are probably even less well-defined).

I think you’re right. I think you would need at least n+1 degree smoothing to use an order n solver.

How do you make sure the nth order derivative doesn’t blow up? What smoothing/interpolation algorithm ensures that at least the 3rd (or 5th) order derivative will be smooth?

I noticed that if I use s=0 and k=5 for the UnivariateSpline in scipy then it will go through all the input points (an interpolating spline), whereas if I set s=None, it will adjust the smoothing so that the 2nd and 3rd derivatives are more continuous (as @ChrisRackauckas suggested). For certain input time series that are more noisy than others I also noticed I needed to Gaussian smooth with higher sigma than ~10, more like sigma~50, so that the fifth order smoothing spline has continuous 3rd-5th derivatives.

A PI adaptive algorithm would likely be more sensitive to an actually underlying discontinuous higher derivative since it attempts to smooth such changes, so if it truly is discontinuous that could effect the integrator. Or another thing could just be default tolerances.

1 Like

Sorry but what is a PI adaptive algorithm? Is there an existing implementation in python/julia?

I played around with smoothing splines some more today. And it’s so strange that even a 5th order smoothing spline leads to a continuous 3rd order derivative but the 4th and 5th order derivatives show wild oscillations (see attached). And RK23 is still taking very small timesteps despite the 3rd order derivative being smooth.

All of the Julia methods default to PI adaptive time stepping.