Solving log-transformed ODE system without overflow error

Question

I have a system of ODEs where my state variables and independent variable span many orders of magnitude (initial values are around 0 at t=0 and are expected to become about 10¹⁰ by t=10¹⁷). I also want to ensure that my state variables remain positive.

According to this Stack Overflow post, one way to enforce positivity is to log-transform the ODEs to solve for the evolution of the logarithm of a variable instead of the variable itself. However when I try this with my ODEs, I get an overflow error probably because of the huge dynamic range / orders of magnitude of my state variables and time variable. Am I doing something wrong or is log-transform just not applicable in my case?

Here is a minimal working example that is successfully solved by scipy.integrate.solve_ivp:

import numpy as np
from scipy.interpolate import interp1d
from scipy.integrate import solve_ivp

# initialize times at which we are given certain input quantities/parameters
# this is seconds corresponding to the age of the universe in billions of years
times = np.linspace(0.1,10,500) * 3.15e16 

# assume we are given the amount of new mass flowing into the system in units of g/sec
# for this toy example we will assume a log-normal distribution and then interpolate it for our integrator function
mdot_grow_array = np.random.lognormal(mean=0,sigma=1,size=len(times))*1.989e33 / 3.15e7
interp_grow = interp1d(times,mdot_grow_array,kind='cubic')

# assume there is also a conversion efficiency for some fraction of mass to be converted to another form
# for this example we'll assume the fractions are drawn from a uniform random distribution and again interpolate
mdot_convert_array = np.random.uniform(0,0.1,len(times)) / 3.15e16 # fraction of M1 per second converted to M2
interp_convert = interp1d(times,mdot_convert_array,kind='cubic')

# set up our integrator function
def integrator(t,y):
    print('Working on t=',t/3.15e16) # to check status of integration in billions of years
    
    # unpack state variables
    M1, M2 = y
    
    # get the interpolated value of new mass flowing in at this time
    mdot_grow_now = interp_grow(t)
    mdot_convert_now = interp_convert(t)    
    
    # assume some fraction of the mass gets converted to another form 
    mdot_convert = mdot_convert_now * M1 
    
    # return the derivatives
    M1dot = mdot_grow_now - mdot_convert
    M2dot = mdot_convert
    
    return M1dot, M2dot
    
# set up initial conditions and run solve_ivp for the whole time range 
# should start with M1=M2=0 initially but then solve_ivp does not work at all, so just use [1,1] instead
initial_conditions = [1.0,1.0] 

# note how the integrator gets stuck at very small timesteps early on
sol = solve_ivp(integrator,(times[0],times[-1]),initial_conditions,dense_output=True,method='RK23')

And here is the same example but now log-transformed following the Stack Overflow post referenced above (since dlogx/dt = 1/x * dx/dt, we simply replace the LHS with x*dlogx/dt and divide both sides by x to isolate dlogx/dt on the LHS; and we make sure to use np.exp() on the state variables – now logx instead of x – within the integrator function):

import numpy as np
from scipy.interpolate import interp1d
from scipy.integrate import solve_ivp

# initialize times at which we are given certain input quantities/parameters
# this is seconds corresponding to the age of the universe in billions of years
times = np.linspace(0.1,10,500) * 3.15e16 

# assume we are given the amount of new mass flowing into the system in units of g/sec
# for this toy example we will assume a log-normal distribution and then interpolate it for our integrator function
mdot_grow_array = np.random.lognormal(mean=0,sigma=1,size=len(times))*1.989e33 / 3.15e7
interp_grow = interp1d(times,mdot_grow_array,kind='cubic')

# assume there is also a conversion efficiency for some fraction of mass to be converted to another form
# for this example we'll assume the fractions are drawn from a uniform random distribution and again interpolate
mdot_convert_array = np.random.uniform(0,0.1,len(times)) / 3.15e16 # fraction of M1 per second converted to M2
interp_convert = interp1d(times,mdot_convert_array,kind='cubic')

# set up our integrator function
def integrator(t,logy):
    print('Working on t=',t/3.15e16) # to check status of integration in billions of years
    
    # unpack state variables
    M1, M2 = np.exp(logy)
    
    # get the interpolated value of new mass flowing in at this time
    mdot_grow_now = interp_grow(t)
    mdot_convert_now = interp_convert(t)    
    
    # assume some fraction of the mass gets converted to another form 
    mdot_convert = mdot_convert_now * M1 
    
    # return the derivatives
    M1dot = (mdot_grow_now - mdot_convert) / M1 
    M2dot = (mdot_convert) / M2
    
    return M1dot, M2dot
    
# set up initial conditions and run solve_ivp for the whole time range 
# should start with M1=M2=0 initially but then solve_ivp does not work at all, so just use [1,1] instead
initial_conditions = [1.0,1.0] 

# note how the integrator gets stuck at very small timesteps early on
sol = solve_ivp(integrator,(times[0],times[-1]),initial_conditions,dense_output=True,method='RK23')
    

Answer 1

[…] is log-transform just not applicable in my case?

I don’t know where your transform went wrong, but it will certainly not achieve what you think it does. Log-transforming as a means to avoid negative values makes sense and works if and only if the following two conditions hold:

1. If the value of a dynamical variable approaches zero (from above), its derivative also approaches zero (from above) in your model.
2. Due to numerical noise, your derivative may turn negative though it actually isn’t.

Conversely, it is not necessary or doesn’t work in the following cases:

• If Condition 1 fails because your derivative never approaches zero in your model, but is strictly positive, you have no problem to begin with, as your derivative should not become negative in any reasonable implementation of your model. (You might make it happen by implementing some spectacular numerical annihilation, but that’s quite a difficult feat to achieve and not what I would consider a reasonable implementation.)

• If Condition 1 fails because your derivative becomes truly negative in your model, logarithms won’t save you, because the dynamics wants to push the derivative below zero and the logarithms cannot represent this. You usually get an overflow error due to the logarithms becoming extremely negative or the adaptive integration fails.

• Even if Condition 1 applies, Condition 2 can usually be handled by avoiding numerical annihilations and similar when implementing your model.

Unless I am mistaken, your model falls into the first category. If M1 goes to zero, mdot_convert goes towards zero and thus M1dot = mdot_grow_now - mdot_convert is strictly positive, because mdot_grow_now is. M2dot is strictly positive anyway. Thus, you gain nothing from log-transforming. In fact, in the vast majority of cases, your dynamical variables will quickly increase.

With all that being said, some things you might want to look into are:

• Normalising your variables to be in the order of magnitude of 1.
• Stochastic differential equations.