求解无溢出错误的对数变换 ODE 系统

Question

我有一个 ODE 系统，其中状态变量和自变量跨越多个数量级（在 t=0 时初始值约为 0，预计到 t=1017 时将变为大约 1010）。我还想确保我的状态变量保持正值。

根据此 Stack Overflow 帖子 ，增强正性的一种方法是对常微分方程进行对数变换来求解变量的对数演化，而不是变量本身。然而，当我用常微分方程尝试此操作时，我收到溢出错误，可能是因为状态变量和时间变量的动态范围/数量级很大。我做错了什么或者对数转换不适用于我的情况吗？

这是一个由 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')

这是相同的示例，但现在按照上面引用的 Stack Overflow 帖子进行了日志转换（因为 dlogx/dt = 1/x * dx/dt，我们只需将 LHS 替换为 x*dlogx/dt 并将两边除以 x 即可隔离LHS 上的 dlogx/dt；并且我们确保在状态变量上使用 np.exp() – 现在使用 logx 而不是x ——在积分器函数内）：

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

[…]对数变换不适用于我的情况吗？

我不知道你的变换哪里出了问题，但它肯定达不到你想象的效果。当且仅当满足以下两个条件时，对数转换作为避免负值的一种方法才有意义并且有效：

1. 如果动态变量的值接近零（从上面），则其模型中的导数也接近零（从上面） 。
2. 由于数值噪声，您的导数可能会变为负值，但实际上并非如此。

相反，在以下情况下没有必要或不起作用：

• 如果条件 1 失败，因为模型中的导数永远不会接近零，而是严格为正数，那么一开始就没有问题，因为您的导数在模型的任何合理实施中不应变得消极。 （你可能会通过实施一些壮观的数字歼灭来实现这一目标，但这是一项相当困难的壮举，而且不是我认为合理的实施。）

• 如果条件 1 失败，因为您的导数在模型中变得真正为负，对数不会拯救你，因为动力学想要将导数推到零以下，而对数不能代表这一点。由于对数变得非常负或自适应积分失败，您通常会收到溢出错误。

• 即使条件 1 适用，条件 2 通常也可以通过在实现模型时避免数值歼灭和类似情况来处理。

除非我弄错了，否则你的模型属于第一类。如果 M1 趋于零，mdot_convert 趋于零，因此 M1dot = mdot_grow_now - mdot_convert 严格为正数，因为 mdot_grow_now > 是。无论如何，M2dot 是严格正数。因此，您不会从对数转换中获得任何好处。事实上，在绝大多数情况下，你的动态变量会迅速增加。

话虽如此，您可能需要研究的一些事情是：

• 将变量标准化为 1 的数量级。
• 随机微分方程。