具有大量数字的 SciPy 集成代码答案

【问题标题】：SciPy integrate.ode with large numbers具有大量数字的 SciPy 集成代码
【发布时间】：2018-09-18 00:43:09
【问题描述】：

我正在尝试使用 SciPy 的集成.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 的元素，这是不应该的。

我该如何解决这个问题？

【问题讨论】：

标签： python scipy ode

【解决方案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 函数中重新调整来做到这一点。

【讨论】：