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

问题描述

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

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.

一个相关的概念(几乎相反)是密集输出"，即所采取的步长与步进器要执行的步长一样大，但是函数的值是插值的(通常其精度可与步进器的精度相提并论) )到您想要的任何内容.底层scipy.integrate.ode的fortran显然可以做到这一点，但是ode没有该接口.另一方面，odeint基于不同的代码，并且显然可以进行密集输出. (您可以在每次调用右侧时输出，以查看发生这种情况的时间，并查看与输出时间无关.)

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.

显然，scipy 1.0.0引入了对通过新接口进行密集输出的支持.特别是，他们建议从scipy.integrate.odeint转到 scipy.integrate.solve_ivp ，作为关键字dense_output.如果设置为True，则返回的对象具有属性sol，您可以使用一个时间数组调用该属性，然后在该时间返回集成函数值.仍然不能解决这个问题的问题，但是在很多情况下它是有用的.

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.

推荐答案

由于 SciPy 0.13.0 ，

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

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

结果:

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

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.

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

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屋！