带有大数的SciPy integrate.ode

Question

我正在尝试使用SciPy的integrate.ode来解决一个ODE，但是失败了。我怀疑这与网格的大小有关，所以我尝试了以下两个测试：

from scipy.integrate import ode

def bla(t,x):
    return 0
M0=1
t0,y0 = 0.5*M0,0.5*M0
r = ode(bla).set_integrator('dopri5')
r.set_initial_value(y0,t0)
t1 = M0
dt =0.01*M0
bli = np.array([])
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    bli = np.append(bli,r.y)

数组bli包含50个值为0.5的元素。

from scipy.integrate import ode

def bla(t,x):
    return 0
M0=1e13
t0,y0 = 0.5*M0,0.5*M0
r = ode(bla).set_integrator('dopri5')
r.set_initial_value(y0,t0)
t1 = M0
dt =0.01*M0
bli = np.array([])
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    bli = np.append(bli,r.y)

本例中的数组bli包含一个值为5e12和r.successful() = False的元素，因为它不应该是这样的。

我该如何解决这个问题呢？

Answer 1

如果在bla中插入print t,x语句，您可以观察到，作为初始化的一部分，求解器使用硬编码的步长h=1e-5执行积分步骤，以获得第一个“实际”步骤的最佳步长。计算的第二个值位于偏移量0.1*h=1e-6处。在t=5e12中，t与t+0.1*h或t+h在浮点数上没有区别。这就是你得到这个错误的原因。

除非更改ode类的代码，或者使用更动态地处理初始化的不同版本或包，否则无法解决此问题。

一般建议是，在使用常规解算器软件包重新缩放问题时，使解算器看到的状态向量和步长处于明显的范围内，介于1e-3和1e6之间。您可以通过在理论上选择适当的单位或在将状态向量转换为模型变量时在ODE函数中进行重新缩放来实现这一点。