在 scipy.integrate.ode 中使用自适应步长 [英] Using adaptive step sizes with scipy.integrate.ode

问题描述

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

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

密集输出

一个相关的概念(几乎相反)是密集输出"，即所采取的步骤与步进器所采取的一样大，但函数的值是内插的(通常精度与步进器的精度相当) 到任何你想要的.Fortran 底层 scipy.integrate.ode 显然能够做到这一点，但 ode 没有接口.另一方面，odeint 基于不同的代码，并且显然可以进行密集输出.(你可以在你的右手边每次被调用时输出，看看什么时候发生，看看它与输出时间无关.)

所以我仍然可以利用适应性，只要我可以提前决定我想要的输出时间步长.不幸的是，对于我最喜欢的系统，在我运行集成之前，我什至不知道作为时间函数的大概时间尺度是什么.所以我必须将采取一个积分器步骤的想法与密集输出的概念结合起来.

编辑 2:再次密集输出

显然，scipy 1.0.0 通过新界面引入了对密集输出的支持.特别是，他们建议从 scipy.integrate.odeint 转向

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

注意:有一个已知的 错误，如果您在之后设置 set_solout 则它不起作用设置初始值.自 SciPy 0.17.0 起已修复.

The (brief) documentation for scipy.integrate.ode says that two methods (dopri5 and dop853) have stepsize control and dense output. Looking at the examples and the code itself, I can only see a very simple way to get output from an integrator. Namely, it looks like you just step the integrator forward by some fixed dt, get the function value(s) at that time, and repeat.

My problem has pretty variable timescales, so I'd like to just get the values at whatever time steps it needs to evaluate to achieve the required tolerances. That is, early on, things are changing slowly, so the output time steps can be big. But as things get interesting, the output time steps have to be smaller. I don't actually want dense output at equal intervals, I just want the time steps the adaptive function uses.

EDIT: Dense output

A related notion (almost the opposite) is "dense output", whereby the steps taken are as large as the stepper cares to take, but the values of the function are interpolated (usually with accuracy comparable to the accuracy of the stepper) to whatever you want. The fortran underlying scipy.integrate.ode is apparently capable of this, but ode does not have the interface. odeint, on the other hand, is based on a different code, and does evidently do dense output. (You can output every time your right-hand-side is called to see when that happens, and see that it has nothing to do with the output times.)

So I could still take advantage of adaptivity, as long as I could decide on the output time steps I want ahead of time. Unfortunately, for my favorite system, I don't even know what the approximate timescales are as functions of time, until I run the integration. So I'll have to combine the idea of taking one integrator step with this notion of dense output.

EDIT 2: Dense output again

Apparently, scipy 1.0.0 introduced support for dense output through a new interface. In particular, they recommend moving away from scipy.integrate.odeint and towards scipy.integrate.solve_ivp, which as a keyword dense_output. If set to True, the returned object has an attribute sol that you can call with an array of times, which then returns the integrated functions values at those times. That still doesn't solve the problem for this question, but it is useful in many cases.

解决方案

Since SciPy 0.13.0,

The intermediate results from the dopri family of ODE solvers can now be accessed by a solout callback function.

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()

Result:

The result seems to be slightly different from Tim D's, although they both use the same backend. I suspect this having to do with FSAL property of dopri5. In Tim's approach, I think the result k7 from the seventh stage is discarded, so k1 is calculated afresh.

Note: There's a known bug with set_solout not working if you set it after setting initial values. It was fixed as of SciPy 0.17.0.

这篇关于在 scipy.integrate.ode 中使用自适应步长的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！