1

我有一段代码增加了所谓的Lorenz95 模型（由 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 并不慢，但不幸的是，我的代码需要对它进行非常大量的调用。有没有办法我可以编写时间步长以使其运行得更快？

非常感谢。

4

1 回答 1

3

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

  • 在算法方面，您使用的是Euler 方法，这是一种一阶方法（因此全局误差与步长成正比）并且具有较小的稳定区域。也就是说，它不是很有效。

  • 另一方面，如果你使用标准的 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)，你的可以不同）非常混乱。初始条件中的小扰动（仅在他的一个坐标中！）产生非常不同的解决方案： Lorenz-95动力系统

但是，它真的比欧拉时间步长更快吗？因为dt = 0.01它似乎快了近三倍，但除了一开始之外，解决方案并不匹配。 欧拉 vs odeint

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

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

于 2013-05-08T17:23:40.793 回答