解决日志转换的ODE系统没有溢出错误

Question

我有一个常微分方程系统，其中我的状态变量和自变量跨越许多数量级(初始值在t=0时约为0，预计约为10⁰通过t=10⁷).我还想确保我的状态变量保持为正。

根据这篇Stack Overflow的帖子，一种增强正性的方法是对数变换ODE，以解决变量对数的演变，而不是变量本身。然而，当我尝试使用ODE时，我会出现溢出错误，这可能是因为我的状态变量和时间变量的巨大动态范围/数量级。我是做错了什么，还是日志转换不适用于我的情况？

以下是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的数量级
• 随机微分方程