ODE 求解器的协程

计算科学 软件
2021-12-19 18:52:33

是否有任何使用协程并产生结果而不是函数并返回的 ODE 求解器包?

简而言之,编程语言中的子程序进行一些计算,将结果返回给调用它的人,然后就完成了。协程进行一些计算,为调用它的人产生结果,然后可以在它上次产生值的地方恢复。在 Python 中,协程用于实现生成器。

对于一些涉及求解 ODE 系统的科学应用,关于积分多长时间的逻辑可能非常复杂。例如,如果您正在编写一些代码来创建矢量场的流图,那么如果 (1) 它还不够长并且 (2) 它不太接近您已经计算的所有其他积分曲线。第二个停止条件要复杂得多,涉及某种空间数据结构,如 KD 树。大多数 ODE 求解例程(例如 SciPy 的solve_ivp函数)为您提供的选项数量非常有限,用于何时停止积分,您必须从一开始就决定您的标准是什么。

相比之下,很容易编写一个自适应积分方案,它会产生积分曲线上的下一个点、该点的向量场值以及它使用的时间步长。ODE 集成例程不关心调用者希望它运行多长时间。它可能是屏幕保护程序的一部分,并且会一直运行,直到用户摆动鼠标。

3个回答

您实际上并不需要协同程序——您可以通过在用户代码中使用常规回调来实现相同的目的。例如,在 C 编程语言中,ODE 求解器例程将有一个额外的参数来表示指向函数的指针。用户代码将使用指向用户函数的指针调用积分器,然后 ODE 积分器在每个时间步结束时调用该用户函数(将当前解决方案作为参数发送)。如果用户函数返回带有true值的 ODE 积分器,则积分将继续。如果用户函数返回false,则积分器终止,释放内存,并返回给 ODE 积分器的调用者。

在诸如 C++ 之类的语言中,您有更多选择:重载基类,或者signal在用户代码可以附加到的 ODE 集成器对象中使用 boost-style s。重点是,所有这些都可以使用传统语言轻松实现。与您的示例相比,它只是颠倒了谁驾驶谁的逻辑。

您展示的示例很有趣,因为您需要大量的积分器和一些顶级逻辑来决定这些积分器中的哪一个下一步要走。通过您提到的协同例程最容易做到这一点,但使用回调也不是那么困难。如果每个集成器都在自己的线程上运行,那么使用回调实现可能最容易,并且回调访问一些共享状态,只要其中一个线程调用回调,它们就会修改这些共享状态。这似乎与协同例程方法完全不同,但值得记住的是,协同例程通常也使用每个协同例程一个堆栈来实现——实际上与单独的线程并没有太大区别。

有些集成器向用户公开了更细粒度的控制,例如scipy.integrate.ode

from scipy.integrate import ode

def f(t,y):
    return y[0]

initial = [1]
ODE = ode(f)
ODE.set_integrator("dopri5")
ODE.set_initial_value(initial)

for time in range(1,10):
    print(ODE.integrate(time))

代替协程,这提供了作为ODE类 ( ode) 的实例 () 的集成器,该类配备了一个方法,该方法执行与未来 () 的小集成integrate并返回结果。

我不确定您是否从这里使用协程中获得任何好处,但如果您真的愿意,可以直接在此周围包装一个协程:

from scipy.integrate import ode

def f(t,y):
    return y[0]

def integrator(f,initial):
    ODE = ode(f)
    ODE.set_integrator("dopri5")
    ODE.set_initial_value(initial)

    target_time = 0
    while True:
        target_time = yield ODE.integrate(target_time)

my_coroutine = integrator(f,[1.0])
next(my_coroutine)
for time in range(1,10):
    print(my_coroutine.send(time))

不完全是一个,但Revised6 Report on the Algorithmic Language Scheme末尾的示例实现了一个经典的 Runge-Kutta 积分器作为返回“无限状态流”的函数。

这是在几行代码中完成的,并且使用它的协程、生成器、itertools 等很容易移植到例如 Python。

我认为这是一个令人信服的想法。自从在上个千年末的 Scheme 报告的早期修订中第一次看到它以来,我一直使用基于此的东西来集成 ODE 以及隐式 ODE 系统(例如,由没有“质量”的有限元离散化产生集总”,因此具有非对角“质量”矩阵)和微分代数方程,这些方程不能由上面提到的 SciPy 等标准 ODE 包处理。到目前为止,这都是内部实现,太临时而无法发布。

在过去的二十年里,我有好几次寻找一个包(特别是 Python 包),它实现了类似 R6RS 的方法来为动态系统(特别是比 ODE 更通用)生成轨迹,但还没有找到。GarlicSim非常好,但在我找到它时它已经被遗弃了。最近,我认为基于StreamzScikit FiniteDiff的流接口看起来很有前途,特别是对于更复杂的(例如耦合或分层)问题。