通过 scipy.integrate.ode 使用自适应步长
Posted
技术标签:
【中文标题】通过 scipy.integrate.ode 使用自适应步长【英文标题】:Using adaptive step sizes with scipy.integrate.ode 【发布时间】:2012-10-07 05:50:12 【问题描述】:scipy.integrate.ode
的(简要)文档说两种方法(dopri5
和 dop853
)具有步长控制和密集输出。查看示例和代码本身,我只能看到一种从集成器获取输出的非常简单的方法。即,看起来您只是将积分器向前移动了某个固定的 dt,获取当时的函数值,然后重复。
我的问题有相当多变的时间尺度,所以我想在需要评估的任何时间步骤获取值以达到所需的容差。也就是说,在早期,事情正在缓慢变化,因此输出时间步长可能很大。但随着事情变得有趣,输出时间步长必须更小。我实际上并不想要等间隔的密集输出,我只想要自适应函数使用的时间步长。
编辑:密集输出
一个相关的概念(几乎相反)是“密集输出”,即所采取的步数与步进器一样大,但函数的值是插值的(通常精度与步进器的精度相当) 到任何你想要的。底层scipy.integrate.ode
的fortran 显然可以做到这一点,但ode
没有接口。另一方面,odeint
是基于不同的代码,并且显然会进行密集输出。 (您可以在每次调用右侧时输出以查看何时发生,并查看它与输出时间无关。)
所以我仍然可以利用自适应性,只要我可以提前决定我想要的输出时间步长。不幸的是,对于我最喜欢的系统,我什至不知道作为时间函数的大致时间尺度是什么,直到我运行集成。因此,我必须将采取积分器步骤的想法与密集输出的概念结合起来。
编辑 2:再次密集输出
显然,scipy 1.0.0 通过一个新接口引入了对密集输出的支持。特别是,他们建议从scipy.integrate.odeint
转向scipy.integrate.solve_ivp
,将其作为关键字dense_output
。如果设置为True
,则返回的对象具有属性sol
,您可以使用时间数组调用该属性,然后返回这些时间的集成函数值。这仍然不能解决这个问题的问题,但它在很多情况下很有用。
【问题讨论】:
【参考方案1】:自从SciPy 0.13.0,
dopri
系列 ODE 求解器的中间结果可以 现在可以通过solout
回调函数访问。
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 的略有不同,尽管他们都使用相同的后端。我怀疑这与dopri5
的FSAL 属性有关。在Tim的方法中,我认为第七阶段的结果k7被丢弃了,所以重新计算k1。
注意:有一个已知的bug with set_solout not working if you set it after setting initial values。从SciPy 0.17.0 开始,它已修复。
【讨论】:
请记住,现在(0.16 版)您必须调用set_solout()
之前 set_initial_value()
否则将不会调用solout。
@syockit 我已经尝试过上述方法(dopri5 和 vode),我发现结果不是确定性的。这意味着对于我的一组僵硬的 ODE,算法有时会成功完成,有时会失败(给定相同的输入)。这对我不起作用。你知道如何强制它成为确定性的,或者我可以在哪里寻找确定性的解决方案?谢谢。
@theJollySin 抱歉回复晚了。不幸的是,稳定算法的选择很大程度上取决于问题的领域,所以我不能给你一个好的一刀切的建议。 scipy.integrate.ode 的当前文档告诉您在处理棘手的问题时使用哪一个或做什么。【参考方案2】:
仅供参考,虽然答案已经被接受,但我应该指出历史记录,PyDSTool 原生支持沿计算轨迹的任何地方的密集输出和任意采样。这还包括求解器内部使用的所有自适应确定的时间步长的记录。它与dopri853 and radau5 接口,并自动生成与它们接口所需的 C 代码,而不是依赖于(慢得多的)python 函数回调来进行右侧定义。据我所知,任何其他专注于 python 的求解器都没有原生或有效地提供这些功能。
【讨论】:
【参考方案3】:这里还有一个选项也应该适用于dopri5
和dop853
。基本上,求解器会根据需要经常调用logistic()
函数来计算中间值,以便我们存储结果:
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
sol = []
def logistic(t, y, r):
sol.append([t, y])
return r * y * (1.0 - y)
r = .01
t0 = 0
y0 = 1e-5
t1 = 5000.0
# Maximum number of steps that the integrator is allowed
# to do along the whole interval [t0, t1].
N = 10000
#backend = 'vode'
backend = 'dopri5'
#backend = 'dop853'
solver = ode(logistic).set_integrator(backend, nsteps=N)
solver.set_initial_value(y0, t0).set_f_params(r)
# Single call to solver.integrate()
solver.integrate(t1)
sol = np.array(sol)
plt.plot(sol[:,0], sol[:,1], 'b.-')
plt.show()
【讨论】:
这里涉及三个时间步长。最大的是输出时间步长;ode
默认输出什么。然后我们有了积分器的自适应时间步长,这就是我想要输出的地方。最小的是 Runge-Kutta 风格的积分器(例如,dopri5
、dop853
)使用的中间时间步长。也就是说,dopri5
实际上会在组成单个时间步的中间步骤期间多次调用logistic
。我担心的是我相信这些中间步骤的准确性较低。在许多情况下,y
值实际上只是一阶准确的。
“我相信那些中间步骤的精度较低” 糟糕,我认为 logistic
会被积分器的自适应时间步长调用。如果没有,我的回答当然没有帮助。但你确定吗?
这个想法是为正确的增量建立更好的近似值;结合多次调用的结果可以让您取消一些错误术语。事实上,现在我看了一下,***的文章有一节解释"The Runge-Kutta method",其中显示该函数在步中点处用不同的值计算了两次。因此,在某些情况下,您的 sol
实际上会为相同的 t
值存储两个不同的 y
值,因为第一个值不太准确。【参考方案4】:
我一直在研究这个以尝试获得相同的结果。事实证明,您可以通过在 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()
结果:
【讨论】:
这似乎可以完成这项工作,好吧。现在,如果我没有花费数周时间将我的代码移植到 GSL...:) 事实证明,您可以使用这个小技巧solver._integrator.iwork[2] = -1
抑制 Fortran 发出的警告(我将编辑上面的代码以显示这一点)。这会设置一个通过 Fortran 接口传递的标志,该标志禁止打印到标准输出。【参考方案5】:
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()
结果:
【讨论】:
是的,我担心是这样。我希望有一种简单的方法来扩展dopri5
和dop853
,但我的耐心到了fortran 结束,所以我想我会重新实现一个时间步进器。不过,似乎很遗憾,python 没有强大、高效和灵活的集成器......
我在尝试转换使用 ode45 的 MATLAB 脚本时需要同样的功能。我已经在 Scipy.org Ticket#1820 提交了一张票。
@Mike:好吧,即使使用dopri5
和dop853
,您也可以始终将值存储在logistic
函数内部,只需调用integrate
方法即可。 (我将把它作为替代答案发布。)
另一方面:我喜欢step=True
功能isn't documented at all。
@Mike 请参阅下面关于 PyDSTool 使用 dopri 853 对密集输出和内部步骤的本机支持的说明。以上是关于通过 scipy.integrate.ode 使用自适应步长的主要内容,如果未能解决你的问题,请参考以下文章
scipy.integrate.odeint 和 scipy.integrate.ode 有啥区别?
scipy.integrate.ode 与两个耦合的 ODE?