使用scipy.integrate.ode自适应步长

Question

(简要)文档scipy.integrate.ode说两个方法(dopri5和dop853)具有逐步控制和密集输出.看一下示例和代码本身,我只能看到从集成器获取输出的一种非常简单的方法.也就是说,看起来你只是通过一些固定的dt向前推进积分器,得到当时的函数值,然后重复.

我的问题有很多可变的时间尺度,所以我想在需要评估的任何时间步长得到值,以达到所需的公差.也就是说,在早期,事情正在缓慢变化,因此输出时间步长可能很大.但随着事情变得有趣,输出时间步长必须更小.我实际上并不想要等间隔的密集输出,我只想要自适应函数使用的时间步长.

编辑:密集输出

一个相关的概念(几乎相反)是"密集输出",其中所采取的步骤与步进器所采取的步骤一样大,但函数的值被插值(通常具有与步进器精度相当的精度)你要.fortran底层scipy.integrate.ode显然能够做到这一点,但ode没有界面. odeint另一方面,它基于不同的代码,并且显然确实执行密集输出.(您可以在每次调用右侧时输出以查看何时发生,并看到它与输出时间无关.)

所以我仍然可以利用自适应性,只要我能够提前确定我想要的输出时间步长.不幸的是,对于我最喜欢的系统,我甚至不知道大概的时间尺度是什么时间函数,直到我运行集成.因此,我必须将采用一个积分器步骤的想法与这种密集输出的概念结合起来.

编辑2:再次密集输出

显然,scipy 1.0.0通过新接口引入了对密集输出的支持.特别是,他们建议远离scipy.integrate.odeint和朝向scipy.integrate.solve_ivp,作为关键字dense_output.如果设置为True,则返回的对象具有sol可以使用数组调用的属性,然后在这些时间返回集成函数值.这仍然没有解决这个问题的问题,但在许多情况下它很有用.

Answer 1

我一直在考虑这个尝试得到相同的结果.事实证明,你可以通过在ode实例化中设置nsteps = 1来使用hack来获得逐步结果.它将在每一步生成一个UserWarning(这可以被捕获和抑制).

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
import warnings


def logistic(t, y, r):
    return r * y * (1.0 - y)

r = .01
t0 = 0
y0 = 1e-5
t1 = 5000.0

#backend = 'vode'
backend = 'dopri5'
#backend = 'dop853'

solver = ode(logistic).set_integrator(backend, nsteps=1)
solver.set_initial_value(y0, t0).set_f_params(r)
# suppress Fortran-printed warning
solver._integrator.iwork[2] = -1

sol = []
warnings.filterwarnings("ignore", category=UserWarning)
while solver.t < t1:
    solver.integrate(t1, step=True)
    sol.append([solver.t, solver.y])
warnings.resetwarnings()
sol = np.array(sol)

plt.plot(sol[:,0], sol[:,1], 'b.-')
plt.show()

结果:

Answer 2

自SciPy 0.13.0起,

dopri现在可以通过solout回调函数访问ODE求解器族的中间结果.

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt


def logistic(t, y, r):
    return r * y * (1.0 - y)

r = .01
t0 = 0
y0 = 1e-5
t1 = 5000.0

backend = 'dopri5'
# backend = 'dop853'
solver = ode(logistic).set_integrator(backend)

sol = []
def solout(t, y):
    sol.append([t, *y])
solver.set_solout(solout)
solver.set_initial_value(y0, t0).set_f_params(r)
solver.integrate(t1)

sol = np.array(sol)

plt.plot(sol[:,0], sol[:,1], 'b.-')
plt.show()

结果:

结果似乎与Tim D's略有不同,尽管他们都使用相同的后端.我怀疑这与FSAL属性有关dopri5.在Tim的方法中,我认为第七阶段的结果k7被丢弃,因此重新计算k1.

注意:如果在设置初始值后设置set_solout,则set_solout无效.它固定为SciPy 0.17.0.

Answer 3

该integrate方法接受一个布尔参数step,该参数告诉该方法返回单个内部步骤.然而,似乎'dopri5'和'dop853'解算器不支持它.

以下代码显示了在使用"vode"解算器时如何获取解算器所采取的内部步骤:

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt


def logistic(t, y, r):
    return r * y * (1.0 - y)

r = .01

t0 = 0
y0 = 1e-5
t1 = 5000.0

backend = 'vode'
#backend = 'dopri5'
#backend = 'dop853'
solver = ode(logistic).set_integrator(backend)
solver.set_initial_value(y0, t0).set_f_params(r)

sol = []
while solver.successful() and solver.t < t1:
    solver.integrate(t1, step=True)
    sol.append([solver.t, solver.y])

sol = np.array(sol)

plt.plot(sol[:,0], sol[:,1], 'b.-')
plt.show()

结果: