python - 如何在Python中优化时间步长算法？

Question

我有一段代码，它增加了所谓的Lorenz95 model的时间步长。 (由埃德·洛伦茨于 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)，你的可能不同)非常困惑。初始条件中的小扰动(仅在一个坐标中!)会产生非常不同的解决方案:

但是，它真的比欧拉时间步长更快吗？对于dt = 0.01，它看起来几乎快了三倍，但除了一开始之外，解决方案并不匹配。

如果减少dt，欧拉方法提供的解决方案会变得越来越类似于odeint解决方案，但需要更长的时间。请注意较小的 dt，后来的欧拉解决方案如何失去对 odeint 解决方案的跟踪。最精确的欧拉解计算 t=6 之前的解所需的时间比计算 t=10 之前的 odeint 的时间长 600 倍。查看完整脚本here 。

最后，这个系统非常不稳定，我猜甚至 odeint 解决方案在所有绘制的时间内都不是准确的。