求解无溢出错误的对数转换 ODE 系统
Solving log-transformed ODE system without overflow error
我有一个 ODE 系统，其中我的状态变量和自变量跨越多个数量级（初始值在 t=0 时大约为 0，预计将变为大约 10¹⁰通过 t=10¹⁷)。我还想确保我的状态变量保持正数。
，强制正性的一种方法是对 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 post 进行对数转换（因为 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')
[…] is log-transform just not applicable in my case?
不知道你的transform哪里错了，但肯定达不到你想的那样。当且仅当满足以下两个条件时，Log-transforming 作为一种避免负值的方法才有意义并且有效：
- 如果动态变量的值接近零（从上方），则其导数在您的模型中也接近零（从上方）。
- 由于数值噪声，您的导数可能会变为负数，但实际上并非如此。
反之，在下列情况下则不需要或不起作用：
如果条件 1 失败是因为您的导数在您的模型中永远不会接近零，而是严格为正，那么您一开始就没有问题，因为您的导数在您的任何合理实施中都不应该变为负数模型。 （你可能会通过实施一些壮观的数值湮灭来实现它，但这是一项很难实现的壮举，而不是我认为合理的实施。）
如果条件 1 失败，因为你的导数在你的模型中变成了真正的负数，对数将救不了你，因为动力学想要将导数推到零以下，而对数不能表示这一点。由于对数变得非常负或自适应积分失败，您通常会遇到溢出错误。
即使条件 1 适用，条件 2 通常可以通过在实现模型时避免数值湮灭和类似操作来处理。
如果我没有记错的话，你的模型属于第一类。如果 M1 趋于零，mdot_convert 趋于零，因此 M1dot = mdot_grow_now - mdot_convert 是严格正的，因为 mdot_grow_now 是。 M2dot 无论如何都是严格正的。因此，您从 log-transforming 中一无所获。事实上，在绝大多数情况下，你的动力变量会迅速增加。
综上所述，您可能想要了解的一些内容是：
. - 随机微分方程。
我有一个 ODE 系统，其中我的状态变量和自变量跨越多个数量级（初始值在 t=0 时大约为 0，预计将变为大约 10¹⁰通过 t=10¹⁷)。我还想确保我的状态变量保持正数。
这是一个最小的工作示例，已被 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 post 进行对数转换（因为 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')
[…] is log-transform just not applicable in my case?
不知道你的transform哪里错了，但肯定达不到你想的那样。当且仅当满足以下两个条件时，Log-transforming 作为一种避免负值的方法才有意义并且有效：
- 如果动态变量的值接近零（从上方），则其导数在您的模型中也接近零（从上方）。
- 由于数值噪声，您的导数可能会变为负数，但实际上并非如此。
反之，在下列情况下则不需要或不起作用：
如果条件 1 失败是因为您的导数在您的模型中永远不会接近零，而是严格为正，那么您一开始就没有问题，因为您的导数在您的任何合理实施中都不应该变为负数模型。 （你可能会通过实施一些壮观的数值湮灭来实现它，但这是一项很难实现的壮举，而不是我认为合理的实施。）
如果条件 1 失败，因为你的导数在你的模型中变成了真正的负数，对数将救不了你，因为动力学想要将导数推到零以下，而对数不能表示这一点。由于对数变得非常负或自适应积分失败，您通常会遇到溢出错误。
即使条件 1 适用，条件 2 通常可以通过在实现模型时避免数值湮灭和类似操作来处理。
如果我没有记错的话，你的模型属于第一类。如果 M1 趋于零，mdot_convert 趋于零，因此 M1dot = mdot_grow_now - mdot_convert 是严格正的，因为 mdot_grow_now 是。 M2dot 无论如何都是严格正的。因此，您从 log-transforming 中一无所获。事实上，在绝大多数情况下，你的动力变量会迅速增加。
综上所述，您可能想要了解的一些内容是：
- 随机微分方程。