如何在python中优化时间步进算法？

Question

我有一段代码，它为所谓的Lorenz95 model (由Ed Lorenz在1995年发明)增加一个时间步长。它通常作为40个变量的模型实现，并显示混沌行为。我已经对算法的时间步长进行了如下编码：

class Lorenz:
    '''Lorenz-95 equation'''

    global F, dt, SIZE
    F = 8
    dt = 0.01
    SIZE = 40

    def __init__(self):
        self.x = [random.random() for i in range(SIZE)]

    def euler(self):
        '''Euler time stepping'''
        newvals = [0]*SIZE
        for i in range(SIZE-1):
            newvals[i] = self.x[i] + dt * (self.x[i-1] * (self.x[i+1] - self.x[i-2]) - self.x[i] + F)
        newvals[SIZE-1] = self.x[SIZE-1] + dt * (self.x[SIZE-2] * (self.x[0] - self.x[SIZE-3]) - self.x[SIZE-1] + F)
        self.x = newvals

这个函数euler并不慢，但不幸的是，我的代码需要对它进行大量的调用。有没有一种方法可以编码时间步进，让它运行得更快？

非常感谢。

Answer 1

至少有两种可能的优化:以更智能的方式工作(算法改进)和更快地工作。

在算法方面，你使用的是Euler method，这是一种一阶方法(所以全局误差与步长成正比)，并且具有较小的稳定区域。也就是说，效率不是很高，在另一端，如果你使用标准的CPython实现，这种代码会很慢。要解决这个问题，您可以简单地尝试在PyPy下运行它。它的即时编译器可以使数值代码的运行速度提高100倍。您还可以编写自定义的C或Cython扩展。

但是有一个更好的方法。求解常微分方程组非常常见，因此Python中的核心科学库之一scipy包装了快速的、经过战斗测试的Fortran库来解决这些问题。通过使用scipy，您可以获得算法改进(因为积分器将具有更高的阶数)和快速实现。

对于一组扰动的初始条件，求解Lorenz 95模型如下所示：

import numpy as np


def lorenz95(x, t):
    return np.roll(x, 1) * (np.roll(x, -1) - np.roll(x, 2)) - x + F

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    from scipy.integrate import odeint
    SIZE = 40
    F = 8
    t = np.linspace(0, 10, 1001)
    x0 = np.random.random(SIZE)
    for perturbation in 0.1 * np.random.randn(5):
        x0i = x0.copy()
        x0i[0] += perturbation
        x = odeint(lorenz95, x0i, t)
        plt.plot(t, x[:, 0])
    plt.show()

并且输出(设置np.random.seed(7)，您的可以不同)非常混乱。初始条件中的小扰动(仅在一个he坐标中！)产生截然不同的解决方案：

​

但是，它真的比Euler时间步长更快吗？对于dt = 0.01来说，它似乎快了近三倍，但除了最开始的时候，其他的解决方案都不匹配。

​

如果减少了dt，Euler方法提供的解决方案将变得越来越类似于odeint解决方案，但需要的时间要长得多。请注意，较小的dt，较晚的Euler解决方案如何放松对odeint解决方案的跟踪。最精确的Euler解决方案计算到t=6的时间是odeint到t=10的600倍。请参阅完整的脚本here。

​

最后，这个系统是如此不稳定，我猜即使是odeint解决方案也不是在所有绘制的时间上都是准确的。