python - 具有大量数据的 SciPyIntegrate.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 在 float 中没有区别点数。这就是您收到错误的原因。

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

一般建议是，当使用通用求解器包重新调整问题时，使求解器看到的状态向量和步长处于可观的范围内，介于 1e-3 和 1e6。您可以在理论方面通过选择适当的单位或在将状态向量转换为模型变量并转换回来时在 ODE 函数中进行重新调整来实现此目的。